Numerical simulation method for oil reservoirs based on hybrid conductivity of meshless method and grid method

Abstract

The application discloses a numerical simulation method for oil reservoirs based on a hybrid conductivity of a meshless method and a grid method, and belongs to the technical field of numerical simulation of oil reservoirs. In view of the problem that the representative conductivity of the well circumference, the vicinity of a fault and a strong heterogeneous zone in a corner point grid model is insufficient, the application adopts a hybrid calculation scheme: the grid method conductivity of a reservoir simulator is used in a main body area of a reservoir; a point cloud is constructed by taking the center of the corner point grid as the center in a local high-precision demand area, the generalized finite difference method is used to reconstruct the in-layer connection, and the extended finite volume method is used to calculate the meshless conductivity; the local meshless connection and the grid method connection of the area boundary, interlayer and conventional area form a global connection table, and the numerical solution is completed by calling the reservoir simulator. The SNARK corner point grid example verifies the applicability of the method. The application has the advantages of calculation efficiency and local precision, does not need to develop a solver again, and is suitable for large-scale actual reservoir simulation.

Description

Technical Field

[0001] This invention relates to the technical field of reservoir numerical simulation, and in particular to a reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods. Background Technology

[0002] The accuracy of reservoir numerical simulations is often constrained by complex computational domains. Fault cutting alters grid connectivity paths, interlayer displacement causes misalignment between adjacent grids, pinch-out boundaries compress flow channels, and strong gradients around the well and heterogeneous interlayers amplify local conductivity errors. Although these areas account for a small percentage of reserves, they often control injection-production response, breakthrough time, and local pressure distribution.

[0003] To preserve these geological features, geological modeling and reservoir simulation typically use corner grids to describe reservoir geometry. Corner grids express twisted meshes using bar coordinates and corner depths, while carrying information about complex boundaries and fault displacements under regular grid numbering.

[0004] After corner meshes are incorporated into the simulator, the discrete flow equations consist of the mesh, mesh properties, well constraints, connectivity, and connectivity conductivity. Connectivity determines which meshes are allowed to exchange fluid, while conductivity determines the flow capacity between mesh pairs. Conventional mesh methods construct conductivity based on mesh adjacency relationships, offering advantages such as computational stability, clear organization, and suitability for large-scale solutions. However, flow paths in complex local regions depend more on real geometry, mesh misalignment, and abrupt property changes. Simply relying on mesh adjacency connections may underestimate or overestimate local exchange capacity, necessitating the introduction of more flexible computational methods in local applications.

[0005] Meshless methods treat the grid center as a point cloud node and establish local connectivity based on the spatial distribution of nodes. The extended finite volume method combines the point cloud connectivity coefficient, permeability, and node control volume into a conserved conductivity expression. This approach is suitable for high-precision local regions in corner mesh models: the corner mesh provides geometric, volumetric, and physical property data, while the meshless method updates the connectivity and conductivity within the selected region. However, solving for the node control volume through local equations significantly increases the costs of neighborhood search, matrix solving, and data assembly when dealing with meshes of tens of thousands, millions, or even tens of millions of grids. Corner mesh volumes are rigorously recovered from corner geometry. In industrial models where the mesh distortion ratio is finite or local deformation is controllable, directly using corner volumes as node control volumes preserves corner geometric information while controlling computational load.

[0006] While existing global meshless discretization schemes can flexibly handle complex geometries, they suffer from a surge in the number of connections, high computational overhead, and the need to redevelop a complete nonlinear solution system (including core modules such as nonlinear solvers, well control processing, and phase state updates). This results in high implementation and verification costs and significant challenges in engineering deployment. Therefore, there is an urgent need for a technical solution that introduces flexible connectivity and conductivity calculations in local high-precision regions while retaining the large-scale solution capabilities of existing simulators. Summary of the Invention

[0007] The purpose of this invention is to provide a reservoir numerical simulation method based on a hybrid of meshless and meshed methods for conductivity. By retaining the conductivity of the meshed method in the main reservoir region and reconstructing the connection relationship in local high-precision areas using point clouds to calculate the conductivity of the meshless method, this invention solves the problem of insufficient representativeness of conductivity in corner mesh models around wells, near faults, and in strongly heterogeneous zones in existing technologies. At the same time, it overcomes the shortcomings of the global meshless discretization scheme, which leads to a surge in the number of connections and the need to redevelop the nonlinear solution system, resulting in high modification costs.

[0008] To achieve the above objectives, this invention provides a reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed approaches, comprising the following steps: S1. Obtain the reservoir corner mesh model and determine the effective mesh set. The coordinates of the eight corner points of each effective grid are recovered from the corner grid model, and the grid center coordinates and grid volume are calculated. S2. Determine the local meshless conductivity calculation area based on the well perimeter, fault-affected zone, or permeability gradient zone. The effective set of grids outside the region is denoted as ; S3, In the local meshless conductivity calculation region Within, the planar coordinates of the effective grid centers are extracted layer by layer to form the target point cloud. and geometric support point set The generalized finite difference method is used to obtain the discrete coefficients within the influence domain and establish the connection relationships of the point cloud, thus obtaining the local meshless layer connection set. ; S4. Based on the set of connections within the local meshless layer. The extended finite volume method is employed, using the corner volumes of the two end grids as the nodal control volumes. Combined with the generalized finite difference divergence coefficients and the permeability of the two end grids, the conductivity contributions at both ends are calculated. After harmonic merging and unit conversion, the conductivity of each connection is obtained. ; S5. Connect the sets within the local meshless layer. and its conductivity Connection within the regional boundary layer Conventional area grid connection and interlayer connections The organization is a global join table. And expand according to the input format of the reservoir simulator; S6. Invoke the reservoir simulator based on the global connection table. Complete the numerical solutions for pressure, saturation, composition, and well dynamics.

[0009] Preferably, in step S1, the coordinates of the grid center are taken as the arithmetic mean of the coordinates of the eight corner points: ; in For the first The coordinates of the center of each grid cell. For the first The first grid The coordinates of the corner points.

[0010] Preferably, in step S1, the mesh volume is calculated using the trilinear mapping integral of the corner geometry, and any point within the corner mesh is written as: ; in For standard parametric unit coordinates, For a trilinear function, the mapping Jacobian matrix is... for: ; The volume is: ; in For the first Volume of each grid; using Gaussian quadrature calculation.

[0011] Preferably, in step S2, the local meshless conductivity calculation region Select by grid number radius: ; in Number The grid, For an effective set of grids, For the central grid set, The radius of the selected area is used; when multiple neighboring areas are in contact or overlap, the union is taken, and the same grid is only entered into the local area once.

[0012] Preferably, in step S3, the method for constructing the point cloud connection relationship is as follows: for the target node In the domain of influence The radius of the influence domain for searching geometric support points within the layer is determined by both the scale of the point cloud at the same level and the minimum number of support points. ; ; in The point cloud scale coefficient. Geometric support point Distance to the nearest support point on the same floor This refers to the set of grids in the original reservoir grid model. The initial influence domain radius, For the target node The actual radius of influence. To minimize the number of support points, For the target node to the The distance to the nearest support point To ensure that the smallest positive number of the boundary neighbors is included in the influence domain, support points located outside the influence domain do not participate in the node selection. Local fitting; for any geometric support point The formula is as follows: ; ; in, and Support nodes The horizontal and vertical coordinates of the plane, and Target nodes The horizontal and vertical coordinates of the plane, and Support nodes Relative to the target node The difference between the coordinates along the x and y directions.

[0013] Preferably, in step S3, within the influence domain Perform a second-order Taylor expansion on the target node: ; in and They are nodes and nodes Pressure or pressure-related variables on the surface. For u pair First-order partial derivatives 、 For u pair First-order partial derivatives 、 For u pair Second-order partial derivatives 、 For u pair Second-order partial derivatives 、 Let u be the mixed second-order partial derivatives of x and y, respectively; combine the expansions of all neighboring nodes in distance order, and assign node weights. and weight matrix Defined as: ; ; ; in, For the target node With neighboring nodes The planar Euclidean distance between them Let be the planar coordinates of the neighboring nodes. Let be the planar coordinates of the target node. It is the Euclidean norm; For the target node α For neighboring nodes β The weights; Let be the weight matrix of the target node α; This is an operator that forms a diagonal matrix from the elements within the parentheses; For the target node α, its first to second... The weights of each neighboring node; The number of neighboring nodes; nodes The first and second derivative vectors at point Obtained by weighted least squares: ; ; in, For the target node α and its A column vector consisting of the variable differences between neighboring nodes The discrete coefficient matrix is ​​obtained from the weighted least squares problem. For u pair First-order partial derivatives 、 For u pair First-order partial derivatives 、 For u pair Second-order partial derivatives 、 For u pair Second-order partial derivatives 、 Let be the mixed second-order partial derivatives of u with respect to x and y, respectively.

[0014] Preferably, in step S3, the... The combination of discrete coefficients corresponding to the second derivatives in the two directions is: ; in For nodes With neighboring nodes The connection coefficient between them and The coefficient matrix is ​​respectively Rows 3 and 4 correspond to the column elements of node β. For nodes exist The column index in the data is retained to satisfy the following conditions. and Neighborhood column coefficients; nodes The discrete term driven by the pressure difference between each neighboring node is written as: ; in For nodes α Discrete terms driven by pressure difference; the corresponding locally meshless layer internal connectivity set is represented as: ; and This represents any two meshes involved in the connection calculation. , The one-dimensional number of the effective grid.

[0015] Preferably, in step S4, for the connection ( ), calculate the contribution of the control volumes of the two end nodes to the original conductivity of the connection, and calculate the meshless conductivity as follows: ; ; in and They are nodes End and node The original conductivity contribution of the end, superscript This indicates an intermediate quantity that has not yet been converted to a unit. To be 、 The result is obtained by combining the discrete coefficients corresponding to the second derivatives in two directions. To increase penetration rate, and They are nodes and nodes The corresponding corner mesh volume; The harmonic average of the permeability at both ends of the grid is used, and the formula is: ; in and They are nodes and nodes Penetration rate; for Inner and original grid Connections with consistent directional adjacency. and Take the grids at both ends respectively For the original mesh For connections with consistent directional adjacency, the values ​​of the two end meshes are taken respectively. For newly added diagonal connections and non-adjacent connections in the point cloud, first calculate the horizontal equivalent penetration rate of each node: ; ; in and They are nodes α and nodes β Horizontal equivalent penetration rate, For nodes α exist direction , For nodes α exist directional permeability component For nodes β exist directional permeability component For nodes β exist The permeability component in each direction; then, substituting it into the harmonic averaging formula, the bidirectional contribution is used to obtain the meshless conductivity through harmonic averaging: ; in This is the meshless conductivity that is ultimately written into the global join table.

[0016] Preferred, global join table The organization method is as follows: ; The join set is decomposed into: ; in For the global join set, For a locally meshless layer, For connectivity within the region boundary layer, For regular area mesh connections, For interlayer connections; the corresponding conductivity is written as: ; in For connection The final conductivity, The conductivity is calculated using the mesh method in the reservoir simulator; the global connectivity table is expanded according to the reservoir simulator input format, and the non-adjacent connectivity is organized as follows: ; in For the first Non-adjacent connected records, Number the connection record. and The position number is used to connect the two grid ends. The conductivity values ​​are input to the reservoir simulator; for adjacent connections within a local region, they are represented as directional conductivity arrays: ; ; in for Connect adjacent elements of the conductance array in the direction of the connection. for Connect adjacent elements of the conductivity array in the direction of the conductivity; conductivity unit conversion: ; in This is a conversion factor for the unit system.

[0017] Preferably, in step S6, the global join table generated in step S5 is... As the connection relationship input to the reservoir simulator, the adjacent connections in the local area are input through MULTX and MULTY, replacing the corresponding adjacent connection conductivity calculated by default by the simulator. Based on this global connection table, the reservoir simulator combines grid properties, physical property parameters and well control constraints to complete the numerical solution of pressure, saturation, composition and well dynamics in nonlinear iteration.

[0018] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the process of the hybrid meshless and meshed reservoir numerical simulation method of the present invention; Figure 2 This is a map of the local meshless conductivity calculation area selected according to the grid number radius. The solid line grid belongs to the effective grid set, the gray dashed line grid is the invalid grid, and the light gray grid and its center point constitute the local meshless conductivity calculation area. Figure 3This is a schematic diagram of the influence domain connection of a point cloud. The hybrid method only constructs local connections within the influence domain for the target nodes in the local meshless conductivity calculation area. Figure 4 This is a map showing the calculation area of ​​local meshless conductivity around multi-layer wells in the SNARK example of this invention. Gray grids represent inactive grids, black grids represent regular effective grids, blue areas represent local grids around the wells, and black dots are point cloud nodes, where (a) is layer 1, (b) is layer 3, (c) is layer 6, (d) is layer 7, (e) is layer 10, and (f) is layer 11. Figure 5 The first layer of the pressure field comparison diagram at the end of the SNARK example of this invention is shown, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 6 The image shows a comparison of the pressure field at the end of the SNARK calculation in this invention, specifically the third layer, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 7 The image shows a comparison of the pressure field at the end of the SNARK calculation in this invention, specifically the sixth layer, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 8 The image shows a comparison of the pressure field at the end of the SNARK calculation in this invention, specifically the seventh layer, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 9 This is a comparison diagram of the pressure field at the end of the SNARK example of this invention, layer 10, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 10 This is a comparison diagram of the pressure field at the end of the SNARK example of this invention, layer 11, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 11 The image shows a comparison of the oil phase saturation field at the end of the first layer in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 12 The image shows a comparison of the gas phase saturation field at the end of the first layer in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 13 The above is a comparison diagram of the water phase saturation field at the end of the first layer in the SNARK example of this invention, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 14 The image shows a comparison of the oil phase saturation field at the end of the third layer in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 15 The image shows a comparison of the gas phase saturation field at the end of the third layer in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 16 The above is a comparison diagram of the water phase saturation field at the end of the third layer in the SNARK example of this invention, including (a) the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 17 This is a comparison diagram of the oil phase saturation field at the end of layer 10 in the SNARK example of this invention, where (a) is the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 18 The image shows a comparison of the gas phase saturation field at the end of layer 10 in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 19 The image shows a comparison of the water saturation field at the end of layer 10 in the SNARK example of this invention, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 20 This is a comparison diagram of the pressure field of the first layer of the SNARK example of the present invention, which includes (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 21 This is a comparison diagram of the pressure field of the first layer of the SNARK example of this invention, which is 2435 days. Among them, (a) full field grid method, (b) well perimeter local mixing method, and (c) full field horizontal gridless control scheme. Figure 22 This is a comparison diagram of the pressure field of the first layer of the SNARK example of this invention, which is 3652 days. Among them, (a) full field grid method, (b) well perimeter local mixing method, and (c) full field horizontal gridless control scheme. Figure 23 This is a comparison diagram of the oil phase saturation field of the first layer of the SNARK example of this invention after 1217 days, including (a) the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 24 This is a comparison diagram of the oil phase saturation field of the first layer of the SNARK example of this invention, after 2435 days, including (a) the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 25 This is a comparison diagram of the oil phase saturation field in the first layer of the SNARK example of this invention, after 3652 days, including (a) the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 26 This is a comparison diagram of the saturation field of weather phase 1217 in the first layer of the SNARK example of this invention, where (a) is the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 27 This is a comparison diagram of the saturation field of weather phase 2435 in the first layer of the SNARK example of this invention, where (a) is the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 28 This is a comparison diagram of the saturation field of weather phase 3652 in the first layer of the SNARK example of this invention, where (a) is the whole field grid method, (b) the well perimeter local mixing method, and (c) the whole field horizontal gridless control scheme. Figure 29 The image shows a comparison of the water saturation field in the first layer of the SNARK example of this invention, after 1217 days, including (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 30 The image shows a comparison of the water saturation field in the first layer of the SNARK example of this invention, which includes (a) the full-field grid method, (b) the well perimeter local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 31 The image shows a comparison of the water saturation field in the first layer of the SNARK example of this invention, which is 3652 days in total. The comparison includes (a) the full-field grid method, (b) the well-circumferential local mixing method, and (c) the full-domain horizontal gridless control scheme. Figure 32 The SNARK example of this invention shows the pressure and production dynamics curves, where (a) is the average pressure across the entire field, (b) is the daily oil production across the entire field, (c) is the daily water production across the entire field, and (d) is the daily gas production across the entire field. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Specific model specifications need to be selected and determined according to the actual specifications of the device, etc. The specific selection calculation method adopts existing technology in the art, and therefore will not be described in detail.

[0021] Example like Figure 1 As shown, this invention provides a reservoir numerical simulation method based on a hybrid of meshless and meshed methods for conductivity: The flowchart of the hybrid meshless and meshed reservoir numerical simulation method starts from corner grid data, proceeds through local region selection, intra-layer point cloud connection calculation, conductivity calculation, global connection table organization, and finally calls the reservoir simulator to complete the numerical solution.

[0022] I. Model and Calculation Parameter Settings: This embodiment uses the SNARK corner grid model to verify the hybrid conductivity calculation method proposed in this invention. The model has dimensions of 24×25×12 and contains 7200 grids, of which 4328 are effective grids. Multiple production wells are distributed across different stratigraphic levels. The embodiment sets up a full-field grid method, a well-perimeter local hybrid method, and a full-domain horizontal gridless control scheme for comparative verification. Specifically, the connectivity and conductivity of the full-field grid method are calculated by the reservoir simulator based on the original corner grid; the well-perimeter local hybrid method uses the method of this invention to reconstruct intra-layer connectivity in a local area around the well and calculate the gridless conductivity; the full-domain horizontal gridless control scheme constructs horizontal gridless connections within all effective grids, generating a total of 15938 intra-layer connection records, including 4030 adjacent connections in the x-direction, 4118 adjacent connections in the y-direction, and 7790 non-adjacent connections. The full-field grid method and the full-domain horizontal gridless control scheme are used to compare and verify the results of the embodiment.

[0023] In this embodiment, the values ​​for each parameter are: point cloud scale coefficient. =1.75, minimum number of support points =6, influence domain boundary envelope =0.001. Wherein =1.75 makes the initial influence domain cover the axial and oblique nearest neighbors at the current point cloud scale; the two-dimensional second-order expansion contains five derivative components to be solved, and six support points are taken to form a weighted least squares problem with a redundancy condition; Only for the purpose of the first The nearest neighbor is stably included in the influence domain without changing the main scale of the point cloud.

[0024] The perforation grid is read using the COMPDAT keyword in the Eclipse simulator and assembled into a central grid set. Take the radius of the grid number. =3, select a local region within the layer containing each perforation grid. Union the well-perimeter neighborhoods that are in contact, overlapping, or contained with each other to form a local meshless conductivity calculation region around the well. .

[0025] S1. Obtain the reservoir corner grid model, determine the effective grid set, recover the coordinates of the eight corner points of each effective grid from the corner grid model, and calculate the grid center coordinates and grid volume.

[0026] Assume the reservoir grid model consists of a set of grids. Composition, in which , , These represent the grid numbers in the three directions. The effective grid set is defined as follows: To avoid mixing up the 3D mesh number with the numbers at both ends of the connection, use... and This represents any two meshes involved in the connection calculation. 、 The one-dimensional number of the effective grid.

[0027] The corner mesh reconstructs the eight corners of each grid using cylindrical coordinates and corner depths. For any valid grid... In this embodiment, the arithmetic mean of the coordinates of the eight corner points is defined as the point cloud node coordinates: ; in For the first The coordinates of the center of each grid cell. For the first The first grid The coordinates of the corner points. The center coordinates are used for subsequent point cloud connection searches.

[0028] The mesh volume is calculated using trilinear mapping integrals of corner geometry. Let the coordinates of the standard parametric element be... Any point within the corner grid is written as: ; in It is a trilinear function. If the first... The local coordinates of the corner points on the standard parameter unit are labeled as follows: ,and Then the trilinear shape function is the product of one-dimensional linear shape functions in the three coordinate directions: ; The mapping Jacobian matrix is: ; use The Gaussian quadrature calculation yields the volume as follows: ; in: ; This quadrature form is applicable to trilinear corner mapping and is used for calculating the volume of corner meshes with distortion.

[0029] In this embodiment, the difference in the total conductivity obtained by the two volume calculation methods (Ponting / Gauss corner volume method and hexahedral partitioning method) in the SNARK example is small. Ponting / Gauss volume is used as the default setting, and the calculation time is about 0.027s.

[0030] S2, see also Figure 2 For the local meshless conductivity calculation region selected by grid number radius, the solid line grid belongs to the effective grid set. The gray dashed grid represents an invalid grid, and the light gray grid and its center point constitute a local meshless conductivity calculation region. .

[0031] The local meshless conductivity calculation area is determined based on the well perimeter, fault-affected zone, or permeability gradient zone. The effective set of grids outside the region is denoted as .

[0032] Local Meshless Conductivity Calculation Region The determination is based on the well perimeter, fault influence range, permeability gradient, sensitive areas for historical data fitting, or numerical accuracy requirements. In this embodiment, the local area is determined based on the well perimeter. Let the central grid set requiring focused processing be... The radius of the selected area is If there are 1 grid, then the local region of the same layer is defined as: ; in Number The grid, For an effective set of grids, For the central grid set, The radius of the selected area is denoted by . This formula selects a two-dimensional neighborhood within the layer containing the central grid based on the grid number distance, facilitating direct connection between local areas and the original grid.

[0033] This embodiment processed a total of 18 perforated well sections, and took... Select a local region within the layer containing each perforation grid according to the above formula. When different well perimeter regions are in contact, overlap, or contain each other, The union of these elements is taken, ensuring that the same mesh enters the local region only once, and the same connection is recorded only once in the local connectivity table. Ultimately, 838 locally effective meshes are selected as... Because the well sections are distributed across multiple stratigraphic levels, Different local morphologies are presented at different layers; see [link / reference] Figure 4 Gray grids represent inactive grids, black grids represent regular active grids, blue areas represent local grids around the well, and black dots represent point cloud nodes. Figure 4 (a) is the first layer, Figure 4(b) is the 3rd floor. Figure 4 (c) is the 6th floor. Figure 4 (d) refers to the 7th floor. Figure 4 (e) is the 10th floor. Figure 4 (f) is the 11th floor.

[0034] The set of effective grids outside the region is denoted as After constructing the point cloud layer by layer, both ends are located in The intralayer connections are included in the meshless conductivity calculation; one end is located at The other end is located The conductivity of intralayer boundary connections follows the simulator mesh method; the conductivity of interlayer connections follows the simulator calculation based on the original corner mesh.

[0035] S3, In the local meshless conductivity calculation region Within, the planar coordinates of the effective grid centers are extracted layer by layer to form the target point cloud. and geometric support point set The generalized finite difference method is used to obtain the discrete coefficients within the influence domain and establish the connection relationships of the point cloud, thus obtaining the local meshless layer connection set. .

[0036] Within the selected area, the planar coordinates of the local grid centers are used to construct the target point cloud, while all grid centers in the same layer are retained as geometric support points. ; in , Only participating in local polynomial fitting, the two ends of the final flow connection are both taken from... This embodiment only constructs horizontal point cloud connections within the same layer.

[0037] For the target node In the domain of influence The geometric support points are searched internally, and local discrete coefficients are constructed using the generalized finite difference method. The geometric neighborhood of the target node is written as: ; The radius of the influence domain is determined by the scale of the point cloud at the same layer and the minimum number of support points. Let... Geometric support point The distance to the nearest support point in the same layer is used to characterize the typical spacing of the point cloud in that layer, with the median of the nearest neighbor distances, and an initial radius is defined: ; in For point cloud scale coefficients, this embodiment takes... =1.75. This value ensures that the initial influence domain covers both axial and diagonal neighbors at the current point cloud scale.

[0038] set up This represents the minimum number of support points required for a local second-order fit. For the target node to the The actual influence zone radius is determined by the distance from the nearest support point: ; in For the target node The actual radius of influence. To minimize the number of support points, this embodiment takes... =6. The two-dimensional second-order expansion contains five derivative components, and six support points are used to form a weighted least squares problem with one redundant condition; To ensure that the boundary neighbors are included in the influence domain by a small positive number, this embodiment takes... =0.001, only used to set the first The nearest neighbors are stably included in the influence domain without changing the main scale of the point cloud. External support points do not participate in the node. Local fitting.

[0039] make Represents a node The pressure or pressure-related variables on the surface. For any geometric support point. ,remember: ; ; At the node to Performing a second-order Taylor expansion, we get: ; in, For u pair First-order partial derivatives 、 For u pair First-order partial derivatives 、 For u pair Second-order partial derivatives 、 For u pair Second-order partial derivatives 、 Let be the mixed second-order partial derivatives of u with respect to x and y, respectively; Combine the expansions of all neighboring nodes in distance order. Node weights and the weight matrix are defined as follows: ; ; ; in Let be the number of neighboring nodes. The first and second derivative vectors at node α are obtained by weighted least squares: ; ; in It is obtained from the corresponding weighted least squares problem.

[0040] In the spatial discretization of the reservoir flow equation, the pressure difference contributes to the internode flow through the second-order spatial term, therefore... , The combination of discrete coefficients corresponding to the second derivatives in the two directions is: ; in For nodes With neighboring nodes The connection coefficient between them and The coefficient matrix is ​​respectively The third and fourth rows correspond to the node. β Column elements, For nodes exist The column number in the [database name]. When constructing flow connections, only those satisfying [the specified condition] are retained. and The neighborhood column coefficients; node The discrete term driven by the pressure difference between each neighboring node is written as: ; The corresponding set of interconnected local meshless layers is represented as follows: ; This formula indicates that connections are determined by the neighborhood relationships of the point cloud, and nodes outside the influence domain are not connected. Layer-by-layer point cloud search ensures that the connectivity relationships are consistent with the intra-layer connectivity organization of layered geological models and commonly used reservoir simulators. See also Figure 3 The diagram illustrates the connection of the influence domain of a point cloud. The hybrid method constructs local connections within the influence domain for the target nodes only in the local meshless conductivity calculation area.

[0041] exist In the process, the meshless method generated 2742 intra-layer connection records; of which 736 were... Adjacent connections in the direction, 724 in total There are 1282 non-adjacent connections.

[0042] S4. Based on the set of connections within the local meshless layer. The extended finite volume method is employed, using the corner volumes of the two end grids as the nodal control volumes. Combined with the generalized finite difference divergence coefficients and the permeability of the two end grids, the conductivity contributions at both ends are calculated. After harmonic merging and unit conversion, the conductivity of each connection is obtained. .

[0043] The flow rate between reservoir grids is determined by both the pressure difference at both ends and the connectivity conductivity. The extended finite volume method will use coefficients that reflect the spatial relationships between nodes. This, combined with permeability and node control volume, yields the connectivity conductivity between point cloud nodes. For connectivity... First, calculate the contribution of the control volumes at both ends of the node to the original conductivity of the connection.

[0044] node Controlled volume Directly taking the volume of the corresponding corner mesh, and following the extended finite volume method, solving for the node control volume through a system of local equations would significantly increase the costs of neighborhood search, local matrix solving, and data assembly. Therefore, this invention uses the corner mesh volume directly as the node control volume. In this embodiment, the volume used is the Ponting / Gauss corner volume.

[0045] For connection ( α,β The original conductivity contributions at both ends are: ; ; in and They are nodes End and node The original conductivity contribution of the end, superscript This indicates an intermediate quantity that has not yet been converted to a unit. To increase penetration rate, and They are nodes and nodes The corresponding corner mesh volume; Harmonic average of the permeability at both ends of the grid is used: ; in and They are nodes and nodes Penetration rate.

[0046] for Inner and original grid Connections with consistent directional adjacency. and Take the grids at both ends respectively For the original mesh For connections with consistent directional adjacency, the values ​​of the two end meshes are taken respectively. .

[0047] For newly added diagonal connections and other non-adjacent connections in the point cloud, first calculate the horizontal equivalent penetration rate of each node: ; ; in and These represent the horizontal equivalent permeability of nodes α and β, respectively. The corresponding nodes are respectively direction and The permeability component in each direction is then substituted into the harmonic averaging formula. The bidirectional contribution is then used to obtain the meshless conductivity using harmonic averaging: ; in This is the meshless conductivity that is ultimately written into the global join table.

[0048] This embodiment compares three oblique connection horizontal equivalent permeability modes: geometric mean, arithmetic mean, and directional projection. 、 Adjacent connections are used with the harmonic mean of the permeability of the corresponding direction. The number of NNC, TRANX, and TRANY connections remains consistent across the three modes, at 7790, 4030, and 4118 connections, respectively. The total conductivity of oblique connections (NNC) is 12512.21, 12513.12, and 12512.96 (FIELD Eclipse input units), with minimal differences. In this embodiment, the geometric mean is used as the default setting.

[0049] S5. Connect the sets within the local meshless layer. and its conductivity Connection within the regional boundary layer Conventional area grid connection and interlayer connections The organization is a global join table. And expand according to the input format of the reservoir simulator.

[0050] The core output of the hybrid method is a global join table: ; in To join sets, This corresponds to the conductivity of the connections. To clarify the origin of different connections, the connection set is decomposed into: ; in for Internal connections are derived from the point cloud connectivity. To cross and Intralayer boundary connection, for Adjacent internal grids are connected. This refers to interlayer connections. In this embodiment, there are a total of 1155 connections within the boundary layer, which are calculated by the reservoir simulator based on the original corner grid.

[0051] The corresponding conductivity is written as: ; in For connection The final conductivity, The conductivity is calculated using the grid method in the reservoir simulator.

[0052] The global join table is expanded according to the reservoir simulator input format. Taking the Eclipse class input as an example, the non-adjacent joins are organized as follows: ; in For the first Non-adjacent connected records, Number the connection record. and The position number is used to connect the two grid ends. The conductivity values ​​are input to the reservoir simulator; in this embodiment, 1282 non-adjacent connections are written into the NNC class connection table. Adjacent connections within a local region are written as a directional conductivity array: ; ; in for Connect adjacent elements of the conductance array in the direction of the connection. for Adjacent elements in the conductivity array are connected in a specific direction. In this embodiment... 736 adjacent connections were written to the TRANX array. 724 directional adjacency connections are written to the TRANY array. Local adjacency connections are established using MULTX and MULTY. Internal corresponding position input, region boundary and The internal connectivity maintains the conductivity calculated by the simulator using the mesh method. Conductivity units are converted to: ; in This is a unit conversion factor. This embodiment uses the FIELD unit system. METRIC unit preparation .

[0053] S6. Call the reservoir simulator and complete the numerical solutions for pressure, saturation, composition and well dynamics based on the global connection table.

[0054] The global connection table generated in step S5 The connectivity data is input to the reservoir simulator. Adjacent connections within a local region are input via MULTX and MULTY, replacing the simulator's default calculation of the corresponding adjacent connection conductivity. Based on this global connectivity table, and combined with mesh attributes, physical parameters, and well control constraints, the reservoir simulator completes the numerical solutions for pressure, saturation, composition, and well dynamics in nonlinear iterations. The core modules of mature reservoir simulators, such as the nonlinear solver, well control processing, and phase update, remain unchanged; only the connectivity data is replaced with the global connectivity table calculated in this invention.

[0055] This embodiment simulates 3652 days, with the final report step being step 126, corresponding to January 1, 2008. The Eclipse CPU time was 7.04 seconds, and the total runtime was 7.04 seconds. Corner grid lines were preserved during field plotting; the saturation variable was only clipped to the physical range [0,1] during the plotting phase to eliminate 10 in the solution output. -4 If a small value exceeds the limit, the original simulation result file will not be modified.

[0056] like Figures 5-10 As shown, in the pressure field comparison of multiple layers at the end of the well, layers 1, 3, 6, 7, 10, and 11 cover the lower, middle, and upper layers with more obvious distribution around the well. Multilayer results are used to examine whether the selection of local areas and the response of the field map are concentrated near the corresponding perforated layer, while avoiding conclusions based on only a single layer.

[0057] like Figures 11-19 As shown, the three sets of three-phase saturation field diagrams of the 1st, 3rd, and 10th layers at the end of the period cover the multiphase distribution of the lower, middle, and upper layers. This is used to supplement the pressure field evidence and to check whether the influence of the local connection table on the three-phase front of oil, gas, and water remains consistent with the layer changes.

[0058] like Figures 20-22As shown, in the pressure field comparison at three time points in the first layer, all three schemes used the same corner geometry and well control system. The differences stemmed from local connectivity and conductivity calculation methods. The well-perimeter local hybrid method introduces local connectivity reconstruction based on the full-field grid method response, causing local adjustments to the pressure propagation pattern near the well perimeter region. The full-domain horizontal gridless comparison scheme shows the overall field map changes when the intra-layer connectivity of all effective grids is given by point cloud relationships.

[0059] like Figures 23-25 As shown, the oil phase saturation field of the first layer is compared. The saturation of the oil, gas, and water phases together reflects the influence of the local connectivity table on the multiphase front propulsion and wellbore displacement capacity. The changes in the oil phase field are complementary to the pressure and water phase fields, and can be used to determine whether the local mixing conductivity has changed the wellbore displacement response.

[0060] like Figures 26-28 As shown, the gas phase saturation field of the first layer is compared. The local mixing scheme preserves the transmissibility of the grid method outside the region and around the well. An internal meshless connectivity table is introduced, so the main differences in gas phase distribution are concentrated around the well and in the adjacent flow channels. This result reflects the basic positioning of the hybrid method: most reservoir regions maintain the grid-based discrete structure of the original simulator, while local high-precision areas are responsible for connectivity and conductivity reconstruction.

[0061] like Figures 29-31 As shown, the first layer of water phase saturation field comparison is presented. The water phase field and the gas phase field together reflect the impact of local connectivity remodeling on frontal propulsion and wellbore exchange capacity. Due to The external connectivity still uses the grid method for conductivity, while the field map changes of the local hybrid method are mainly concentrated near the well perimeter, and the overall computational domain shape remains consistent with the original corner grid boundary.

[0062] like Figure 32 As shown, the production dynamics and overall pressure curves demonstrate that the local hybrid method can incorporate local meshless permeability while maintaining well control in the reservoir simulator and the nonlinear solution process. Further sensitivity analyses can be conducted focusing on well radius, local morphology, and oblique permeability patterns to quantify the impact of local connectivity reconstruction on well dynamics.

[0063] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A numerical simulation method for oil reservoirs based on a hybrid conductivity method combining meshless and meshed methods, characterized in that: Includes the following steps: S1. Obtain the reservoir corner mesh model and determine the effective mesh set. The coordinates of the eight corner points of each effective grid are recovered from the corner grid model, and the grid center coordinates and grid volume are calculated. S2. Determine the local meshless conductivity calculation area based on the well perimeter, fault-affected zone, or permeability gradient zone. The effective set of grids outside the region is denoted as ; S3, In the local meshless conductivity calculation region Within, the planar coordinates of the effective grid centers are extracted layer by layer to form the target point cloud. and geometric support point set The generalized finite difference method is used to obtain the discrete coefficients within the influence domain and establish the connection relationships of the point cloud, thus obtaining the local meshless layer connection set. ; S4. Based on the set of connections within the local meshless layer. The extended finite volume method is employed, using the corner volumes of the two end grids as the nodal control volumes. Combined with the generalized finite difference divergence coefficients and the permeability of the two end grids, the conductivity contributions at both ends are calculated. After harmonic merging and unit conversion, the conductivity of each connection is obtained. ; S5. Connect the sets within the local meshless layer. and its conductivity Connection within the regional boundary layer Conventional area grid connection and interlayer connections The organization is a global join table. And expand according to the input format of the reservoir simulator; S6. Invoke the reservoir simulator based on the global connection table. Complete the numerical solutions for pressure, saturation, composition, and well dynamics.

2. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 1, characterized in that: In step S1, the coordinates of the grid center are calculated using the arithmetic mean of the coordinates of the eight corner points: ； in For the first The coordinates of the center of each grid cell. For the first The first grid The coordinates of the corner points.

3. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 2, characterized in that: In step S1, the mesh volume is calculated using the trilinear mapping integral of the corner geometry, and any point within the corner mesh is written as: ； in For standard parametric unit coordinates, For a trilinear function, the mapping Jacobian matrix is... for: ； The volume is: ； in For the first Volume of each grid; using Gaussian quadrature calculation.

4. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 1, characterized in that: In step S2, the local meshless conductivity calculation region Select by grid number radius: ； in Number The grid, For an effective set of grids, For the central grid set, The radius of the selected area is used; when multiple neighboring areas are in contact or overlap, the union is taken, and the same grid is only entered into the local area once.

5. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 1, characterized in that: In step S3, the method for constructing the point cloud connectivity is as follows: for the target node In the domain of influence The radius of the influence domain for searching geometric support points within the layer is determined by both the scale of the point cloud at the same level and the minimum number of support points. ； ； in The point cloud scale coefficient. Geometric support point Distance to the nearest support point on the same floor This refers to the set of grids in the original reservoir grid model. The initial influence domain radius, For the target node The actual radius of influence. To minimize the number of support points, For the target node to the The distance to the nearest support point To ensure that the smallest positive number of the boundary neighbors is included in the influence domain, support points located outside the influence domain do not participate in the node selection. Local fitting; for any geometric support point The formula is as follows: ； ； in, and Support nodes The horizontal and vertical coordinates of the plane, and Target nodes The horizontal and vertical coordinates of the plane, and Support nodes Relative to the target node The difference between the coordinates along the x and y directions.

6. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 5, characterized in that: In step S3, within the influence domain Perform a second-order Taylor expansion on the target node: ； in and They are nodes and nodes Pressure or pressure-related variables on the surface. For u pair First-order partial derivatives 、 For u pair First-order partial derivatives 、 For u pair Second-order partial derivatives 、 For u pair Second-order partial derivatives 、 Let u be the mixed second-order partial derivatives of x and y respectively; combine the expansions of all neighboring nodes in distance order, and let the node weights be... and weight matrix Defined as: ； ； ； in, For the target node With neighboring nodes The planar Euclidean distance between them Let be the planar coordinates of the neighboring nodes. Let be the planar coordinates of the target node. It is the Euclidean norm; For the target node α For neighboring nodes β The weights; Let be the weight matrix of the target node α; This is an operator that forms a diagonal matrix from the elements within the parentheses; For the target node α, its first to second... The weights of each neighboring node; The number of neighboring nodes; nodes The first and second derivative vectors at point Obtained by weighted least squares: ； ； in, For the target node α and its A column vector consisting of the variable differences between neighboring nodes The discrete coefficient matrix is ​​obtained from the weighted least squares problem. For u pair First-order partial derivatives 、 For u pair First-order partial derivatives 、 For u pair Second-order partial derivatives 、 For u pair Second-order partial derivatives 、 Let be the mixed second-order partial derivatives of u with respect to x and y, respectively.

7. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 6, characterized in that: In step S3, The combination of discrete coefficients corresponding to the second derivatives in the two directions is: ； in For nodes With neighboring nodes The connection coefficient between them and The coefficient matrix is ​​respectively Rows 3 and 4 correspond to the column elements of node β. For nodes exist The column index in the data is retained to satisfy the following conditions. and Neighborhood column coefficients; nodes The discrete term driven by the pressure difference between each neighboring node is written as: ； in For nodes α Discrete terms driven by pressure difference; the corresponding locally meshless layer internal connectivity set is represented as: ； in and This represents the index of any two meshes involved in the connection calculation. , The one-dimensional number of the effective grid.

8. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 7, characterized in that: In step S4, for the connection ( ), calculate the contribution of the control volumes of the two end nodes to the original conductivity of the connection, and calculate the meshless conductivity as follows: ； ； in and They are nodes End and node The original conductivity contribution of the end, superscript This indicates an intermediate quantity that has not yet been converted to a unit. To be 、 The result is obtained by combining the discrete coefficients corresponding to the second derivatives in two directions. To increase penetration rate, and They are nodes and nodes The corresponding corner mesh volume; The harmonic average of the permeability at both ends of the grid is used, and the formula is: ； in and They are nodes and nodes Penetration rate; for Inner and original grid Connections with consistent directional adjacency. and Take the grids at both ends respectively For the original mesh For connections with consistent directional adjacency, the values ​​of the two end meshes are taken respectively. For newly added diagonal connections and non-adjacent connections in the point cloud, first calculate the horizontal equivalent penetration rate of each node: ； ； in and They are nodes α and nodes β Horizontal equivalent penetration rate, For nodes α exist directional permeability component For nodes α exist directional permeability component For nodes β exist directional permeability component For nodes β exist The permeability component in each direction; then, substituting it into the harmonic averaging formula, the bidirectional contribution is used to obtain the meshless conductivity through harmonic averaging: ； in This is the meshless conductivity that is ultimately written into the global join table.

9. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 1, characterized in that: Global join table The organization method is as follows: ； in For connection The final conductivity, and This represents any two meshes involved in the connection calculation. , One-dimensional numbering of the effective grid. The global join set; the join set is decomposed into: ； in For a locally meshless layer, For connectivity within the region boundary layer, For regular area mesh connections, For interlayer connections; the corresponding conductivity is written as: ； in The conductivity is calculated using the mesh method in the reservoir simulator; the global connectivity table is expanded according to the reservoir simulator input format, and the non-adjacent connectivity is organized as follows: ； in For the first Non-adjacent connected records, Number the connection record. and The position number is used to connect the two grid ends. The conductivity values ​​are input to the reservoir simulator; for adjacent connections within a local region, they are represented as directional conductivity arrays: ； ； in for Connect adjacent elements of the conductance array in the direction of the connection. for Connect adjacent elements of the conductivity array in the direction of orientation; conductivity units are converted to: ； in This is a conversion factor for the unit system.

10. The reservoir numerical simulation method based on a hybrid conductivity method combining meshless and meshed methods according to claim 1, characterized in that: In step S6, the global join table generated in step S5 is... As the connection relationship input to the reservoir simulator, the adjacent connections in the local area are input through MULTX and MULTY, replacing the corresponding adjacent connection conductivity calculated by default by the simulator. Based on this global connection table, the reservoir simulator combines grid properties, physical property parameters and well control constraints to complete the numerical solution of pressure, saturation, composition and well dynamics in nonlinear iteration.

Images