EFVM-based meshless numerical simulation method for wellbore-reservoir coupling
The EFVM method is used to construct a meshless numerical simulation of wellbore-reservoir coupling, which solves the problems of insufficient prediction accuracy and complex mesh generation caused by wellbore simplification in traditional methods. It realizes a dynamic feedback mechanism between wellbore and reservoir, improves simulation accuracy and robustness, and is suitable for efficient simulation of multi-layered strongly heterogeneous reservoirs.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- YANGTZE UNIVERSITY
- Filing Date
- 2026-02-28
- Publication Date
- 2026-06-02
AI Technical Summary
Traditional reservoir numerical simulation methods simplify the wellbore to constant pressure or constant flow boundary conditions, ignoring the complex multiphase flow processes within the wellbore, resulting in insufficient prediction accuracy. Furthermore, existing well-reservoir coupling methods suffer from problems such as complex grid generation, poor geometric adaptability, and unstable coupling convergence, making it difficult to meet the simulation requirements of multi-layered, highly heterogeneous reservoirs and complex well networks.
A meshless numerical simulation method based on EFVM for wellbore-reservoir coupling is adopted. By constructing the reservoir computational domain and discretizing the point cloud, the simulation parameters are initialized. The meshless extended finite volume method (EFVM) is used to simulate reservoir seepage. Combined with Beggs-Brill wellbore multiphase flow calculation, iterative convergence and relaxation update of the layered bottom hole flowing pressure are realized, and a dynamic feedback mechanism between wellbore multiphase flow and reservoir seepage is established.
It improves the accuracy of pressure and production prediction, enhances geometric adaptability, can flexibly handle irregular boundary reservoirs and complex well networks, ensures robust convergence under complex conditions, and is suitable for simulating multi-layered, highly heterogeneous reservoirs.
Smart Images
Figure CN122133394A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of reservoir numerical simulation technology, and in particular to a meshless numerical simulation method for wellbore-reservoir coupling based on EFVM. Background Technology
[0002] Reservoir numerical simulation is a core technology for dynamic analysis and production prediction in oil and gas development, and an important means to achieve refined reservoir management and efficient development. By analyzing fluid seepage, phase changes, and reservoir pressure variations, it provides a scientific basis for reservoir productivity evaluation, scheme optimization, and long-term development decisions.
[0003] Traditional reservoir numerical simulations simplify the wellbore to constant pressure or constant flow boundary conditions, neglecting the complex multiphase flow processes within the wellbore and their dynamic feedback to formation pressure and production, resulting in insufficient prediction accuracy. Existing well-reservoir coupling methods mostly rely on grid discretization techniques, which suffer from problems such as complex grid generation, poor geometric adaptability, and unstable coupling convergence, making it difficult to meet the simulation requirements of multi-layered, highly heterogeneous reservoirs and complex well networks. Summary of the Invention
[0004] The purpose of this invention is to provide a meshless numerical simulation method for wellbore-reservoir coupling based on EFVM, which solves the problems of grid dependence, low coupling accuracy, and poor convergence in existing technologies.
[0005] To achieve the above objectives, this invention provides a meshless numerical simulation method for wellbore-reservoir coupling based on EFVM, comprising the following steps: S1. Construct the reservoir computational domain and discretize the point cloud; S2. Initialize simulation parameters; S3, Meshless Extended Finite Volume Method (EFVM) reservoir seepage simulation; S4, Beggs-Brill wellbore multiphase flow calculation; S5. Coupled iterative convergence and relaxation update; S6, Phase Advancement and Simulation Termination.
[0006] Preferably, in S1, a gridless point cloud discrete computational domain is generated based on the actual geological parameters of the reservoir. The point cloud contains a total of The point cloud consists of nodes, including internal nodes and boundary nodes of the reservoir. Virtual nodes are added to the boundary nodes to optimize the local point cloud quality. The reservoir includes regular rectangular reservoirs, irregular multi-layered reservoirs, and actual oilfield blocks. The actual geological parameters of the reservoir include reservoir geometry, layered permeability, porosity, and layer thickness.
[0007] Preferably, the virtual node coordinates in S1 are calculated using the boundary node coordinates, the weighted average distance between neighboring nodes, and the boundary unit outward normal vector, as expressed in the following expression: ; In the formula, Represents virtual nodes The coordinates; Represents boundary nodes The coordinates; Represents boundary nodes The unit outward normal vector at that location; Represents the distance between neighboring nodes. The weighted average distance.
[0008] Preferably, in S2, the total simulation time, time stage division, wellhead pressure-time function, wellbore parameters, fluid parameters, and convergence threshold are set, and the predicted value of the layered bottom hole flowing pressure for the first time stage is set. The predicted value of the stratified bottom-hole flowing pressure at each subsequent time stage is the stratified bottom-hole flowing pressure after convergence in the previous time stage; the wellhead pressure time function includes the production wellhead pressure time function and the injection wellhead pressure time function, the production wellhead pressure decreases with time, and the injection wellhead pressure increases with time; the wellbore parameters include the wellbore diameter, wellbore length and wellbore roughness; the fluid parameters include the fluid viscosity, fluid density and fluid phase characteristics.
[0009] Preferably, in S3, based on the point cloud of S1, a local support domain and equivalent control volume are generated for each node. The reservoir two-phase control equations are discretized using the EFVM method, and the predicted values of the bottom-hole flowing pressure of the stratified well are obtained. As boundary conditions, the reservoir pressure field and saturation field are solved using Newton-Raphson iteration to output the total flow rate of each well. and tiered traffic .
[0010] Preferably, the specific process of discretizing the reservoir two-phase control equations using the EFVM method in S3 is as follows: A nodal generalized difference operator is constructed using Taylor expansion and weighted least squares method. Combined with the nodal control volume and two-point flux scheme, a discrete equation satisfying local mass conservation is obtained. The expression of the discrete equation is: ; In the formula, and Both represent node numbers. , ,and ; Represents a node The number of nodes in a local point cloud; Indicates in Time, node With nodes The relative permeability of oil or water between them; Indicates in Time, node With nodes The volume coefficient between oil and water; Indicates in Time, node With nodes The viscosity between oil and water; Indicates in Time, node With nodes Transfer rate between; Indicates in Time, node With nodes The pressure between oil or water; Indicates in Time, node Pressure of oil or water; Represents a node Source / collection; Indicates the volume controlled by the node; Indicates the time step; Indicates in Time, node Porosity; Indicates in Time, node Porosity; Indicates in Time, node The saturation level of oil or water; Indicates in Time, node The saturation level of oil or water.
[0011] Preferably, in S4, the wellhead pressure of S2 and the total flow rate of S3 are used as inputs. Based on the wellbore geometry and fluid properties, a segmented advancement and intra-segment iteration solution strategy is adopted to sequentially perform flow pattern determination, frictional pressure drop calculation, gravity term calculation, and acceleration term calculation, thereby obtaining the calculated values of the pressure distribution along the wellbore and the stratified bottom hole flowing pressure. The wellbore inclination angle is set differently depending on the well type; the wellbore inclination angle of the production well is taken as... The fluid flows upward from the formation to the wellhead, and the wellbore inclination angle of the water injection well is taken as... Fluid is injected into the formation from the wellhead; well types include vertical wells, deviated wells, horizontal wells, and multi-well combinations.
[0012] Preferably, in S5, the predicted value of the bottom-hole flowing pressure of the stratified well is calculated. Calculated values of bottom-hole flowing pressure in stratified wells The absolute value of the error is calculated and tested independently for injection wells and production wells; if the absolute value of the error is not less than the threshold... If the layered bottom-hole flowing pressure is not found, the coupling update strategy of the layered relaxation factor is used to update the layered bottom-hole flowing pressure, and the new predicted value of the layered bottom-hole flowing pressure is obtained, then return to S3; otherwise, the coupling convergence is successful, and the next step is continued; the update expression is: ; In the formula, Indicates the first The predicted values of the bottom-hole flowing pressure of the stratified well are obtained after the next iteration update; Indicates the first The predicted value of the layered bottom flow pressure set in the next iteration; Represents the global relaxation factor; The layer weight matrix represents the layer weight matrix that satisfies the following conditions: , Indicates the first Layer weight coefficients; This represents the layered bottom flow pressure value calculated by the wellbore multiphase flow model.
[0013] Preferably, in S5, the convergence of the iteration is determined by the dimensionless or absolute change in the bottom hole pressure, expressed as: ; In the formula, Indicates all layers Take the maximum value.
[0014] Preferably, in S6, the EFVM reservoir seepage simulation calculation for the current time stage is completed based on the converged stratified bottom hole flowing pressure; it is determined whether the total simulation time has been reached. If not, the reservoir state (pressure, saturation) is updated as the initial value for the next time stage, and the process returns to S2; if the total simulation time has been reached, the simulation results are output. The simulation results include wellhead pressure, stratified bottom hole pressure, stratified flow rate, total flow rate, pressure distribution along the wellbore, reservoir pressure field, and saturation field.
[0015] Therefore, the present invention employs the above-mentioned meshless numerical simulation method for wellbore-reservoir coupling based on EFVM, which has the following beneficial effects: (1) A strong coupling between the Beggs–Brill model and EFVM was achieved: Through bidirectional iteration of bottom-hole flow pressure and flow rate, a dynamic feedback mechanism between multiphase flow in the wellbore and reservoir seepage was established, which solved the problem of ignoring the dynamic influence of the wellbore in traditional unidirectional coupling and improved the accuracy of pressure and production prediction. (2) Meshless characteristics enhance geometric adaptability: EFVM is based on point cloud discretization, which does not require complex mesh generation. It can flexibly handle irregular boundary reservoirs, fractured reservoirs and complex well networks, reducing the mesh effect and discretization difficulty. (3) Layered weighted relaxation factor ensures convergence stability: For multi-layered heterogeneous and strongly nonlinear problems, the method avoids coupled iterative oscillation and divergence by dual regulation of layer weight and global relaxation factor, thereby improving the robustness of the method under complex conditions.
[0016] 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
[0017] Figure 1 This is an overall flowchart of the wellbore-reservoir coupled meshless numerical simulation method based on EFVM of the present invention; Figure 2 This is a schematic diagram of the reservoir computational domain point cloud in Embodiment 2 of the present invention; where (a) is the reservoir computational domain and (b) is the point cloud; Figure 3 The figures are the pressure variation curves at the wellheads of the production well and the injection well in Embodiment 2 of the present invention; where (a) is the injection well and (b) is the production well. Figure 4 This is a simulated 150-day oil saturation distribution map of an oil reservoir in Example 2 of the present invention; where (a) is the first layer, (b) is the second layer, (c) is the third layer, (d) is the fourth layer, and (e) is the fifth layer; Figure 5 This is a simulated reservoir pressure distribution map for 150 days according to Embodiment 2 of the present invention; where (a) is the first layer, (b) is the second layer, (c) is the third layer, (d) is the fourth layer, and (e) is the fifth layer; Figure 6 The curves showing the total liquid extraction volume of the injection well and the production well as the change over time in Embodiment 2 of the present invention are shown; where (a) is the injection well, (b) is production well 1, and (c) is production well 2. Figure 7 The curves showing the variation of the stratified flow rates of the injection well and the production well over time in Embodiment 2 of the present invention are shown; where (a) is the injection well, (b) is production well 1, and (c) is production well 2. Figure 8 The curves showing the variation of bottom-hole flowing pressure of the injection well and the production well over time in Embodiment 2 of the present invention are shown; where (a) is the injection well, (b) is production well 1, and (c) is production well 2. Figure 9 This is a pressure profile of the injection well and the production well in the 30th time stage of Embodiment 2 of the present invention; wherein (a) is the pressure profile of the injection well, (b) is the pressure profile of production well 1, and (c) is the pressure profile of production well 2. Figure 10 This is a reservoir distribution map of Embodiment 3 of the present invention; wherein (a) is a reservoir connectivity map and well distribution, (b) is a three-dimensional permeability distribution map of a single reservoir, and (c) is a three-dimensional map of the entire reservoir and well group distribution; Figure 11 The pressure change curves at the wellheads of the production well and the injection well in Embodiment 3 of the present invention are shown; (a) is the injection well, and (b) is the production well. Figure 12 This is a distribution map of oil saturation in an oil reservoir according to Example 3 of the present invention; where (a) is the second layer, (b) is the fourth layer, (c) is the sixth layer, and (d) is the eighth layer; Figure 13 This is a reservoir pressure distribution diagram of Embodiment 3 of the present invention; where (a) is the second layer, (b) is the fourth layer, (c) is the sixth layer, and (d) is the eighth layer; Figure 14 The curves showing the total extraction volume of each well in Example 3 of the present invention change over time; where (a) is production well P10, (b) is production well P20, (c) is production well P3, (d) is production well P15, (e) is production well P1, (f) is water injection well P9, (g) is water injection well P23, and (h) is water injection well P11. Figure 15 The curves showing the variation of flow rate of each well over time in Embodiment 3 of the present invention are shown; where (a) is production well P10, (b) is production well P20, (c) is production well P3, (d) is production well P15, (e) is production well P1, (f) is water injection well P9, (g) is water injection well P23, and (h) is water injection well P11. Figure 16 The curves showing the variation of bottom-hole flowing pressure of each well in Embodiment 3 of the present invention over time are shown; where (a) is production well P10, (b) is production well P20, (c) is production well P3, (d) is production well P15, (e) is production well P1, (f) is water injection well P9, (g) is water injection well P23, and (h) is water injection well P11. Detailed Implementation
[0018] 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.
[0019] Example 1 Please see Figure 1 The meshless numerical simulation method for wellbore-reservoir coupling based on EFVM includes the following steps: S1. Construct the reservoir computational domain and discretize the point cloud; based on the reservoir geometry, layer permeability, porosity, and layer thickness, generate a meshless point cloud within the computational domain, covering all areas of the reservoir; for boundary nodes, optimize the symmetry and density of the local point cloud by adding virtual nodes. Calculate the virtual node coordinates using the boundary node coordinates, the weighted average distance to neighboring nodes, and the boundary unit outward normal vector to ensure accurate application of boundary conditions; the expression for the virtual node coordinates is: ; In the formula, Represents virtual nodes The coordinates; Represents boundary nodes The coordinates; Represents boundary nodes The unit outward normal vector at that location; Represents the distance between neighboring nodes. The weighted average distance.
[0020] S2. Initialize simulation parameters; set the total simulation time, time stage division, wellhead pressure-time function, wellbore parameters, fluid parameters, and convergence threshold; set the predicted value of the bottom hole flowing pressure for the first time stage. The predicted value of the stratified bottom-hole flowing pressure at each subsequent time stage is the stratified bottom-hole flowing pressure after convergence in the previous time stage; the wellhead pressure time function includes the production wellhead pressure time function and the injection wellhead pressure time function, the production wellhead pressure decreases with time, and the injection wellhead pressure increases with time; the wellbore parameters include the wellbore diameter, wellbore length and wellbore roughness; the fluid parameters include the fluid viscosity, fluid density and fluid phase characteristics.
[0021] S3, EFVM reservoir seepage simulation; based on the point cloud of S1, a local support domain and equivalent control volume are generated for each node. The reservoir two-phase control equations are discretized using the EFVM method, and the predicted values of the layered bottom hole flowing pressure are obtained. As boundary conditions, the reservoir pressure field and saturation field are solved using Newton-Raphson iteration to output the total flow rate of each well. and tiered traffic The specific process of discretizing the two-phase governing equations of the reservoir using the EFVM method is as follows: Nodal generalized difference operators are constructed using Taylor expansion and weighted least squares method. Combined with the nodal control volume and two-point flux scheme, discrete equations satisfying local mass conservation are obtained. The expression of the discrete equations is: ; In the formula, and Both represent node numbers. , ,and ; Represents a node The number of nodes in a local point cloud; Indicates in Time, node With nodes The relative permeability of oil or water between them; Indicates in Time, node With nodes The volume coefficient between oil and water; Indicates in Time, node With nodes The viscosity between oil and water; Indicates in Time, node With nodes Transfer rate between; Indicates in Time, node With nodes The pressure between oil or water; Indicates in Time, node Pressure of oil or water; Represents a node Source / collection; Indicates the volume controlled by the node; Indicates the time step; Indicates in Time, node Porosity; Indicates in Time, node Porosity; Indicates in Time, node The saturation level of oil or water; Indicates in Time, node The saturation level of oil or water.
[0022] S4, Beggs-Brill wellbore multiphase flow calculation: Using the wellhead pressure of S2 and the total flow rate of S3 as inputs, and based on the wellbore geometry and fluid properties, a segmented advancement-intra-segment iteration solution strategy is employed. Flow pattern determination, frictional pressure drop calculation, gravity term calculation, and acceleration term calculation are performed sequentially to obtain the calculated values of the pressure distribution along the wellbore and the layered bottomhole flowing pressure. .
[0023] S5. Coupled iterative convergence and relaxation update; calculation of predicted values of bottom-hole flowing pressure in stratified wells. Calculated values of bottom-hole flowing pressure in stratified wells The absolute value of the error is calculated and tested independently for injection wells and production wells; if the absolute value of the error is not less than the threshold... If the layered bottom hole pressure is not found, the layered relaxation factor is used to update the layered bottom hole pressure through a coupling update strategy to obtain the new predicted value of the layered bottom hole pressure, and the process returns to S3; otherwise, the coupling convergence is successful, and the next step continues. The updated expression is: ; In the formula, Indicates the first The predicted values of the bottom-hole flowing pressure of the stratified well are obtained after the next iteration update; Indicates the first The predicted value of the layered bottom flow pressure set in the next iteration; Represents the global relaxation factor. The current situation is relatively stable; for production wells, , For water injection wells, , ; The layer weight matrix represents the layer weight matrix that satisfies the following conditions: , Indicates the first Layer weight coefficients; This represents the layered bottom-hole flowing pressure value calculated using a multiphase flow model of the wellbore. Iterative convergence is determined by the dimensionless or absolute change in bottom-hole pressure, expressed as follows: ; In the formula, Indicates all layers Take the maximum value.
[0024] S6, Stage Advancement and Simulation Termination: Based on the converged bottom-hole flowing pressure of each layer, complete the EFVM reservoir seepage simulation calculation for the current time stage; determine whether the total simulation time has been reached. If not, update the reservoir state as the initial value for the next time stage and return to S2; if the total simulation time has been reached, output the simulation results; the simulation results include wellhead pressure, bottom-hole pressure of each layer, layer flow rate, total flow rate, pressure distribution along the wellbore, reservoir pressure field, and saturation field.
[0025] Example 2 An irregular boundary reservoir model was adopted, vertically divided into five layers: the first layer is 1000m to 1010m, the second layer is 1010m to 1020m, the third layer is 1020m to 1025m, the fourth layer is 1025m to 1030m, and the fifth layer is 1030m to 1035m. One injection well and two production wells, all vertical, were defined. One injection well was located at (0,60) and another at (460,160), with a bottomhole flowing pressure of 20MPa. One production well was located at (0,0), with a bottomhole flowing pressure of 10MPa. Table 1 shows the reservoir rock and fluid properties.
[0026] Table 1. Reservoir rock and fluid physical properties
[0027] This embodiment adopts Figure 2 The point cloud shown is used for calculation, simulating a production period of 150 days, with each stage lasting 5 days, for a total of 30 stages. The method of this invention performs a wellbore-reservoir coupling simulation once in each stage. Please refer to... Figure 3-9 , Figure 6 The total extraction volume is the total flow rate. In this embodiment, the second layer exhibits a higher flow rate than the first layer, indicating that the stratified injection and production behavior is no longer determined by the layer sequence, but by the combined effects of permeability, porosity, layer thickness, and near-well connectivity. The wellbore pressure of all wells gradually decreases along the depth direction, reflecting the dynamic coupling characteristics of wellbore flow and formation seepage.
[0028] Example 3 This invention was applied to an actual block in an oilfield. The reservoir is a medium-to-high porosity, high-permeability, and highly heterogeneous clastic rock reservoir. The reservoir porosity is generally concentrated, with an average porosity of 0.215 and an average permeability of 322.544 mD. The average single-layer thickness is 1.14 m, which is generally thin, but the multiple layers stacked together form an effective reservoir, giving the reservoir a certain degree of vertical continuity while also exhibiting significant differences in physical properties between layers. For reservoir distribution details, please refer to [link to relevant documentation]. Figure 10 Reservoir connectivity diagram and well distribution as follows Figure 10 As shown in (a) above, the three-dimensional permeability distribution map of the third reservoir is as follows: Figure 10 As shown in (b), there are 8 wells in the block, including 3 water injection wells and 5 production wells. The overall reservoir and well group distribution is shown in the three-dimensional diagram. Figure 10 (c) The thin-layered, high-permeability, and highly heterogeneous characteristics of this reservoir determine that fluids are prone to interlayer differential flow during reservoir development, which places higher demands on wellbore reservoir coupling calculations, layered production capacity allocation, and injection-production control.
[0029] This embodiment selects 8 vertical wells for wellbore reservoir coupling calculations. Please refer to [link to relevant documentation]. Figure 11-16This includes five production wells (P10, P20, P3, P15, P1) and three injection wells (P9, P23, P11). Decreasing and increasing pressure curves were selected for the injection and production wells, respectively. Figure 11 The pressure change curves at the wellheads of injection wells and production wells are shown separately. Figure 12 The simulation shows the oil saturation distribution of four oil reservoir layers during the final 150-day stage. Figure 13 The simulation shows the pressure distribution of four reservoir layers during the final 150-day stage. Figure 14 The curves showing the total fluid extraction volume of injection wells and production wells over time are presented. Figure 15 The flow rate variation curves of injection and production wells are presented. The flow rate varies significantly and non-monotonicly with each layer, with different layers alternating as the main injection and production contributing layers at different times, reflecting the influence of strong heterogeneity and complex connectivity in real reservoirs. This result indicates that under actual reservoir conditions, the flow rate of each layer is jointly controlled by multiple physical property parameters, and therefore does not follow a simple sequence gradient law. Figure 16 The stratified bottomhole flowing pressure variation curves of injection and production wells are shown. Compared with Example 2, it can be found that although the stratified flow rate and stratified bottomhole flowing pressure in this example no longer follow a simple sequence progression law, the injection and production wells maintain good consistency in terms of the total flow rate evolution trend, the wellbore pressure transmission process to the formation, and the temporal characteristics of the stratified response. This indicates that the established wellbore-reservoir coupling model can work stably under different well type combinations and strongly heterogeneous actual reservoir conditions, and reasonably characterizes the stratified seepage and wellbore pressure response characteristics under the combined action of multiple physical parameters, thus verifying the applicability and reliability of the model in actual reservoir engineering analysis.
[0030] Therefore, this invention adopts the above-mentioned meshless numerical simulation method based on EFVM for wellbore-reservoir coupling. By achieving strong coupling between the Beggs–Brill model and EFVM, the meshless feature is used to improve the geometric adaptability of complex reservoirs. The layered weighted relaxation factor is used to ensure the coupling convergence stability. It has both high computational accuracy and strong engineering applicability, providing reliable numerical support for injection-production optimization, production capacity evaluation and production system adjustment of various reservoirs.
[0031] 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. An EFVM-based meshless numerical simulation method for wellbore-reservoir coupling, characterized in that, Includes the following steps: S1. Construct the reservoir computational domain and discretize the point cloud; S2. Initialize simulation parameters; S3, EFVM reservoir seepage simulation; S4, Beggs-Brill wellbore multiphase flow calculation; S5. Coupled iterative convergence and relaxation update; S6, Phase Advancement and Simulation Termination.
2. The EFVM-based wellbore-reservoir coupled meshfree numerical simulation method of claim 1, wherein, In S1: According to the actual geological parameters of the oil reservoir, a meshless point cloud discrete calculation domain is generated, the point cloud includes a total of nodes, the point cloud includes internal nodes and boundary nodes of the oil reservoir, and a virtual node is added to the boundary nodes; Actual geological parameters of an oil reservoir include its geometry, layered permeability, porosity, and layer thickness.
3. The EFVM-based meshfree numerical simulation method of wellbore-reservoir coupling according to claim 2, characterized in that, The coordinates of the virtual nodes in S1 are calculated using the coordinates of the boundary nodes, the weighted average distance of the neighboring nodes, and the unit outward normal vector of the boundary. The expression is as follows: ； wherein denotes the coordinates of a virtual node ; denotes the coordinates of a boundary node ; denotes the unit exterior normal vector at a boundary node ; denotes the weighted average distance of a neighborhood node to a node .
4. The EFVM-based wellbore-reservoir coupled meshfree numerical simulation method of claim 3, wherein, In S2: Setting the total simulation time, time stage division, wellhead pressure time function, wellbore parameters, fluid parameters and convergence threshold, setting the predicted value of the first time layered bottom hole flowing pressure , the predicted value of the layered bottom hole flowing pressure of each subsequent time stage is the layered bottom hole flowing pressure after convergence of the previous time stage; The wellhead pressure-time function includes the production wellhead pressure-time function and the injection wellhead pressure-time function. The production wellhead pressure decreases with time, while the injection wellhead pressure increases with time. Wellbore parameters include wellbore diameter, wellbore length, and wellbore roughness; Fluid parameters include fluid viscosity, fluid density, and fluid phase characteristics.
5. The EFVM-based meshfree numerical simulation method of wellbore-reservoir coupling according to claim 4, characterized in that, In S3: Based on the S1-based point cloud, a local support domain and an equivalent control volume are generated for each node, the two-phase control equation of the reservoir is discretized by the EFVM method, and the predicted value of the bottom hole pressure of the layered well is obtained As a boundary condition, the oil reservoir pressure field and saturation field are solved by Newton iteration, and the total flow rate of each well is output And the layered flow rate .
6. The EFVM-based wellbore-reservoir coupled meshfree numerical simulation method of claim 5, wherein, The specific process of discretizing the reservoir two-phase governing equations using the EFVM method in S3 is as follows: A nodal generalized difference operator is constructed using Taylor expansion and weighted least squares method. Combined with the nodal control volume and two-point flux scheme, a discrete equation satisfying local mass conservation is obtained. The expression of the discrete equation is: ； In the formula, and Both represent node numbers. , ,and ; Represents a node The number of nodes in a local point cloud; Indicates in Time, node With nodes The relative permeability of oil or water between them; Indicates in Time, node With nodes The volume coefficient between oil and water; Indicates in Time, node With nodes The viscosity between oil and water; Indicates in Time, node With nodes Transfer rate between; Indicates in Time, node With nodes The pressure between oil or water; Indicates in Time, node Pressure of oil or water; Represents a node Source / collection; Indicates the volume controlled by the node; Indicates the time step; Indicates in Time, node Porosity; Indicates in Time, node Porosity; Indicates in Time, node The saturation level of oil or water; Indicates in Time, node The saturation level of oil or water.
7. The meshless numerical simulation method for wellbore-reservoir coupling based on EFVM according to claim 6, characterized in that, In S4: Using the wellhead pressure of S2 and the total flow rate of S3 as inputs, and based on the wellbore geometry and fluid properties, a segmented advancement and intra-segment iteration solution strategy is employed. This involves sequentially performing flow pattern determination, frictional pressure drop calculation, gravity term calculation, and acceleration term calculation to obtain the calculated values of the pressure distribution along the wellbore and the layered bottom-hole flowing pressure. .
8. The meshless numerical simulation method for wellbore-reservoir coupling based on EFVM according to claim 7, characterized in that, In S5: Calculate the predicted value of bottom-hole flowing pressure in a stratified well. Calculated values of bottom-hole flowing pressure in stratified wells The absolute value of the error is calculated and tested independently for injection wells and production wells; if the absolute value of the error is not less than the threshold... Then, the layered bottom hole flowing pressure is updated through the coupling update strategy of the layered relaxation factor to obtain the new predicted value of the layered bottom hole flowing pressure, and then returned to S3; Conversely, if the coupling converges successfully, proceed to the next step; the updated expression is: ； In the formula, Indicates the first The predicted values of the bottom-hole flowing pressure of the stratified well are obtained after the next iteration update; Indicates the first The predicted value of the layered bottom flow pressure set in the next iteration; Represents the global relaxation factor; The layer weight matrix represents the layer weight matrix that satisfies the following conditions: , Indicates the first Layer weight coefficients; This represents the layered bottom flow pressure value calculated by the wellbore multiphase flow model.
9. The meshless numerical simulation method for wellbore-reservoir coupling based on EFVM according to claim 8, characterized in that, In S5, iterative convergence is determined by the dimensionless or absolute change in bottom hole pressure, expressed as follows: ； In the formula, Indicates all layers Take the maximum value.
10. The meshless numerical simulation method for wellbore-reservoir coupling based on EFVM according to claim 9, characterized in that, In S6: Based on the converged bottom-hole flowing pressure of each layer, the EFVM reservoir seepage simulation calculation for the current time stage is completed; it is determined whether the total simulation time has been reached. If not, the reservoir state is updated as the initial value for the next time stage, and the process returns to S2; if the total simulation time has been reached, the simulation results are output. The simulation results include wellhead pressure, bottom-hole pressure of each layer, layer flow rate, total flow rate, pressure distribution along the wellbore, reservoir pressure field, and saturation field.