Oil reservoir meshless streamline simulation method

Abstract

The invention discloses an oil reservoir meshless streamline simulation method, which belongs to the technical field of oil reservoir numerical simulation and comprises the following steps: constructing meshless point cloud of an oil reservoir computational domain and calculating node control volume; solving node pressure distribution based on an extended finite volume method, and calculating a seepage velocity field by adopting a partition strategy; obtaining a streamline and flight time through piecewise linear tracking, and reducing the dimension of the saturation equation to the streamline for solving; the node average water saturation is obtained through saturation mapping; by the adoption of the reservoir meshless streamline simulation method, limitation of a traditional grid method and numerical dissipation of a meshless method are effectively overcome, saturation calculation precision and simulation efficiency are remarkably improved through innovative speed calculation and streamline tracking, and excellent adaptability to complex reservoir conditions is shown.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of oil reservoir numerical simulation, and in particular to an oil reservoir meshless streamline simulation method. BACKGROUND

[0002] Oil reservoir numerical simulation is a core technical means for studying the seepage law of underground multiphase fluid, and the key lies in high-precision and high-efficiency discrete solution of the seepage control equation set. After long-term development, the finite volume discrete method based on the two-point flux approximation (TPFA) format which meets the local mass conservation, has compactness and monotonicity, has been widely used in various commercial and academic oil reservoir numerical simulators. However, the TPFA format has significant numerical errors when dealing with full tensor permeability fields or complex non-structured grids. Therefore, more advanced discrete methods such as multi-point flux approximation (MPFA) format, mixed finite element method (MFE) and simulated finite difference method (MFD) have been developed to improve the calculation accuracy under complex geological conditions. On the other hand, with the increase of the complexity of the geometry of the oil reservoir, the traditional Cartesian grid is difficult to apply, and the complex grid such as the corner grid and the non-structured grid can improve the discrete precision, but when dealing with complex oil reservoir geometric boundaries, it often faces problems such as difficulty in grid generation and difficulty in ensuring grid quality. The meshless method is discretized by point cloud, which effectively avoids the complex grid topology construction process and shows significant advantages in adaptability in complex regions. Therefore, it has important engineering application value to develop a meshless oil reservoir simulator based on point cloud discretization.

[0003] Among them, the generalized finite difference method (GFDM) as a typical meshless method directly approximates the spatial derivative through node arrangement, Taylor expansion and weighted least squares method, and has been successfully applied to the solution of many scientific and engineering problems. In recent years, GFDM has been further combined with the concept of node control volume to develop the extended finite volume method (EFVM), which can directly extract the conductivity between nodes while ensuring local mass conservation, and thus compatible with the nonlinear solution process of existing oil reservoir simulators, laying a foundation for the construction of a practical meshless oil reservoir simulator.

[0004] Although the oil reservoir numerical simulation method is relatively mature in the solution of pressure equation, when dealing with saturation or component concentration equation dominated by convection, both traditional grid methods and emerging meshless methods face significant numerical dissipation problems. Although this problem can be alleviated by greatly increasing the grid or node, this will lead to a sharp increase in calculation cost, which is difficult to meet the needs of large-scale practical applications in the field.

[0005] In traditional grid simulators, two main approaches are typically employed to mitigate numerical dissipation caused by convection terms: First, higher-order schemes such as the discontinuous Galerkin method (DG) are introduced to improve the accuracy of saturation calculations. Examples include the developed MFE-DG and MFD-DG hybrid methods. However, these methods have stringent stability requirements, often necessitating smaller time steps and significantly increasing computational burden. Second, streamline simulation methods are developed. These methods effectively improve the resolution and computational efficiency by solving the saturation equation along streamlines. Since Pollock proposed the three-dimensional streamline tracing algorithm, the streamline method has been gradually extended to consider complex physical processes such as gravity, capillary forces, component simulation, and thermal recovery. It has been widely applied and continuously developed in structured grids, unstructured grids, and various fractured reservoir models.

[0006] However, it should be noted that although the streamline method has formed a relatively complete system in grid-based reservoir simulation, no streamline simulation technology has yet emerged to be combined with it in the emerging gridless reservoir numerical simulation methods. Given the potential advantages of gridless methods in terms of adaptability to complex regions and simulation accuracy, this paper introduces the streamline method into the gridless simulation framework, constructing the first gridless reservoir streamline simulation method, abbreviated as EFVM-SL. This has significant theoretical and engineering value for significantly reducing numerical dissipation, improving the theoretical system of gridless reservoir simulation, and promoting its practical industrial application. Summary of the Invention

[0007] The purpose of this invention is to provide a meshless streamline simulation method for oil reservoirs to solve the problems of calculating seepage velocity distribution, streamline tracing, time-of-flight calculation, and high-resolution saturation calculation under a meshless framework.

[0008] To achieve the above objectives, the present invention provides a meshless streamline simulation method for oil reservoirs, comprising the following steps: S1. Generate a meshless point cloud in the reservoir computational domain, and calculate the control volume of each node in the point cloud based on the extended finite volume method. S2. Based on the control volume, the extended finite volume method is used to discretize the two-phase flow equation of porous media and solve for the nodal pressure distribution in the meshless point cloud. S3. Calculate the node seepage velocity based on the node pressure distribution. For nodes containing injection and production wells, use the radial seepage velocity model. For ordinary nodes, use the generalized finite difference method to calculate the pressure gradient to obtain the seepage velocity. S4. Perform piecewise linear streamline tracing based on seepage velocity in a gridless point cloud to determine the streamline trajectory and flight time along the streamline. S5. Based on the definition of flight time, the two-dimensional / three-dimensional convection transport problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate. The water saturation equation is solved along each streamline to obtain the saturation distribution on the streamline. S6. Based on the relationship between streamline segments and node control areas, map streamline saturation to node control volumes and calculate the average water saturation within the node control areas.

[0009] Preferably, the specific steps of S1 include: S11. Discretize the reservoir computational domain using a meshless method to obtain the node set. This forms an initial gridless point cloud, in which, The total number of nodes; S12 represents each node in the point cloud. The local point cloud is determined by using a circular influence domain. and neighboring node set ,in, Based on the generalized finite difference method, each node is defined using the extended finite volume method. A dedicated control area; ; in, For nodes The control domain, For nodes Controlling volume, The volume of the computational domain; S13, For nodes Determine whether it is a node Connectable point clouds, if they satisfy This constructs a connectable point cloud. Based on the connectable point cloud, a pair is created for each node. Construct a linear equation concerning the control volume; ; in, For nodes Controlling volume, , , , For nodes and nodes The relevant discrete coefficients of the generalized finite difference method; S14, if node Add virtual nodes outside the computation domain to serve as boundary nodes. The control volume of boundary nodes is corrected based on virtual nodes; The coordinates of the virtual node are: ; in, and They are nodes and nodes coordinates Boundary nodes The unit outward normal vector at that location, For nodes Local point cloud Node-to-node Weighted average distance, i.e. , For nodes Relative to node The weight function value; The revised control volume formula is: ; S15. Based on the nodal characteristic angles of the boundary nodes. The control volume equations of the boundary nodes are modified, and the control volume equations of all node pairs and the boundary nodes are combined to form an overdetermined linear system of equations. ; in, To calculate the total number of node pairs that can be formed within the domain, These are weighting coefficients. , ; S16. Solve the overdetermined linear equations using the least squares method to obtain the control volume value for each node.

[0010] Preferably, the specific steps of S2 include: S21. Obtain the pressure field using the generalized finite difference method. exist First spatial derivative at , and second-order spatial derivative , , ; ; ; ; ; ; in, , , , , For discrete coefficients, For nodes Pressure value, As the central node Pressure value; S22. Integrating the two-phase flow equations for porous media over the nodal control volumes and discretizing them in time using an implicit scheme, the specific formulas obtained are as follows: ; in, For divergence operators, For oil phase fluidity, For water phase mobility, For absolute penetration rate, For pressure gradient, For nodes The spatial region corresponding to the control volume is the node. With nodes The harmonic average of absolute permeability. For time steps At that time, node With nodes The relative permeability of the oil phase between them For time steps At that time, node With nodes The relative permeability of the water phase between them For time steps At that time, node With nodes The viscosity of the oil phase between them For time steps At that time, node With nodes The viscosity of the aqueous phase between them, For time steps At that time, node Pressure value, For time steps At that time, node Pressure value; S23. The phase permeability and viscosity are calculated using the single-point upstream weighted formula, and the absolute permeability is calculated using the harmonic average formula. The calculation formulas are as follows: ; ; ; S24. Extracting internode conductivity using an extended finite volume discretization scheme based on the two-phase flow equation in porous media. And process source and sink items ; ; ; S25. Based on the extended finite volume method discretization scheme, the closed boundary is processed. In the pressure discretization equation of the boundary node, the conductivity term between it and the virtual node is ignored, and the node control volume is corrected to the real control volume. Discrete equations are established for all real nodes to form a closed linear equation system. The pressure distribution of all nodes is obtained by solving the system.

[0011] Preferably, the extended finite volume method discretization scheme for the two-phase flow equation of porous media in S24 is as follows: ; Introducing conductivity simplifies to: ; in, For time steps At that time, node Source and sink items in the oil phase, For time steps At that time, node Sources and sinks of water phases, This represents the overall compression ratio.

[0012] Preferably, the specific steps of S3 include: S31. Divide the nodes into ordinary nodes that do not contain injection and production wells. Sets of source and sink nodes that include source and sink items such as injection and production wells. ; S32, For ordinary nodes The first spatial derivative of the pressure field and the nodal parameters obtained based on the generalized finite difference method The pressure value reaches the node The first spatial derivative of the pressure function is used to estimate the seepage velocity at the node using Darcy's law, combined with the total mobility distribution. The calculation formula is as follows: ; ; in, Local point cloud scale of nodes , for directional seepage velocity components, for directional seepage velocity components, For nodes absolute penetration rate; S33, For source and sink nodes By controlling its volume to be equivalent to a cylinder, its equivalent radius is calculated. The radial seepage velocity value at this node is obtained. This allows for the calculation of the seepage velocity at the nodes. The calculation formula is as follows: ; ; ; in, To calculate the coordinates of any node to be interpolated within the domain, For nodes The coordinates; S34. Integrate the seepage velocities of all ordinary nodes and source / sink nodes, and obtain the seepage velocity at any point in the computational domain through linear interpolation. This forms a continuous seepage velocity field within the reservoir computational domain, calculated using the following formula: ; in, For nodes At the point to be interpolated Linear interpolation basis functions, For nodes At the point to be interpolated Linear interpolation basis functions.

[0013] Preferably, the specific steps of S4 include: S41. Taking the injection well node as the center, and... On a circle with radius , the starting points of the streamlines are distributed at equal intervals according to the required number of streamlines; S42. Determine the starting point of the streamline as (Note: is missing from the original text). Calculate its seepage velocity Find the node closest to the starting point of the streamline. Obtain the porosity of the node. Set a small distance Calculate flight time using the single tracking distance step size. The streamline position after the movement is obtained. ; ; ; S43, with Repeat step S42 to obtain new nodes. and corresponding flight time The process continues until the streamline reaches the node where the production well is located, completing the single streamline tracing. ; .

[0014] Preferably, the specific steps of S5 include: S51. Transform the two-dimensional / three-dimensional convection transport problem into a one-dimensional flow branch function problem through streamline simulation; ; in, Porosity Water saturation Moisture content, For total phase velocity, The distance along the streamline; S52. Based on the definition of flight time, the above one-dimensional split function problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate. Flight time is defined as: ; The rewritten one-dimensional constant coefficient equation is: ; S53. Solve the above one-dimensional equation using numerical methods to obtain the water saturation at each point on the streamline. ; in, This refers to the local time step when calculating water saturation along the streamline. This represents the flight time step.

[0015] Preferably, the formula for calculating the average water saturation within the node control area in S6 is as follows: ; in, For the nodes that flow through The streamline set of the control area For set The first in A streamline.

[0016] Therefore, the above-mentioned method for simulating reservoir streamlines has the following advantages: (1) This method is the first to effectively combine the meshless method with streamline simulation, overcoming the limitation of the lack of streamline technology in traditional meshless reservoir simulation. This method can be stably calculated in different types of point clouds and has the advantage of flexible computational domain discretization.

[0017] (2) An innovative method for calculating the radial velocity of an equivalent cylinder with injection and production well nodes was proposed, which effectively corrected the calculation deviation of the traditional interpolation method at singular points. Through the piecewise linear streamline tracing algorithm and time-of-flight conversion, the three-dimensional convection problem was transformed into a one-dimensional equation solution, which significantly improved the calculation efficiency and accuracy.

[0018] (3) It can effectively suppress numerical dissipation under various reservoir conditions, significantly reduce the width of the transition zone at the water drive front, and make the saturation distribution more consistent with reality. At the same time, it can accurately capture the streamline changes in heterogeneous reservoirs, adapt to complex boundary morphology, and has natural parallelism, with computational efficiency significantly better than traditional methods.

[0019] 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

[0020] Figure 1 This is a flowchart of a meshless streamline simulation method for oil reservoirs according to the present invention; Figure 2 This invention node A schematic diagram showing the influence area and the nodes it contains; Figure 3 This invention includes nodes for injection and production wells. j Example image of a local point cloud of the Cartesian shape; Figure 4 This is a schematic diagram of streamline tracing in a meshless point cloud according to the present invention; Figure 5 This is a schematic diagram of the reservoir computational domain in Embodiment 1 of the present invention; Figure 6 These are two point cloud diagrams used in Embodiment 1 of the present invention; Figure 7 The diagram shows the streamline distribution and flight time distribution obtained by tracking Cartesian point clouds and irregular point clouds respectively using the meshless streamline simulation method of Embodiment 1 of the present invention. Figure 8 This is a comparison chart of the oil saturation distribution on the streamline at 100 days, 200 days, and 300 days, calculated under different conditions according to Embodiment 1 of the present invention. Figure 9 This is a comparison chart of the average oil saturation distribution of the node control domain / grid calculated under different conditions in Embodiment 1 of the present invention over 200 days; Figure 10 This is a comparison chart of the average oil saturation distribution of the node control domain / grid calculated under different conditions in Embodiment 1 of the present invention over 300 days; Figure 11 This is a permeability distribution diagram of a heterogeneous material in Embodiment 2 of the present invention; Figure 12 This is a diagram showing streamlines and their flight time distribution under different point cloud conditions, as shown in Embodiment 2 of the present invention. Figure 13 This is a map showing the distribution of oil saturation on streamlines calculated under different point cloud conditions in Embodiment 2 of the present invention; Figure 14 This is a distribution map of the average oil saturation of the node control area calculated under different point cloud conditions in Embodiment 2 of the present invention; Figure 15 This is a schematic diagram of the reservoir computational domain and point cloud in Embodiment 3 of the present invention; Figure 16 The streamlines calculated by EFVM-SL in Embodiment 3 of the present invention and the time-of-flight distribution diagram thereon are shown. Figure 17 This is a diagram showing the oil saturation distribution on the streamline at 100 days, 300 days, and 500 days, calculated by EFVM-SL in Embodiment 3 of the present invention. Figure 18 The figures show the average oil saturation distribution of the node control region calculated by EFVM-SL and the average oil saturation distribution of the grid calculated by EFVM in Embodiment 3 of the present invention. Detailed Implementation

[0021] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0022] Example like Figure 1 As shown, this invention provides a meshless streamline simulation method for oil reservoirs, comprising the following steps: S1. Generate a meshless point cloud within the reservoir computational domain, and calculate the control volume of each node in the point cloud based on the extended finite volume method.

[0023] The Extended Finite Volume Method (EFVM) originates from the Generalized Finite Difference Method (GFDM), a meshless difference method based on Taylor expansion and weighted least squares. In this method, each node possesses a local point cloud composed of nearby nodes that participated in constructing the discrete scheme of the unknown function derivative of that node.

[0024] S11. Discretize the reservoir computational domain using a meshless method to obtain the node set. This forms an initial gridless point cloud, in which, This represents the total number of nodes.

[0025] S12 represents each node in the point cloud. The local point cloud is determined by using a circular influence domain. and neighboring node set At the node When constructing the generalized difference operator, except for the nodes Besides itself, the set of nodes participating in the construction is defined as a node. The set of adjacent nodes, denoted as ,therefore, .

[0026] like Figure 2 As shown, let Includes Each node, in the enumeration When including nodes in the global point cloud, there is a specific order, and the corresponding nodes are in the global point cloud. The serial number in the middle is denoted as , At this time, if Includes nodes , record nodes exist The serial number is Then there should be .

[0027] Based on the generalized finite difference method, each node is defined using the extended finite volume method. A dedicated control area; ; in, For nodes The control domain, For nodes Controlling volume, Let be the volume of the computational domain.

[0028] S13. Based on the computational domain point cloud, EFVM further defines a connectable point cloud; for nodes... Determine whether it is a node Connectable point clouds, if they satisfy This constructs a connectable point cloud. Given this connectable point cloud, EFVM, based on the divergence theorem and the two-point flux estimation scheme, pairs each node with... Construct a linear equation concerning the control volume; ; in, For nodes Controlling volume, , , , For nodes and nodes The relevant discrete coefficients of the generalized finite difference method; It is important to note that when and The above formula still holds true when the thickness is different.

[0029] S14, if node Add virtual nodes outside the computation domain to serve as boundary nodes. This is done to improve the quality of the local point cloud at boundary nodes, making the centroid of the local point cloud closer to the node itself. Typically, a virtual node is added to each boundary node; the control volume of the boundary node is then adjusted based on these virtual nodes. Assume there are boundary nodes The coordinates of its corresponding virtual node are: ; in, and They are nodes and nodes coordinates Boundary nodes The unit outward normal vector at that location, For nodes Local point cloud Node-to-node Weighted average distance, i.e. , For nodes Relative to node The weight function value; It should also be noted that after adding virtual nodes outside the boundary nodes, the control volume corresponding to the boundary nodes will no longer be completely contained within the computational domain. In this case, the corrected control volume formula is: ; S15. In the EFVM framework, the control volume of a node on the boundary that has undergone virtual point processing is called the node virtual control volume, to distinguish it from the node's real control volume that is entirely within the computational domain; this is based on the node characteristic angle of the boundary node. The control volume equations of the boundary nodes are modified, and the control volume equations of all node pairs and the boundary nodes are combined to form an overdetermined linear system of equations. ; in, To calculate the total number of node pairs that can be formed within the domain, These are weighting coefficients. , ; Reasoning process: If it exists right Node pairs, then each pair of nodes corresponds to a node such as The equation shown. Combined with... ,available There are 1 linear equation. The total number of unknowns in the nodal control volume is 1. And usually satisfy > Therefore, this constitutes an overdetermined system of linear equations, which can be solved to obtain... Then calculate To overcome the numerical problems caused by uneven node distribution, an empirical method based on weighted least squares is adopted. This method obtains high-precision control volume values ​​by establishing and solving an overdetermined system of equations concerning the node control volume. Step 1: Select an appropriate method to determine the local point cloud for each node, i.e., determine the point cloud within the local point cloud. Node pairs are used to obtain the connection point cloud of the computational domain.

[0030] Step 2: Select the weight function to obtain each The governing volume equations for node pairs are derived, and a new equation is obtained by weighting the equations. and This constitutes an overdetermined system of linear equations.

[0031] S16. Solve the overdetermined linear equations using the least squares method to obtain the control volume value for each node.

[0032] Therefore, EFVM does not require describing the geometry of the node control domain, but only focuses on the numerical size of the node control volume.

[0033] S2. Based on the control volume, the extended finite volume method is used to discretize the two-phase flow equation of porous media and solve for the nodal pressure distribution in the meshless point cloud.

[0034] Two-phase flow equations for porous media For porous media, immiscible two-phase flow, phase The mass conservation equation is: ; in, For the sake of the prime minister seepage velocity; As a source-sink term, it is usually caused by injection-production wells in reservoir simulation; Rock porosity; For the sake of the prime minister saturation; For time.

[0035] If Darcy's law is satisfied and the effect of gravity is ignored, then: ; in, For the sake of the prime minister The pressure; For the sake of the prime minister The flow rate; For the sake of the prime minister The relative permeability is the saturation. The function; For the sake of the prime minister The viscosity of a substance is usually a function of pressure, but in practical applications it is often treated as a constant that does not change with pressure. For penetration rate; This represents the pressure gradient.

[0036] Will Substitution phase From the mass conservation equation, we can obtain: ; Ignoring capillary force, let's record By superimposing the oil phase and water phase equations, we obtain: ; in, .

[0037] S21. Obtain the pressure field using the generalized finite difference method. exist First spatial derivative at , and second-order spatial derivative , , ; In the two-dimensional case, let nodes... The coordinates are , A node in The coordinates are . exist The Taylor expansion at point is: ; in, , remember , , , , Define the weighted error function for: ; in, , It is a weight function exist The value at that location. In GFDM, it is typically used for each node. Define an influence region, denoted as . The extent of the influence domain determines which nodes will be included in the local point cloud. A circular influence domain is the most commonly used, and its radius is denoted as . Weight function Quartic spline functions or inverse cubic functions are often chosen. The weight function has two characteristics compared to the weight function: (1) Once the nodes are determined... i A local point cloud, the relative weight between any two points in the point cloud and the radius of the influence domain. (2) The classic nine-point finite difference scheme can be derived through the weight function.

[0038] The formula for calculating the quartic spline function is as follows: ; The formula for calculating the inverse cubic function is: ; in, .

[0039] When the weighted error function To obtain the minimum value, it is required that it is equal to... The partial derivatives of all components are zero, and after simplification, we can obtain the following formula: The linear equation that must be satisfied to obtain the minimum value is: .

[0040] in, , , , , ,

[0041] Solving the above linear equation, we can obtain: ; in, .

[0042] remember The elements are , exist The generalized finite difference expressions for the first and second spatial derivatives at point A are: ; ; ; ; ; in, , , , , For discrete coefficients, For nodes Pressure value, As the central node Pressure value; S22. Integrating the two-phase flow equations for porous media over the nodal control volumes and discretizing them in time using an implicit scheme, the specific formulas obtained are as follows: ; in, For divergence operators, For oil phase fluidity, For water phase mobility, For absolute penetration rate, For pressure gradient, For nodes The spatial region corresponding to the control volume is the node. With nodes The harmonic average of absolute permeability. For time steps At that time, node With nodes The relative permeability of the oil phase between them For time steps At that time, node With nodes The relative permeability of the water phase between them For time steps At that time, node With nodes The viscosity of the oil phase between them For time steps At that time, node With nodes The viscosity of the aqueous phase between them, For time steps At that time, node Pressure value, For time steps At that time, node Pressure value; S23. The phase permeability and viscosity are calculated using the single-point upstream weighted formula, and the absolute permeability is calculated using the harmonic average formula. The calculation formulas are as follows: ; ; ; S24. Extracting internode conductivity using an extended finite volume discretization scheme based on the two-phase flow equation in porous media. And process source and sink items ; ; ; In the reservoir computational domain, source and sink terms induced by injection and production wells are considered singular source and sink terms. Without the concept of node control volume, it would be difficult to handle such singular source and sink terms, which is one of the key motivations for developing EFVM based on GFDM. In EFVM, the integral value of the source and sink term over the node control volume is the injection or production volume flow rate of the corresponding injection or production well at that node, i.e. Furthermore, the source and sink terms related to well injection and production rates are handled in a manner largely consistent with the traditional finite volume method (FVM), but with one key difference: Since EFVM is a meshless method, although the nodal control volumes are calculated, the specific geometry of the nodal control domains remains unknown. Therefore, when calculating well indices, only the nodal control volumes can be utilized. Consequently, the extended finite volume method discretization scheme for the two-phase flow equations in porous media is as follows: ; To directly utilize the nonlinear solver in the existing FVM-based reservoir numerical simulator, conductivity can be introduced to simplify the process as follows: ; in, For time steps At that time, node Source and sink items in the oil phase, For time steps At that time, node Sources and sinks of water phases, This represents the overall compression ratio.

[0043] S25. Based on the extended finite volume method discretization scheme, the closed boundary is processed. In the pressure discretization equation of the boundary node, the conductivity term between it and the virtual node is ignored, and the node control volume is corrected to the real control volume. Discrete equations are established for all real nodes to form a closed linear equation system. The pressure distribution of all nodes is obtained by solving the system.

[0044] S3. Calculate the node seepage velocity based on the node pressure distribution. For nodes containing injection and production wells, use the radial seepage velocity model. For ordinary nodes, use the generalized finite difference method to calculate the pressure gradient to obtain the seepage velocity.

[0045] S31. Divide the nodes into ordinary nodes that do not contain injection and production wells. Sets of source and sink nodes that include source and sink items such as injection and production wells. ; S32, For ordinary nodes The first spatial derivative of the pressure field and the nodal parameters obtained based on the generalized finite difference method The pressure value reaches the node The first spatial derivative of the pressure function is used to estimate the seepage velocity at the node using Darcy's law, combined with the total mobility distribution. The calculation formula is as follows: ; Streamline simulation is most commonly used in scenarios where reservoir pressure distribution does not change significantly over time. In such scenarios, the total mobility distribution of the oil-water two-phase system does not change significantly with changes in saturation; otherwise, as the water injection process continues, significant changes in the total mobility of the oil-water two-phase system will lead to substantial changes in the pressure distribution. Therefore, nodes... seepage velocity at the location It can be estimated as follows: ; in, Local point cloud scale of nodes , for directional seepage velocity components, for directional seepage velocity components, For nodes absolute penetration rate; S33, Node containing injection and production wells Descartes' local point cloud, such as Figure 3 As shown, the pressure values ​​of nodes 1, 2, 3, and 4 are the same, and the pressure values ​​of nodes 5, 6, 7, and 8 are also the same. Therefore, for the source-sink node... By controlling its volume to be equivalent to a cylinder, its equivalent radius is calculated. The radial seepage velocity value at this node is obtained. This allows for the calculation of the seepage velocity at the nodes. The calculation formula is as follows: ; ; ; in, To calculate the coordinates of any node to be interpolated within the domain, For nodes The coordinates; S34. Integrate the seepage velocities of all ordinary nodes and source / sink nodes, and obtain the seepage velocity at any point in the computational domain through linear interpolation. This forms a continuous seepage velocity field within the reservoir computational domain, calculated using the following formula: ; in, For nodes At the point to be interpolated Linear interpolation basis functions, For nodes At the point to be interpolated Linear interpolation basis functions.

[0046] The specific process of linear interpolation is as follows: a simplex (a triangle in two dimensions and a tetrahedron in three dimensions) composed of neighboring data points is constructed around each interpolation point. Linear fitting is then performed within this simplex to calculate the seepage velocity value at the interpolation point. This interpolation method is fast and does not produce values ​​outside the range of known data.

[0047] S4. Perform piecewise linear streamline tracing based on seepage velocity in a gridless point cloud to determine the streamline trajectory and flight time along the streamline.

[0048] S41. Taking the injection well node as the center, and... On a circle with radius [radius value], distribute the starting points of the streamlines at equal intervals according to the required number of streamlines. Of course... It can also be modified to Other values, such as these, usually have little impact on the final simulation results.

[0049] S42, such as Figure 4 As shown, the starting point of the streamline is defined as denoted as . Calculate its seepage velocity Find the node closest to the starting point of the streamline. Obtain the porosity of the node. Set a small distance Calculate flight time using the single tracking distance step size. The streamline position after the movement is obtained. ; ; ; S43, with Repeat step S42 to obtain new nodes. and corresponding flight time The process continues until the streamline reaches the node where the production well is located, completing the single streamline tracing. ; .

[0050] S5. Based on the definition of flight time, the two-dimensional / three-dimensional convection transport problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate. The water saturation equation is solved along each streamline to obtain the saturation distribution on the streamline.

[0051] S51. Transform the two-dimensional / three-dimensional convection transport problem into a one-dimensional flow branch function problem through streamline simulation; ; in, Porosity Water saturation Moisture content, For total phase velocity, The distance along the streamline; S52. Based on the definition of flight time, the above one-dimensional split function problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate.

[0052] When streamlines pass through different grids, the velocity on the streamlines This usually changes. Under heterogeneous conditions, the porosity of the mesh through which each streamline passes will also differ. Therefore, in the above formula... and All are about The complexity of the function makes it difficult to solve the equation directly on the streamline. The concept of Time of Flight (TOF) is defined, which represents the distance traveled along the streamline from the starting point at the actual speed. / distance traveled Time required . Simultaneously apply the integral expression to both sides Taking the derivative, we get Substituting this into the above one-dimensional split-function problem, we can transform it into a one-dimensional constant-coefficient equation with flight time as the coordinate: Through the streamline tracing process described above, the flight time of each streamline as it passes through each grid can be calculated. It can be observed that... It is a The constant coefficient pure convection equation as spatial coordinates.

[0053] S53. Solve the above one-dimensional equation using numerical methods to obtain the water saturation at each point on the streamline. Numerical methods can include the windward finite difference method, the discontinuous Galerkin method, the weighted essential non-oscillatory (WENO) method, etc. ; in, This refers to the local time step when calculating water saturation along the streamline. For flight time step; and CFL conditions must be met. It is worth noting that... The maximum value depends only on the selected The streamline simulation method is not limited by the grid size and velocity values ​​within the reservoir computational domain. In contrast, the traditional IMPES method imposes strict constraints on the time step, especially in fractured reservoirs where the small size and high velocity of fracture elements make these constraints more pronounced. Therefore, the streamline simulation method effectively addresses this issue. Furthermore, the tracing process of each streamline is independent of the calculation of saturation distribution along the streamline, making it naturally suitable for parallel computing and thus achieving higher computational efficiency.

[0054] S6. Based on the relationship between streamline segments and node control areas, map streamline saturation to node control volumes and calculate the average water saturation within the node control areas.

[0055] In traditional grid-based streamline tracing, the average water saturation of all streamlines passing through a grid segment is typically calculated by weighting the flight time required for that segment and then averaging the average water saturation of that grid segment. In the meshless framework presented in this paper, since each streamline segment is traced to the node closest to its starting point, it can be considered to have passed through the control region of that node. Based on this, a similar approach to traditional methods can be used to calculate the average saturation within the control regions of each node. However, there are some differences: in the meshless streamline simulation method proposed in this paper, the same streamline may require multiple consecutive tracings when passing through the control region of a node, depending on the distance step size of each tracing. The selection of step size. When the size is smaller, a streamline may require more tracking attempts to pass through a node's control area; however, in grid-based streamline simulation methods, a streamline typically only needs to be tracked once within a grid.

[0056] Therefore, node i The average water saturation within the control area can be calculated as follows: ; in, For the nodes that flow through The streamline set of the control area For set The first in A streamline.

[0057] Example 1 like Figure 5 As shown, this embodiment uses a reservoir model with dimensions of 200m × 200m × 10m. A water injection well is located at the center of the model, with a fixed water injection rate of 80m / s. 3 One production well is set up at each of the four corners of the reservoir per day, with a constant fluid production rate of 20m / day. 3 / day. Table 1 lists the core physical properties of reservoir rocks and fluids. To visually compare the differences in numerical dissipation between streamline simulation and traditional methods in saturation calculation, this example uses special parameter settings. First, the viscosity of both oil and water is set to 1 cp. Second, a linear phase permeability curve is used where "the sum of oil phase permeability and water phase permeability is always 1". This setting makes the analytical solution of the water drive front exhibit a pattern "from 1- S or (Residual oil saturation) to The discontinuity of "(bound water saturation)" means that "1-" in the saturation distribution results S or to The width of the transition band can directly reflect the magnitude of the numerical dissipation error.

[0058] Table 1 Some reservoir rock and fluid physical properties

[0059] To verify the adaptability of the streamline simulation method in this embodiment to different types of point clouds, Embodiment 1 employs... Figure 6 (a) Cartesian point cloud and Figure 6 (b) The irregular point cloud is used for calculation. Figure 7 (a) and (b) respectively show the streamline distribution and time-of-flight distribution of the Cartesian point cloud. Figure 7 (c) and (d) respectively show the streamlines and flight time distribution of the irregular point cloud.

[0060] Figure 8 The oil saturation distribution on the streamlines of this embodiment, EFVM-SL, was further compared at 100, 200, and 300 days. To verify the calculation accuracy of this method, Figure 8 (g), (h), and (i) simultaneously present the results of the classical streamline simulation method (i.e., FVM-SL, because the classical streamline simulation method generally uses FVM to calculate pressure and seepage velocity distribution). This method is based on a 10m×10m×10m Cartesian grid and uses the Pollock streamline tracing algorithm. The comparison shows that, regardless of whether Cartesian point clouds or irregular point clouds are used, the calculation results of the meshless streamline simulation method are in high agreement with the classical streamline simulation method.

[0061] Figure 9 , Figure 10 The oil saturation distribution at 200 and 300 days was compared under three computational scenarios (two types of meshless point clouds + Cartesian mesh). Specifically, this included the nodal control domain / mesh-average oil saturation calculated using the meshless Extended Finite Volume Method (EFVM) and the Finite Volume Method (FVM), as well as the nodal control domain / mesh-average oil saturation calculated using the meshless streamline simulation method and the classical streamline simulation method. Figure 9 ,Figure 10 The comparison results clearly show that the meshless streamline simulation method can significantly reduce the numerical dissipation error of the meshless extended finite volume method (EFVM). Table 2 compares the CPU time of EFVM-SL and EFVM under two point cloud conditions, using an Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz (2.59GHz) as the computer processor. It can be seen that due to the inherent parallelism of the streamline simulation method, when simulating 1000 days of reservoir injection-production dynamics, the computation time of EFVM-SL is less than one-fifth of that of EFVM, demonstrating significantly higher computational efficiency.

[0062] Table 2 Comparison of CPU time for different methods in Example 1

[0063] Example 2 This second example uses the same computational domain as the first example, but employs... Figure 11 The heterogeneous permeability distribution shown is used to verify the adaptability of the proposed meshless streamline simulation method EFVM-SL to heterogeneous reservoir models. Streamline simulation calculations for Example 2 are performed using Cartesian point clouds and irregular point clouds from Example 1, respectively. Figure 12 The streamlines calculated for these two point cloud scenarios and their time-of-flight distributions are shown. Figure 13 The distribution of oil saturation on streamlines at 100 days, 300 days, and 500 days was compared under these two point cloud conditions. Figure 14 The average oil saturation of the nodal control area at 100 days, 300 days, and 500 days was compared under these two point cloud conditions. It can be seen that the flight time distribution on streamlines, the oil saturation distribution on streamlines, and the average oil saturation distribution of the nodal control area calculated using these two point cloud conditions are very similar. This indicates that EFVM-SL can be effectively applied to heterogeneous reservoir models under various point cloud conditions. It can also be observed that compared to the homogeneous reservoir model in Example 1, the heterogeneous permeability distribution in Example 2 alters the velocity distribution of the flow field, thus distorting the streamline distribution.

[0064] Example 3 like Figure 15 As shown in (a), this third embodiment is a reservoir model with irregular boundaries. There is a water injection well at location (225, 50) with a fixed injection rate of 60 cubic meters per day. There are also two production wells at locations (0, 60) and (460, 160), both with a fixed extraction rate of 30 cubic meters per day. Other physical properties are the same as in Example 1. This third embodiment uses... Figure 15 The point cloud in (b) is used for calculation. Figure 16The streamline distribution and time-of-flight distribution along the streamlines calculated by EFVM-SL are shown. Figure 17 The distribution of oil saturation on the streamline at 100 days, 300 days, and 500 days, calculated by EFVM-SL, is shown. Figure 18 The average oil saturation distribution of the node control region calculated by EFVM-SL and the average oil saturation distribution of the grid calculated by EFVM were compared. It can be seen that the EFVM-SL proposed in this embodiment can significantly reduce the numerical dissipation error of EFVM, obtain high-precision saturation distribution calculation results, and is applicable to reservoir models with complex morphology.

[0065] Therefore, this invention adopts the above-mentioned meshless streamline simulation method for reservoirs, which for the first time combines meshless and streamline simulation, effectively overcoming the limitations of traditional mesh methods and the numerical dissipation of meshless methods; through innovative velocity calculation and streamline tracing, it significantly improves the accuracy of saturation calculation and simulation efficiency, and shows excellent adaptability to complex reservoir conditions.

[0066] 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 meshless streamline simulation method for oil reservoirs, characterized in that, Includes the following steps: S1. Generate a meshless point cloud in the reservoir computational domain, and calculate the control volume of each node in the point cloud based on the extended finite volume method. S2. Based on the control volume, the extended finite volume method is used to discretize the two-phase flow equation of porous media and solve for the nodal pressure distribution in the meshless point cloud. S3. Calculate the node seepage velocity based on the node pressure distribution. For nodes containing injection and production wells, use the radial seepage velocity model. For ordinary nodes, use the generalized finite difference method to calculate the pressure gradient to obtain the seepage velocity. S4. Perform piecewise linear streamline tracing based on seepage velocity in a gridless point cloud to determine the streamline trajectory and flight time along the streamline. S5. Based on the definition of flight time, the two-dimensional / three-dimensional convection transport problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate. The water saturation equation is solved along each streamline to obtain the saturation distribution on the streamline. S6. Based on the relationship between streamline segments and node control areas, map streamline saturation to node control volumes and calculate the average water saturation within the node control areas.

2. The method for simulating reservoir streamlines without a mesh as described in claim 1, characterized in that, The specific steps of S1 include: S11. Discretize the reservoir computational domain using a meshless method to obtain the node set. This forms an initial gridless point cloud, in which, The total number of nodes; S12 represents each node in the point cloud. The local point cloud is determined by using a circular influence domain. and neighboring node set ,in, Based on the generalized finite difference method, each node is defined using the extended finite volume method. A dedicated control area; ； in, For nodes The control domain, For nodes Controlling volume, The volume of the computational domain; S13, For nodes Determine whether it is a node Connectable point clouds, if they satisfy This constructs a connectable point cloud. Based on the connectable point cloud, a pair is created for each node. Construct a linear equation concerning the control volume; ； in, For nodes Controlling volume, , , , For nodes and nodes The relevant discrete coefficients of the generalized finite difference method; S14, if node Add virtual nodes outside the computation domain to serve as boundary nodes. The control volume of boundary nodes is corrected based on virtual nodes; The coordinates of the virtual node are: ; in, and They are nodes and nodes coordinates Boundary nodes The unit outward normal vector at that location, For nodes Local point cloud Node-to-node Weighted average distance, i.e. , For nodes Relative to node The weight function value; The revised control volume formula is: ; S15. Based on the nodal characteristic angles of the boundary nodes. The control volume equations of the boundary nodes are modified, and the control volume equations of all node pairs and the boundary nodes are combined to form an overdetermined linear system of equations. ； in, To calculate the total number of node pairs that can be formed within the domain, These are weighting coefficients. , ; S16. Solve the overdetermined linear equations using the least squares method to obtain the control volume value for each node.

3. The method for simulating reservoir streamlines without a mesh as described in claim 2, characterized in that, The specific steps of S2 include: S21. Obtain the pressure field using the generalized finite difference method. exist First spatial derivative at , and second-order spatial derivative , , ; ； ； ； ； ； in, , , , , For discrete coefficients, For nodes Pressure value, As the central node Pressure value; S22. Integrating the two-phase flow equations for porous media over the nodal control volumes and discretizing them in time using an implicit scheme, the specific formulas obtained are as follows: ； in, For divergence operators, For oil phase fluidity, For water phase mobility, For absolute penetration rate, For pressure gradient, For nodes The spatial region corresponding to the control volume is the node. With nodes The harmonic average of absolute permeability. For time step At that time, node With nodes The relative permeability of the oil phase between them For time step At that time, node With nodes The relative permeability of the water phase between them For time step At that time, node With nodes The viscosity of the oil phase between them For time step At that time, node With nodes The viscosity of the aqueous phase between them, For time step At that time, node Pressure value, For time step At that time, node Pressure value; S23. The phase permeability and viscosity are calculated using the single-point upstream weighted formula, and the absolute permeability is calculated using the harmonic average formula. The calculation formulas are as follows: ； ； ； S24. Extracting internode conductivity using an extended finite volume discretization scheme based on the two-phase flow equation in porous media. And process source and sink items ; ； ； S25. Based on the extended finite volume method discretization scheme, the closed boundary is processed. In the pressure discretization equation of the boundary node, the conductivity term between it and the virtual node is ignored, and the node control volume is corrected to the real control volume. Discrete equations are established for all real nodes to form a closed linear equation system. The pressure distribution of all nodes is obtained by solving the system.

4. The reservoir streamline simulation method without meshes according to claim 3, characterized in that, The extended finite volume method discretization scheme for the two-phase flow equation in porous media in S24 is as follows: ； Introducing conductivity simplifies to: ； in, For time step At that time, node Source and sink items in the oil phase, For time step At that time, node Sources and sinks of water phases, This represents the overall compression factor.

5. The method for simulating reservoir streamlines without a mesh as described in claim 3, characterized in that, The specific steps of S3 include: S31. Divide the nodes into ordinary nodes that do not contain injection and production wells. Sets of source and sink nodes that include source and sink items such as injection and production wells. ; S32, For ordinary nodes The first spatial derivative of the pressure field and the nodal parameters obtained based on the generalized finite difference method The pressure value reaches the node The first spatial derivative of the pressure function is used to estimate the seepage velocity at the node using Darcy's law, combined with the total mobility distribution. The calculation formula is as follows: ； ； in, Local point cloud scale of nodes , for directional seepage velocity components, for directional seepage velocity components, For nodes absolute penetration rate; S33, For source and sink nodes By controlling its volume to be equivalent to a cylinder, its equivalent radius is calculated. The radial seepage velocity value at this node is obtained. This allows for the calculation of the seepage velocity at the nodes. The calculation formula is as follows: ； ； ； in, To calculate the coordinates of any node to be interpolated within the domain, For nodes The coordinates; S34. Integrate the seepage velocities of all ordinary nodes and source / sink nodes, and obtain the seepage velocity at any point in the computational domain through linear interpolation. This forms a continuous seepage velocity field within the reservoir computational domain, calculated using the following formula: ； in, For nodes At the point to be interpolated Linear interpolation basis functions, For nodes At the point to be interpolated Linear interpolation basis functions.

6. The method for simulating reservoir streamlines without a mesh as described in claim 1, characterized in that, The specific steps of S4 include: S41. Taking the injection well node as the center, and... On a circle with radius , the starting points of the streamlines are distributed at equal intervals according to the required number of streamlines; S42. Determine the starting point of the streamline as (Note: The original text contains a typo and can be left as a typo). Calculate its seepage velocity Find the node closest to the starting point of the streamline. Obtain the porosity of the node. Set a small distance Calculate flight time using the single tracking distance step size. The streamline position after the movement is obtained. ; ; ; S43, with Repeat step S42 to obtain new nodes. and corresponding flight time The process continues until the streamline reaches the node where the production well is located, completing the single streamline tracing. ； 。 7. The method for simulating reservoir streamlines without a mesh as described in claim 6, characterized in that, The specific steps of S5 include: S51. Transform the two-dimensional / three-dimensional convection transport problem into a one-dimensional flow branch function problem through streamline simulation; ； in, Porosity This represents the water saturation level. Moisture content, For total phase velocity, The distance along the streamline; S52. Based on the definition of flight time, the above one-dimensional split function problem is transformed into a one-dimensional constant coefficient equation with flight time as the coordinate. Flight time is defined as: ; The rewritten one-dimensional constant coefficient equation is: ; S53. Solve the above one-dimensional equation using numerical methods to obtain the water saturation at each point on the streamline. ； in, This refers to the local time step when calculating water saturation along the streamline. This represents the flight time step.

8. The method for simulating reservoir streamlines without a mesh as described in claim 7, characterized in that, The formula for calculating the average water saturation within the node control area in S6 is as follows: ； in, For the nodes that flow through The streamline set of the control area For set The first in A streamline.

Images