Novel deep coal bed gas fracture network seepage simulation method based on connection element system

Abstract

The invention discloses a novel deep coal bed gas fracture network seepage simulation method based on a connection element system, which comprises the following steps: constructing a fracture network connection relationship by constructing a point set and a connection screening mode; calculating seepage characteristic parameters of any connection unit based on the constructed fracture network connection relation; and solving a local dual-medium two-phase flow seepage control equation on the basis of the constructed fracture network connection relationship and the obtained seepage characteristic parameters to obtain the gas production rates of different perforation positions. By the adoption of the method, the meshless connection network is constructed according to the connectivity and diversion characteristics in the multi-fracture horizontal well, multiphase flow in an unconventional oil and gas reservoir is efficiently simulated, the discrete scale and calculation time are remarkably reduced while the calculation precision is kept or improved, and the calculation efficiency is improved. And an engineered technical path is provided for hydraulic fracturing design and production optimization of unconventional oil and gas reservoirs.

Description

Technical Field

[0001] This invention relates to the field of unconventional oil and gas field development engineering technology, and in particular to a novel method for simulating seepage in deep coalbed methane fracture networks based on a connecting element system. Background Technology

[0002] Unconventional oil and gas resources, with their vast resource potential found in various strata such as shale, tight sandstone, and tight carbonate rocks, are gradually becoming an important alternative energy option. However, compared with conventional oil and gas, unconventional oil and gas reservoirs generally exhibit low porosity, ultra-low permeability, and significant formation heterogeneity. They are also mostly located in deep strata, accompanied by extreme geological conditions such as high temperature, high pressure, and complex stress fields, thus posing more severe engineering challenges to development and exploitation.

[0003] In the industrial development of unconventional oil and gas, hydraulic fracturing technology is a key means of enhancing the recovery rate of unconventional oil and gas. The complex fracture network formed by it directly determines the fluid migration path and production capacity distribution. At present, the mainstream flow characterization methods mainly include three types: continuous medium model, discrete model, and embedded discrete model. However, the above methods have significant shortcomings: the continuous medium model requires a large number of grids to describe the fracture morphology in detail, which leads to a sharp increase in computation time; the discrete fracture model (DFM) faces problems such as complex grid subdivision and poor grid quality, resulting in insufficient adaptability; although the embedded discrete fracture model (EDFM) does not require local grid refinement, it has errors in the simulation of low conductivity fractures and pressure distribution assumptions, and its accuracy is limited in scenarios such as fracture confluence and matrix coupling.

[0004] In recent years, meshless methods have been increasingly used in fluid mechanics. These methods use point clouds instead of grids to discretize the computational domain, allowing for more flexible characterization of fractures and complex boundary conditions. Common implementations include weighted least squares, elementless Galerkin, and generalized finite difference. While they have advantages in characterizing complex boundaries, existing techniques still have shortcomings in representing inter-well connectivity and controlling computational costs. Although new methods such as the meshless connective element method (CEM) have provided new ideas for reservoir simulation, traditional numerical methods face difficulties in adapting to irregular fractures, controlling computational costs, and improving history fitting efficiency. Therefore, there is an urgent need to provide new methods for simulating seepage in deep coalbed methane fracture networks based on connective element systems. Summary of the Invention

[0005] The purpose of this invention is to provide a novel method for simulating seepage in deep coalbed methane fracture networks based on a connective element system. This method can not only accurately reproduce the pressure distribution and fracture connectivity characteristics, but also significantly improve computational efficiency compared to traditional grid methods. It also possesses good numerical stability and scalability, providing a theoretical basis and decision support for the design of fracturing parameters in unconventional oil and gas wells.

[0006] To achieve the above objectives, this invention provides a novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system, comprising the following steps: S1. Construct the mesh connection relationship by building point sets and using connection filtering methods; S2. Based on the fracture-mesh connection relationship constructed in S1, calculate the seepage characteristic parameters of any connection unit. The seepage characteristic parameters include: the conductivity of the connection unit between matrix nodes, the conductivity between two fracture nodes in the fracture layer, and the conductivity between the matrix and the fracture. S3. Based on the seam-mesh connection relationship of S1 and the seepage characteristic parameters obtained in S2, solve the local dual-medium two-phase flow seepage control equation to obtain the gas production rate at different perforation locations.

[0007] Preferably, S1 is as follows: S11. Construct a fracture-matrix point set to provide basic data for the subsequent construction of the fracture network connection relationship; S12. Select reasonable connection units based on the influence domain radius and the maximum control angle; S13, based on the crack-matrix point set of S11 and the reasonable connection of S12, forms a crack network connection relationship.

[0008] Preferably, the construction of the crack-matrix point set in S11 is specifically as follows: (1) Crack points: Extracted based on the binarization result of the original crack structure image. A value of 1 represents the crack location and is directly converted into crack point coordinates. (2) Matrix points: The Poisson disk sampling method is used to deploy matrix points to obtain a matrix point set with uniform spatial distribution and appropriate spacing between them; (3) When the Euclidean distance between any matrix point and any crack point is less than the set minimum distance threshold, the matrix point is discarded.

[0009] Preferably, S12 is as follows: (1) Establish initial connection units based on the radius of the influence domain: Using any node as the center, set a fixed influence domain radius; traverse all other nodes, and if the Euclidean distance between nodes is less than or equal to the set fixed influence domain radius, then construct an initial connection unit between the two nodes; otherwise, do not connect them. (2) Based on the maximum control angle, reasonable connection units are selected and retained. The selection rules are as follows: Two nodes of any initial connection unit and a third node form a triangle. Calculate the cosine of the included angle between the two sides of the initial connection unit in the triangle. If there is an included angle cosine that is less than the cosine of the maximum control angle, i.e., the included angle is greater than the maximum control angle, then the initial connection unit is discarded. If all included angle cosine values ​​in all the triangles are greater than or equal to the cosine of the maximum control angle, then the initial connection unit is retained as a valid connection unit.

[0010] Preferably, S2 is as follows: S21. Based on the principle of conservation of matter, define the node control domain; S22. Based on the two-phase flow control equation of the local dual-medium matrix system, the conductivity of the inter-matrix node connection unit is obtained. S23. Based on the two-phase flow seepage control equation of a local dual-medium fracture system, the conductivity between two fracture nodes in the fracture layer is defined. S24 characterizes the conductivity between the two-phase flow matrix and the crack.

[0011] Preferably, the conductivity of the inter-matrix node connection unit in S22 is as follows: ; in, It is the conductivity of the inter-node connection unit in the matrix; It is a matrix node , Intermediate and average penetration rates; It is a matrix node Control volume; It is the number of matrix nodes; It is a connection unit At the node The volume distribution coefficient in the controlled volume.

[0012] Preferably, the conductivity between two fracture nodes in the fracture layer in S23 is as follows: ; in, It is the conductivity between two fracture nodes in the fracture layer; It is a node and The mobility values ​​between; It is a node in the crack system. , Harmonized average permeability; The thickness of the gas reservoir; For nodes With nodes The average crack aperture; For connection unit The length of the crack.

[0013] Preferably, the conductivity between the matrix and the fracture in S24 is as follows: ; in, It is the conductivity between the matrix and the crack; It is the height of the crack; It is the flow coefficient of the corresponding phase in a two-phase flow; It is a node The harmonized average permeability of the matrix and fracture flow within the control domain; Taken as all nodes The sum of half of the crack connection elements at the endpoints is specifically expressed as: , Represented by node The set of other end nodes of the crack connection unit at the endpoint; Based on nodes and nodes Crack connection unit at the endpoint Length; It is the equivalent normal distance between the crack and the matrix flow, specifically expressed as: .

[0014] Preferably, S3 is as follows: S31. Combining the node control domain, the two-phase flow seepage control equations of the local dual-medium matrix system and the local dual-medium fracture system in S2 are discretized; based on the fracture network connection unit constructed in S1 and the conductivity obtained in S2, the discretized local dual-medium two-phase flow seepage control equations are solved to obtain the fluid flow rate of the node connection unit. S32. Based on the obtained conductivity, fluid flow rate and reasonable connection volume of the connecting unit, calculate the gas production rate at different perforation positions.

[0015] Preferably, the discretized local dual-medium two-phase flow seepage control equation in S31 is as follows: ; in, , They are nodes The pressure; It is to depict nodes Within the control domain, the amount of mass exchanged between matrix nodes and fracture nodes; , They are t , The porosity at time t, where porosity is a function of pressure, is specifically expressed as: ; ; , They are time points , t Time node Phase saturation; This is the penetration rate, taken as the harmonic average, specifically expressed as: , yes Time Node Absolute penetration rate at the location; yes Time Node j Absolute penetration rate at the location; It is the viscosity, taken as the arithmetic mean, specifically expressed as: , yes Phase at time node Viscosity at that point; yes Phase at time node j Viscosity at that point; It is the relative permeability of the gas phase and the aqueous phase, specifically: , It is the equivalent relative permeability of the corresponding phase on the connecting unit; It is a moment t arrive Nodes within the time period Water saturation; It is a moment t arrive During the time period, the connecting unit , j The pressure of the corresponding phases; It is a moment t arrive During the time period, the connecting unit The pressure corresponding to the upper phase.

[0016] Therefore, the novel method for simulating seepage in deep coalbed methane fracture networks based on the aforementioned interconnected element system has the following beneficial effects: (1) This method generates a uniform matrix point cloud by sampling with a Poisson disk, which effectively eliminates the common element distortion problem in traditional mesh generation; then, the complex crack morphology is directly discretized into a point cloud, preserving the trunk and branch features completely; finally, a point conduction capacity model is constructed based on the crack connection element method (FCEM) theory, which realizes the accurate description of the multiphase flow relationship between cracks and matrix.

[0017] (2) This method can accurately reproduce the fracture connectivity structure in both the conceptual example and the actual case of a dual-horizontal well in a coalbed methane reservoir. It also demonstrates significant advantages in computational efficiency while maintaining high accuracy. In summary, the Fracture Connectivity Element Method (FCEM) combines accuracy and efficiency in the simulation and optimization of unconventional oil and gas hydraulic fracturing. It can effectively cope with extreme geological conditions in deep formations and has good prospects for engineering application.

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

[0019] Figure 1 This is a diagram showing the relationship between the influence domain and the angle limitation in Embodiment 1 of the present invention; Figure 2 This is a schematic diagram of the morphology and permeability distribution characteristics of a three-cluster fracture network in a horizontal well of a coalbed methane reservoir in Embodiments 2 and 3 of the present invention, and the reservoir connection conductivity and connection volume calculated by combining seepage characteristic parameters. In this diagram, (a) is the permeability distribution; (b) is the generated connection; (c) is the connection conductivity; and (d) is the connection volume. Figure 3 This is a graph showing the effect of different coarsening levels on simulation accuracy and computation time in Embodiment 2 of the present invention. Among them, (a) is a connection model diagram with 26 corresponding lines; (b) is a connection model diagram with 24 corresponding lines; (c) is a connection model diagram with 16 corresponding lines; and (d) is a connection model diagram with 8 corresponding lines. Figure 4 This is a graph showing the impact of the number of retained paths on the accuracy and computational efficiency of gas production prediction in Embodiment 2 of the present invention. (a) is a comparison of daily gas production of a single well; (b) is a comparison of the impact of retained paths. Figure 5 This is the effect of different point cloud densities on simulation accuracy and computation time in Embodiment 2 of the present invention, wherein (a) is the point cloud distribution map corresponding to a minimum point cloud distance of 30; (b) is the point cloud distribution map corresponding to a minimum point cloud distance of 25; (c) is the point cloud distribution map corresponding to a minimum point cloud distance of 20; and (d) is the point cloud distribution map corresponding to a minimum point cloud distance of 15. Figure 6 This is a graph showing the impact of the minimum point cloud distance on the accuracy and computational efficiency of gas production prediction in Embodiment 2 of the present invention. (a) shows the daily gas production of a single well; (b) shows a comparison of the impact of the minimum point cloud distance. Figure 7This is a graph showing the effect of different maximum control angles on simulation accuracy and computation time in Embodiment 2 of the present invention. Among them, (a) is the connection distribution diagram corresponding to a maximum control angle of 115°; (b) is the connection distribution diagram corresponding to a maximum control angle of 125°; (c) is the connection distribution diagram corresponding to a maximum control angle of 135°; and (d) is the connection distribution diagram corresponding to a maximum control angle of 145°. Figure 8 This is a graph showing the impact of the maximum control angle on the accuracy and efficiency of gas production prediction in Embodiment 2 of the present invention. In this graph, (a) is the daily gas production of a single well; and (b) is a comparison of the impact of the maximum control angle. Figure 9 This is a comparison diagram of the conceptual calculation examples of daily gas production per well simulated by the method of the present invention and the tNavigator method in Embodiment 3 of the present invention; Figure 10 This is a fracture morphology and connection model of an actual well case in Embodiment 3 of the present invention, wherein (a) is the permeability distribution of the actual gas reservoir model of a coalbed methane double-horizontal well under the same working conditions; and (b) is the modeling result of modeling and simulating the gas reservoir using the method of the present invention. Figure 11 This is a single-well daily gas production curve diagram in Embodiment 3 of the present invention, wherein (a) is a comparison diagram of the results obtained by the first simulation of well A using this method and the actual production data; (b) is a comparison diagram of the results obtained by the first simulation of well B using this method and the actual production data. Figure 12 This is a curve of the historical fitting result of daily gas production of a single well in Embodiment 3 of the present invention. Among them, (a) is a comparison chart of the result obtained by historical fitting of well A using this method and the actual production data; (b) is a comparison chart of the result obtained by historical fitting of well B using this method and the actual production data. Detailed Implementation

[0020] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0021] This invention presents a novel method for simulating seepage in deep coalbed methane fracture networks based on a interconnected metasystem, comprising the following steps: S1. Construct the mesh connection relationship by building point sets and using connection filtering methods; S2. Based on the fracture-mesh connection relationship constructed in S1, calculate the seepage characteristic parameters of any connection unit. The seepage characteristic parameters include: the conductivity of the connection unit between matrix nodes, the conductivity between two fracture nodes in the fracture layer, and the conductivity between the matrix and the fracture. S3. Based on the seam-mesh connection relationship of S1 and the seepage characteristic parameters obtained in S2, solve the local dual-medium two-phase flow seepage control equation to obtain the gas production rate at different perforation locations.

[0022] Example 1 A novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity-based system includes the following steps: S1. Construct the connection relationship of the mesh by building point sets and using connection filtering methods.

[0023] S11. Construct a crack-matrix point set to provide basic data for subsequent connection relationship construction.

[0024] (1) Crack points: Extracted based on the binarization result of the original crack structure image (i.e., Excel data file with 0-1 distribution). The 1 value represents the crack position and is directly converted into crack point coordinates. This method can ensure that the crack morphology is consistent with the real image.

[0025] (2) Matrix points: The Poisson disk sampling method is used to deploy matrix points to obtain a spatially uniform point set with appropriate spacing between them; the specific expression is shown below: ; in, It is a matrix point The spatial coordinate vector; It is a matrix point The spatial coordinate vector; It is the point set obtained by sampling from the Poisson disk; It is the radius of the Poisson disk; Therefore The process of sampling Poisson disk for control parameters.

[0026] (3) In order to avoid the matrix point and the crack point overlapping or being too close to interfere with the crack connectivity analysis, a minimum distance threshold is set. When the Euclidean distance between any matrix point and any crack point is less than the threshold, the matrix point will be removed.

[0027] By employing the above-mentioned point placement strategy, a crack-matrix point set with a clear structure and reasonable physical properties is obtained, providing support for the construction of connection units.

[0028] S12. Select reasonable connection units based on the radius of influence and the maximum control angle.

[0029] (1) Establish initial connection units based on the radius of the influence domain.

[0030] like Figure 1 As shown, for each node, a fixed influence radius is set with that node as the center. ; Traverse all other nodes, if the Euclidean distance between each node and the central node is less than or equal to If the two nodes are connected, an initial connection unit is established; otherwise, no connection is made. Based on the influence domain, the point... , The connection between them is: ; in, It is a node With nodes The connection criterion; It is a node Spatial coordinate components; It is a node Spatial coordinate components.

[0031] To avoid redundant connections or unreasonable short connections, all initially established connection units are screened.

[0032] (2) Select and retain reasonable connection units based on the maximum control angle.

[0033] Introducing the maximum control angle ( Using this as a criterion, the initial connection units are screened.

[0034] For any initial connection element, the two nodes it connects can form a triangle with a third node. Calculate the cosine of the angle between the two sides of the initial connection element and the other sides of this triangle. If there exists any triangle such that If the angle is greater than the maximum control angle, the connection is deemed unreasonable and discarded; if the cosine of the included angle in all related triangles is greater than or equal to... If so, the connection unit is considered to be reasonably retained.

[0035] S2. Based on the fracture-mesh connection relationship constructed in S1, calculate the seepage characteristic parameters of any connection unit. The seepage characteristic parameters include: the conductivity of the connection unit between matrix nodes, the conductivity between two fracture nodes in the fracture layer, and the conductivity between the matrix and the fracture.

[0036] Since the nodes of the connecting metasystem do not have a substantial control domain, in order to handle physical problems with source and sink terms, S21, based on the principle of matter conservation, defines the node control domain as follows: ; in, It is the entire gas reservoir control area; It is a node The control domain; It is the total volume of the gas reservoir; The node controls the volume; It is a node The control domain; It is the set of indices of all nodes / control volumes involved in the subdivision; the control domains between all nodes do not intersect, and the sum of all control domains is the entire gas reservoir region.

[0037] The following section uses the two-phase flow control equation as an example to derive the calculation method for the seepage characteristic parameters of the connecting unit.

[0038] S22. Based on the two-phase flow control equation of the local dual-medium matrix system, the conductivity of the inter-matrix node connection unit is obtained.

[0039] First, the governing equations for two-phase flow in a local dual-medium matrix system: ; in, and These are the absolute permeability and relative permeability of the matrix system, respectively. ; Indicates the viscosity of gases and water. ; It is the formation volume coefficient, which is dimensionless; It is the Hamiltonian gradient operator; It is the pressure of the matrix layer. ; This refers to the volumetric flow rate (source and sink parameters) under standard ground conditions. ; It is the Dirac function, which is dimensionless; It is the saturation of the matrix layer, which is dimensionless; It is time. ; It represents the amount of mass exchanged between the matrix and the crack, indicating the flow of gas or water between the two systems, with units consistent with source and sink terms. It is a gas phase; It is an aqueous phase.

[0040] Secondly, the pressure diffusion term in the above formula is applied in the nodal control domain. Integral inner, we get: ; in, It is a node Controlling volume, ; It is a control domain Average penetration rate ; It is the second-order spatial variation of the matrix system pressure field within the control volume; within the nodal control domain, this region is considered locally homogeneous, and the permeability is... .

[0041] Then, adopt Formulas and weighted least squares approximation methods, utilizing nodes The partial derivative estimate of the stress function value within the influence domain is obtained for that node, as shown below: ; in, It is the central node of the influence domain The Laplace value of the pressure function at that point; It is the central node Pressure value at the location; It is the first in the influence domain Pressure values ​​at adjacent nodes; yes The coefficients of the operator estimate , Indicates the sequence number of other nodes in this influence domain; Indicates the number of nodes within the affected domain (excluding the central node); The values ​​correspond to the partial derivatives of the pressure function.

[0042] Next, substituting the partial derivative estimate into the integral equation, we transform it into the form of pressure difference between nodes multiplied by conductivity: ; Among them, for a certain influence domain center node The area of ​​influence includes other 1 node Represented by nodes Connection unit calculated for the central node Conductivity, considering nodes The central influence domain (including nodes) This area of ​​influence includes other... From these nodes, we can obtain: ; in, Represented by nodes Connection unit calculated for the central node The conductivity. For typical heterogeneous gas reservoirs, the connecting unit conductivity Therefore, it is necessary to take the harmonic average of the penetration rates of the two, i.e.: ;in, It is a node With nodes The conductivity of the connection units between them; Based on nodes Centered on nodes Connection unit obtained from physical property calculation on one side unilateral conductivity; Similarly; It is a node , The reverse equivalent connection conductivity between them.

[0043] According to the principle of conservation of matter, in the connecting unit The above satisfies ,in, It is a connection unit At the node The volume allocation coefficient in the controlled volume; It is a connection unit At the node The volume distribution coefficient in the control volume, i.e., the harmonic average of conductivity is equivalent to the harmonic average of permeability.

[0044] Based on the principle of material balance connecting elements and the fact that the sum of the control volumes of all nodes equals the gas reservoir volume. d Based on the principle, the control volume of each node is derived, thereby allowing the definition of the connection elements between matrix nodes. Its conductivity is: .

[0045] S23. Based on the two-phase flow control equation of the local dual-medium fracture system, the conductivity between two fracture nodes in the fracture layer is defined.

[0046] The governing equations for two-phase flow seepage in a local dual-medium fracture system are as follows: ; in, and These are the absolute and relative permeability of the crack, respectively. ; Indicates the viscosity of gases and water. ; It is the pressure of the matrix layer. ; It is the saturation of the matrix layer, which is dimensionless; It is the volumetric flow rate (source and sink) under standard ground conditions. .

[0047] The conductivity between two fracture nodes in a fracture layer is defined as: ; in, It is the conductivity between two fracture nodes in the fracture layer; It is a node in the crack system. , Harmonized average permeability between ; The thickness of the gas reservoir; For nodes With nodes The average crack aperture; For connection unit The length of the crack; It is a node and The flow rate values ​​between.

[0048] S24 characterizes the conductivity between the two-phase flow matrix and the crack.

[0049] First, define the interface area between the matrix and the crack.

[0050] Since the connecting element system does not describe the specific shape of the control volume of the nodes, the total length of the cracks within the node control domain is used. With crack height The product is used to characterize the interfacial area between the crack and the matrix. .

[0051] For a certain central node , Taken as all nodes The sum of half of the crack connection elements at the endpoints is shown below: ; in, Represented by node The set of other end nodes of the crack connection unit at the endpoint; Based on nodes and nodes Crack connection unit at the endpoint Length, .

[0052] Secondly, the mass exchange rate between matrix nodes and fracture nodes is defined as follows: ; in, It is a node The harmonized average permeability of the matrix and fracture flow within the control domain, ; It is the flow coefficient of the corresponding phase in a two-phase flow; It is a node The pressure of the matrix layer within the control domain; It is a node The pressure of the fracture layer within the control domain; It is a node The exchange conductivity coefficient of the matrix fractures; It is the equivalent normal distance between the crack and the matrix flow. Considering that the fracture extends through the entire thickness of the gas reservoir, i.e., the fracture height is equal to the gas reservoir thickness. The same, therefore, .

[0053] Finally, the conductivity between the matrix and the crack is defined as follows: .

[0054] S3. Based on the seam-mesh connection relationship of S1 and the seepage characteristic parameters obtained in S2, solve the local dual-medium two-phase flow seepage control equation to obtain the gas production rate at different perforation locations.

[0055] S31. Combining the node control domain, the two-phase flow control equations of the local dual-medium matrix system and the local dual-medium fracture system in S2 are discretized: Based on the fracture network connection unit constructed in S1 and the conductivity obtained in S2, the discretized local dual-medium two-phase flow control equations are solved to obtain the fluid flow rate of the node connection unit. The discretized local dual-medium two-phase flow seepage control equations are as follows: ; in, Connecting units between matrix nodes conductivity, ; , They are nodes Pressure ; It is to depict nodes Within the control domain, the amount of mass exchanged between matrix nodes and fracture nodes; , They are t , The porosity at time t, where porosity is a function of pressure, is specifically expressed as: ; ; , They are time points , t Time node Phase saturation; This is the penetration rate, taken as the harmonic average, specifically expressed as: , It is the node at that moment. Absolute penetration rate at the location; It is the node at that moment. j Absolute penetration rate at the location; It is the viscosity, taken as the arithmetic mean, specifically expressed as: , Is the current phase state at node Viscosity at that point; Is the current phase state at node j Viscosity at that point; It is the relative permeability of the gas phase and the aqueous phase, specifically expressed as: , It is the equivalent relative permeability of the corresponding phase on the connecting unit; It is a node within this period of time. Water saturation; During this period, the connecting unit , j The pressure of the corresponding phases; During this period, the connecting unit The pressure corresponding to the upper phase.

[0056] S32. Based on the obtained conductivity, fluid flow rate and reasonable connection volume of the connecting unit, calculate the gas production rate at different perforation positions.

[0057] Example 2 To avoid the exponential expansion of connectivity due to an excessive number of nodes during the initial discretization of complex crack morphology, this study introduces a simpler descriptive structure while maintaining the basic geometric characteristics of the crack.

[0058] Specifically, this method simplifies complex fracture networks by employing morphological skeletonization and combining it with Dijkstra's shortest path algorithm. This method significantly reduces the number of nodes while preserving the fracture backbone topology, thus improving the efficiency of subsequent flow simulations while maintaining computational accuracy. Furthermore, the discretization of the matrix reservoir also requires careful consideration: node density and the number of connections directly affect computational results and speed. Therefore, ablation experiments are conducted focusing on three key variables—fracture morphology coarsening, matrix point cloud density, and maximum control angle—to determine the optimal modeling criteria.

[0059] Crack coarsening effect: Based on Dijkstra's algorithm, all possible paths to the crack endpoint can be searched sequentially from the perforation point. The coarsening criterion of this method is to prioritize retaining several of the longest paths among all paths. For example... Figure 2 The actual fracture morphology shown in Figure (a) was used to select 26, 24, 16, and 8 paths respectively as retained objects, and a comparative analysis was conducted under the same matrix point cloud density and maximum connection radius. The evaluation indexes were the RMSE of the gas production curve simulated by this method and the gas production curve simulated by tNavigator, as well as the simulation speed of the two methods. The results are shown in Table 1.

[0060] Table 1. Effects of different coarsening levels on simulation accuracy and computation time.

[0061] like Figure 4 As shown, 26 paths resulted in the lowest error (9.14%), but the computation time was relatively long (12.9s). When the number of paths was reduced to 8, the computation speed was the fastest (5.6s), but this was accompanied by the highest error (21.17%), making it difficult to guarantee simulation accuracy. In contrast, retaining 24 paths achieved a better balance between error (9.97%) and computation time (6.3s), significantly shortening the computation time while maintaining high accuracy. Therefore, considering both accuracy and efficiency, 24 paths were determined to be the optimal number of paths to retain in this method.

[0062] Influence of matrix point cloud density: The parameter affecting the density of the matrix point cloud is the minimum distance between two points in the generated point cloud. A point in the point cloud is retained if its Euclidean distance to all other matrix points is greater than the given minimum distance. Similarly, Figure 2 The actual crack morphology shown in (a) was used as the minimum distance limit between point clouds, with values ​​of 30, 25, 20, and 15. A comparative analysis was conducted under the same crack morphology and maximum connection radius. The evaluation index was the RMSE of the gas production curve simulated by tNavigator and the calculation speed. The results are shown in Table 2.

[0063] Table 2. The impact of different point cloud densities on simulation accuracy and computation time.

[0064] like Figure 6 As shown, when the minimum distance is 30, although the calculation speed is the fastest (4.6s), the error is relatively large (18.69%). When the distance decreases to 15, the error is the lowest (8.45%), but the calculation time is extremely long (42.7s), which is difficult to meet the requirements of efficient simulation. In contrast, the calculation error is only 9.97% and the calculation time is 6.3s when the distance is 25, which is significantly better than the calculation efficiency of higher density while ensuring high accuracy. Therefore, considering both accuracy and efficiency, 25 is determined to be the optimal minimum distance for matrix point clouds.

[0065] The effect of the maximum control angle: In this method, the maximum control angle determines whether connectivity is established between nodes, thus directly affecting the number of connections generated by the model. To study the impact of this parameter on the simulation results, 115°, 125°, 135°, and 145° were selected as the maximum control angles, and the fracture morphology was kept unchanged for comparative analysis. The evaluation indicators are the RMSE of the gas production curve simulated by tNavigator and the computation speed, and the results are shown in Table 3.

[0066] Table 3. Effects of different maximum control angles on simulation accuracy and computation time.

[0067] like Figure 8 As shown, the calculation speed is fastest (4.9s) at an angle of 115°, but the error is relatively high (11.25%). When the angle increases to 145°, although the error is lowest (9.12%), it requires a longer calculation time (9.5s), and the efficiency drops significantly. In contrast, the calculation error at 135° is only 9.97%, and the calculation time is 6.3s, achieving an optimal balance between accuracy and efficiency. Therefore, considering both simulation accuracy and calculation speed, 135° is determined to be the optimal maximum control angle for this method.

[0068] Therefore, considering both computational accuracy and efficiency, the optimal parameter combination is selected as 24 retained paths, 25 minimum point cloud distances, and 135° maximum control angle. This combination can effectively improve the overall performance of flow simulation while maintaining the rationality of the crack network geometry and topology.

[0069] Example 3 like Figure 2 The three-cluster fracture network morphology and permeability distribution characteristics of a fractured horizontal well section in coalbed methane reservoirs shown in (a) and (b) are used to calculate the reservoir connectivity conductivity and connectivity volume based on seepage characteristic parameters. Figure 2 As shown in (c) and (d) in the table. Then, a flow simulation study of the seam mesh was conducted, and the basic parameters of the model are shown in Table 4.

[0070] Table 4 Parameters of the Actual Calculation Example

[0071] Assuming the fracture segment is simulated under a constant water production regime for 900 days, with a time step of 30 days, the gas production rate is calculated using the tNavigator numerical simulation software. Then, based on the obtained conductivity and connectivity volume, the gas production rate under the same water production regime is calculated using the fracture connective element method (FCEM) proposed in this study, and compared with tNavigator. Figure 9 As shown in the figure. The results show that, under the same conditions, the simulation time of this method is 6.3 seconds, while that of tNavigator is 369 seconds. The weighted relative error (WRE) of the gas production curves of the two methods is 9.97%, indicating that FCEM significantly improves computational efficiency while ensuring computational accuracy. Furthermore, it can accurately characterize fracture morphology using a node-based approach, enabling relatively fast and efficient simulation of fracture network flow.

[0072] To verify the stability of the model, based on Figure 10 Numerical tests were conducted on the actual coalbed methane reservoir model with dual horizontal wells shown in (a) under the same operating conditions.

[0073] Based on the model parameters in Table 5, this study uses the proposed fractured cohesive element method (FCEM) to model and simulate this gas reservoir. The modeling results are as follows: Figure 10 As shown in (b), using constant water production, the daily gas and water production curves of simulated wells A and B are as follows. Figure 6 As shown. Using one day as a time step, the simulation lasted 579 days, taking 84 seconds to complete, and yielded the simulated dynamic curves. The initial simulation showed that the daily gas production curves of wells A and B matched the actual gas production curves with 72.32% and 71.67% respectively. Figure 11 As shown.

[0074] Table 5 Parameters of the Actual Calculation Example

[0075] like Figure 11 As shown, for a real-world example with densely packed fractures, the Fracture Connector Element Method (FCEM) deviates significantly from the actual results. Therefore, this paper introduces the Synchronous Perturbation Stochastic Approximation (SPSA) algorithm for automatic history fitting. The objective function for history fitting is constructed using historical production dynamics data such as daily gas production and bottom hole flowing pressure, with the conductivity and control volume of each connector element as the independent variables. The SPSA algorithm parameters are set as follows: 100 iterations, 0.5 iteration step size, 5 perturbations, and 0.1 perturbation step size. The fitted results are shown below. Figure 12 As shown.

[0076] Based on the connectivity parameters calculated from the physical property parameters, the daily gas production curves of single wells obtained by FCEM through simulation and historical data fitting show a high degree of consistency with the measured trends, with fitting rates of 87.22% and 85.13% for well A and well B, respectively. Compared with traditional numerical simulations (such as ECLIPSE), which require discretizing approximately 1,360,000 grid cells for a gas reservoir of 1600m × 850m, the FCEM in this study describes the system in a meshless manner using only 502 connectivity element nodes and 1,748 connectivity paths. While maintaining simulation accuracy, the discretization scale is reduced by nearly 1,800 times, and the simulation time is reduced from several hours to several minutes, achieving significant computational acceleration and scale advantages.

[0077] Therefore, the novel method for simulating seepage in deep coalbed methane fracture networks based on the aforementioned interconnected element system has the following beneficial effects: (1) This method generates a uniform matrix point cloud by sampling with a Poisson disk, which effectively eliminates the common element distortion problem in traditional mesh generation; then, the complex crack morphology is directly discretized into a point cloud, preserving the trunk and branch features completely; finally, a point conductivity model is constructed based on FCEM theory, realizing an accurate description of the multiphase flow relationship between cracks and matrix.

[0078] (2) This method can accurately reproduce the fracture connectivity structure in both the conceptual example and the actual case of a dual-horizontal well in a coalbed methane reservoir. It also demonstrates significant advantages in computational efficiency while maintaining high accuracy. In summary, the FCEM method combines accuracy and efficiency in the simulation and optimization of unconventional oil and gas hydraulic fracturing. It can effectively cope with extreme geological conditions in deep formations and has good prospects for engineering application.

[0079] 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 novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity-based system, characterized in that, Includes the following steps: S1. Construct the mesh connection relationship by building point sets and using connection filtering methods; S2. Based on the fracture-mesh connection relationship constructed in S1, calculate the seepage characteristic parameters of any connection unit. The seepage characteristic parameters include: the conductivity of the connection unit between matrix nodes, the conductivity between two fracture nodes in the fracture layer, and the conductivity between the matrix and the fracture. S3. Based on the seam-mesh connection relationship of S1 and the seepage characteristic parameters obtained in S2, solve the local dual-medium two-phase flow seepage control equation to obtain the gas production rate at different perforation locations.

2. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 1, characterized in that, S1 specifically refers to: S11. Construct a fracture-matrix point set to provide basic data for the subsequent construction of the fracture network connection relationship; S12. Select reasonable connection units based on the influence domain radius and the maximum control angle; S13, based on the crack-matrix point set of S11 and the reasonable connection of S12, forms a crack network connection relationship.

3. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 2, characterized in that, The specific steps for constructing the fracture-matrix point set in S11 are as follows: (1) Crack points: Extracted based on the binarization result of the original crack structure image. A value of 1 represents the crack location and is directly converted into crack point coordinates. (2) Matrix points: The Poisson disk sampling method is used to deploy matrix points to obtain a matrix point set with uniform spatial distribution and appropriate spacing between them; (3) When the Euclidean distance between any matrix point and any crack point is less than the set minimum distance threshold, the matrix point is discarded.

4. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 3, characterized in that, S12 specifically refers to: (1) Establish initial connection units based on the radius of the influence domain: Using any node as the center, set a fixed influence domain radius; traverse all other nodes, and if the Euclidean distance between nodes is less than or equal to the set fixed influence domain radius, then construct an initial connection unit between the two nodes; otherwise, do not connect them. (2) Based on the maximum control angle, reasonable connection units are selected and retained. The selection rules are as follows: Two nodes of any initial connection unit and a third node form a triangle. Calculate the cosine of the included angle between the two sides of the initial connection unit in the triangle. If there is an included angle cosine that is less than the cosine of the maximum control angle, i.e., the included angle is greater than the maximum control angle, then the initial connection unit is discarded. If all included angle cosine values ​​in all the triangles are greater than or equal to the cosine of the maximum control angle, then the initial connection unit is retained as a valid connection unit.

5. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 4, characterized in that, S2 specifically refers to: S21. Based on the principle of conservation of matter, define the node control domain; S22. Based on the two-phase flow control equation of the local dual-medium matrix system, the conductivity of the inter-matrix node connection unit is obtained. S23. Based on the two-phase flow seepage control equation of a local dual-medium fracture system, the conductivity between two fracture nodes in the fracture layer is defined. S24 characterizes the conductivity between the two-phase flow matrix and the crack.

6. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 5, characterized in that, The conductivity of the inter-matrix node connection elements in S22 is shown below: ； in, It is the conductivity of the inter-node connection unit in the matrix; It is a matrix node , Intermediate and average penetration rates; It is a matrix node Control volume; It is the number of matrix nodes; It is a connection unit At the node The volume distribution coefficient in the controlled volume.

7. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 6, characterized in that, The conductivity between two fracture nodes in the fracture layer of S23 is shown below: ； in, It is the conductivity between two fracture nodes in the fracture layer; It is a node and The mobility values ​​between; It is a node in the crack system. , Harmonized average permeability; The thickness of the gas reservoir; For nodes With nodes The average crack aperture; For connection unit The length of the crack.

8. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 7, characterized in that, The conductivity between the matrix and the fracture in S24 is shown below: ； in, It is the conductivity between the matrix and the crack; It is the height of the crack; It is the flow coefficient of the corresponding phase in a two-phase flow; It is a node The harmonized average permeability of the matrix and fracture flow within the control domain; Taken as all nodes The sum of half of the crack connection elements at the endpoints is specifically expressed as: , Represented by node The set of other end nodes of the crack connection unit at the endpoint; Based on nodes and nodes Crack connection unit at the endpoint Length; It is the equivalent normal distance between the crack and the matrix flow, specifically expressed as: .

9. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system as described in claim 8, characterized in that, S3 specifically refers to: S31. Combining the node control domain, the two-phase flow seepage control equations of the local dual-medium matrix system and the local dual-medium fracture system in S2 are discretized; based on the fracture network connection unit constructed in S1 and the conductivity obtained in S2, the discretized local dual-medium two-phase flow seepage control equations are solved to obtain the fluid flow rate of the node connection unit. S32. Based on the obtained conductivity, fluid flow rate and reasonable connection volume of the connecting unit, calculate the gas production rate at different perforation positions.

10. The novel method for simulating seepage in deep coalbed methane fracture networks based on a connectivity element system according to claim 9, characterized in that, The discretized local dual-medium two-phase flow seepage control equations in S31 are as follows: ； in, , They are nodes The pressure; It is to depict nodes Within the control domain, the amount of mass exchanged between matrix nodes and fracture nodes; , They are t , The porosity at time t, where porosity is a function of pressure, is specifically expressed as: ; ; , They are time points , t Time node Phase saturation; This is the penetration rate, taken as the harmonic average, specifically expressed as: , yes Time Node Absolute penetration rate at the location; yes Time Node j Absolute penetration rate at the location; It is the viscosity, taken as the arithmetic mean, specifically expressed as: , yes Phase at time node Viscosity at that point; yes Phase at time node j Viscosity at that point; It is the relative permeability of the gas phase and the aqueous phase, specifically: , It is the equivalent relative permeability of the corresponding phase on the connecting unit; It is a moment t arrive Nodes within the time period Water saturation; It is a moment t arrive During the time period, the connecting unit , j The pressure of the corresponding phases; It is a moment t arrive During the time period, the connecting unit The pressure corresponding to the upper phase.

Images