A Meshless Numerical Modeling Method for Fractured Reservoirs Based on Extended Finite Volume Method

Rao, Xiang; Zhao, Hui; Liu, Yina

doi:10.2118/210581-PA

Summary

In this paper, a meshless numerical modeling method named mesh-free discrete fracture model (MFDFM) of fractured reservoirs based on the newly developed extended finite volume method (EFVM) is proposed. First, matching and nonmatching point cloud generation algorithms are developed to discretize the reservoir domain with fracture networks, which avoid the gridding challenges of the reservoir domain in traditional mesh-based methods. Then, taking oil/water two-phase flow in fractured reservoirs as an example, MFDFM derives the EFVM discrete scheme of the governing equations, constructs various types of connections between matrix nodes and fracture nodes, and calculates the corresponding transmissibilities. Finally, the EFVM discrete scheme of the governing equations and the generalized finite difference discrete scheme of various boundary conditions form the global nonlinear equations, which do not increase the degree of nonlinearity compared with those in the traditional finite volume method (FVM)-based numerical simulator. The global equations can be solved by the existing nonlinear solver in the FVM-based reservoir numerical simulator by only adding the linear discrete equations of boundary conditions, which reduce the difficulty of forming a general purpose MFDFM-based fractured reservoir numerical simulator. Several numerical test cases are implemented to illustrate that the proposed MFDFM can achieve good computational performance under matching and nonmatching point clouds, and for heterogeneous reservoirs, complex fracture networks, complex boundary geometry, and complex boundary conditions, by comparing the computational results of MFDFM with embedded discrete fracture model (EDFM). Thus, MFDFM retains the computational performances of the traditional mesh-based methods and can avoid the difficulties of handling complex geometry and complex boundary conditions of the computational domain, which is the first meshless numerical framework to model fractured reservoirs in parallel with the mesh-based discrete fracture model (DFM) and EDFM.

Introduction

In recent years, great breakthroughs have been made in the exploration and development of unconventional oil and gas reservoirs in the world. Horizontal well fracturing technology is the core technology of unconventional oil and gas reservoir development and has been widely used in the development of low permeability/tight oil and gas, shale oil and gas reservoirs. Efficient and accurate fracturing production evaluation and dynamic analysis technology have an important application background. Numerical simulation is the most widely used production prediction method for the fractured reservoir. At present, the numerical simulation technology of fractured oil and gas reservoirs is mainly divided into two categories: the continuous medium model and the discrete fracture model. Among them, the continuous medium model mainly includes the equivalent continuous medium model (Snow 1968) and the dual or more medium model (Warren and Root 1963; Pruess and Narasimhan 1985; Wu et al. 2004). Due to the large differences in height, length, opening, and distribution of reservoir fractures in space, the uniformity assumption of the continuous medium model is significantly different from reality. To describe the fracture geometry more accurately, some models explicitly characterize the large-scale fractures to obtain high-accuracy results by reducing the dimension of discrete fractures, mainly including the DFM based on a matching grid and EDFM based on a structured grid.

DFM based on a matching grid generates unstructured grids to match the geometry of large-scale fractures, so that the fractures are located at the intersection of adjacent grids (Slough et al. 1999). At present, many scholars have developed a variety of DFMs based on different numerical methods to simulate the single-phase and multiphase flow problems in fractured media, including:

Galerkin finite element method (FEM): Noorishad and Mehran (1982) used the Galerkin FEM based on the variational principle to simulate the solute dispersion and convection in the fracture media. Karimi-Fard and Firoozabadi (2001) proposed an FEM-based two-phase flow model considering the gravity and capillary force. However, in the multiphase flow problem, due to the difference between the phase permeability curves of the matrix and the fracture, the corresponding node values of the saturation of the two media are not equal at the common node of the fracture and the matrix. However, in the Galerkin FEM, the degree of freedom of the fracture and the matrix at the common node needs to be consistent, so the method does not provide the local mass conservation.
Block-center FVM: Karimi-Fard et al. (2004) coupled the two-phase flow equations based on FVM, which satisfies the local mass conservation.
Control volume FEM: Compared with FEM, control volume FEM satisfies the local mass conservation and does not need to calculate the interflow between the matrix grid and the fracture grid (Monteagudo and Firoozabadi 2004; Matthäi et al. 2005).
Mixed FEM is based on the approximation scheme of flow on grid edge expressed by grid-center pressure and grid-edge pressure in the case of the unstructured grid and combines the conditions of continuous flow and continuous pressure to offset the flow term in the global equations, thus forming a closed equation system to solve (Chavent and Roberts 1991; Zhang et al. 2017a). Zhang et al. (2017b) proposed DFM based on multiscale hybrid FEM, which further improved the computational efficiency of hybrid FEM-based models.
Boundary element methods (BEMs): Fang et al. (2017) used the boundary integral equation to obtain the approximate scheme of the mass transfer between the matrix grid and the fracture grid and then couple it into the finite element discrete format of the matrix system and the fracture system to form a global equation system for solution. BEM is based on the fundamental solution of differential equation that has a semianalytical accuracy. However, because the high accuracy of BEM is generally ensured when handling elliptic differential equations that always have the fundamental solutions (Bogomolny 1985; Cheng and Cheng 2005; Salam et al. 2017; Rao et al. 2019a), such as the Laplace equation of single-phase steady flow and the Helmholtz equation from Laplace transformation of single-phase transient flow equation, BEM is generally used for the study of single-phase flow problems and is difficult to be applied to hyperbolic multiphase flow problems with strong convection.

Although the above numerical methods are different, resulting in differences in the accuracy and efficiency of various DFMs, it is very difficult to generate a high-quality matching grid due to the complex geometry of the actual reservoir fracture network, and the matching grid needs to be updated continuously if the dynamic behaviors of the fracture network are considered, will significantly increase the computational cost.

To avoid the difficulty of generating a high-quality matching unstructured grid in DFM, the EDFM handles the explicitly characterized fractures by dividing the reservoir matrix into a structured grid and embedding the discrete fractures into the structured matrix grid as the source or sink term of the matrix grid. Li and Lee (2008) proposed EDFM for the first time and simulated the multiphase flow equations by using FVM. Moinfar et al. (2013) established 3D EDFM for compositional reservoir model for the first time and analyzed the influence of natural fracture and hydraulic fractures on flow field and production of horizontal fractured wells. In the practical application of the large-scale field, the required large number of grids in EDFM will significantly increase the computational cost. Multiscale methods were used to improve the computation efficiency in this case. For example, Hajibeygi et al. (2011) and Ţene et al. (2016) applied iterative multiscale FVM and algebraic multiscale methods to EDFM, respectively. Because the flow between matrix grid and fracture grid in EDFM generally adopts a two-point linear approximation scheme, the calculation accuracy of EDFM for the early production performance of fractured horizontal wells is low. In this regard, Rao et al. (2019b) added multilayer virtual grids on both sides of the fracture to characterize the mass transfer between matrix and fracture more accurately and to improve the calculation accuracy of EDFM. Wang et al. (2020) constructed the unsteady transmissibility between the fracture grid and matrix grid to improve the early calculation accuracy. In addition to the commonly used EDFM based on FVM, Yan et al. (2016) and Zhang et al. (2017c) established the mimetic finite difference method and multiscale mimetic finite difference methodbased EDFMs, respectively. Wu et al. (2020) used the Green element method based on boundary integral equation to deal with the single-phase flow equation after Laplace transform and combined with Stephest numerical inversion of Laplace transform to realize the high-accuracy transient analysis of single-phase flow, but this method can only be applied to single-phase flow. Ţene et al. (2017) pointed out that EDFM cannot effectively deal with low-conductivity “fractures” (e.g., interlayer, fault, etc.) placed in the reservoir domain and proposed to project the fracture grid to the matrix grid interface to construct additional fracture-matrix connection and establish projection-based EDFM (pEDFM) to solve this problem. Jiang and Younis (2017) pointed out that EDFM will induce significant error in the simulation of multiphase flow across fractures and used an improved pEDFM to effectively resolve this limitation. Rao et al. (2020) improved the transmissibility formula of connections between grids, established the connection between fracture grids projected on the interface of the same matrix grid, and eliminated the significant errors of the original pEDFM and EDFM in the calculation of a class of reservoir flow problems. Olorode et al. (2020) constructed 3D pEDFM for compositional reservoir models. Rao and Liu (2022) pointed out that current pEDFM implementations still cannot effectively handle many cases of flow barriers and demonstrate that a necessary condition for a fracture projection configuration to be physical in pEDFM is that the fracture projection configuration is geometrically homeomorphic to the fracture. Up to now, EDFM family methods (i.e., EDFM and pEDFM) have been widely used in the numerical simulation of fractured oil and gas reservoirs (Li et al. 2018; Dachanuwattana et al. 2018; Shakiba et al. 2018; Rao et al. 2022b; HosseiniMehr et al. 2022). For example, Yu et al. (2018) applied EDFM to the flow simulation of complex fracture networks of shale gas reservoirs, and Ren et al. (2018) and Ren and Younis (2021) proposed coupling extended finite element method and pEDFM to study the dynamic behaviors of fractures in multiphase flow. Based on extended finite element method and EDFM, Yan et al. (2018) improved the simulation accuracy of the dynamic behavior of natural and induced fractures by introducing the MINC (multiple interacting continua) model. Li and Voskov (2021) developed a multilevel discrete fracture model (MLDFM) by adopting a triple continuum model considering the effect of flow direction on the complex local-scale flow response.

Although the EDFM-based methods avoid the difficulty of generating high-quality matching grids, the commonly used EDFM based on the Cartesian grid is generally difficult to obtain efficient and effective treatment of complex reservoir geometry and complex boundary conditions (such as Dirichlet boundary at constant pressure, nonclosed Neumann boundary conditions, etc.). For the discretization of the reservoir domain with complex geometry, Karimi-Fard et al. (2006) presented a systematic multisubregion upscaling approach, which could be regarded as an EDFM based on unstructured grid. Then, Karimi-Fard and Durlofsky (2016) developed a general gridding, discretization, and coarsening methodology to make this approach practically adaptive to the complex geological geometry of the reservoir domain. Awadalla and Voskov (2018) applied this approach to tight reservoirs and compared its computational performances with those of traditional EDFM. Besides, Xu et al. (2019) presented a corner-point grids-based EDFM for reservoir complex geological geometry.

Although the unstructured grid in multisubregion and the corner grid can effectively handle the complex reservoir geometry, the generic gridding methods are still complicated. Meshless methods have a great potential to avoid the complexity and difficulty of gridding, but to our knowledge, up to now, there is nearly no study on general purpose meshless flow simulation of fractured reservoirs. Considering that the meshless methods have been widely used in the field of solid mechanics including fracture propagation simulation and have been proved to have the advantage of computational efficiency compared with the grid-based methods (Rao and Rahman 2000; Zhuang et al. 2012), if the general purpose meshless flow simulation method for fractured reservoirs can be realized, it may also significantly facilitate the hot study on coupled geomechanics and multiphase flow in fractured porous media in recent years (Yan et al. 2018; Ren and Younis 2021). Therefore, this paper aims to develop a new numerical modeling framework for the flow simulation of fractured reservoirs, which is based on meshless collocation discretization of the reservoir domain.

The generalized finite difference method (GFDM) is a domain-type meshless method (Benito et al. 2001), which characterizes the calculation domain by a point cloud and divides the domain into intersecting subdomains (called the node influence domain). Using the Taylor series expansion of multivariate functions and the weighted least squares method in each node influence domain, GFDM approximates spatial derivatives of the unknown function by the difference schemes of the unknown function values of nodes in the subdomain, which overcomes the dependence of the traditional FDM on the mesh division. Benito et al. (2003) proposed an h-adaptive GFDM which demonstrates the advantages of meshless methods compared with mesh-based methods for adaptive analysis. Up to now, GFDM has been widely used to solve various scientific and engineering problems (Cheng and Liu 2002; Gavete et al. 2003; Ureña et al. 2012; Li and Fan 2017; Wang et al. 2018; Gu et al. 2019; Fu et al. 2019, 2020b, 2020a; Gu and Sun 2020; Li 2021; Xia and Gu 2021; Chen et al. 2021; Qu and He 2020). Rao (2022a) and Rao et al. (2022a) introduced the first-order single-point upstream (SPU) scheme in FDM/FVM-based reservoir numerical simulator to GFDM, simulated the convection-dominated heat and mass transfer and two-phase porous flow problem (Edwards 2011; Li and Tchelepi 2015). The simulation results illustrated GFDM could achieve good computational performances, indicating the application potentials of meshless GFDM in the reservoir numerical simulation. Then, to handle the singular source or sink terms resulted by wells in GFDM framework, Rao (2022b) defined and calculated the node control volume for GFDM in the case of the connectable point cloud to form a node control domain-based meshless method. Node control domain-based meshless method can be regarded as an EFVM which extended traditional low-order FVM to a meshless version. Rao (2022b) took the oil/water two-phase flow as an example to derive the EFVM formulation that provides the local mass conservation. To some extent, EFVM was the integral form of GFDM, and GFDM could be regarded as the differential form of EFVM. Compared with the traditional FVM in reservoir numerical simulation, EFVM can flexibly discretize the computational domain with complex geometry by a point cloud, without mesh division, and can handle more complex boundary conditions. In addition, EFVM can directly use the nonlinear solver in the current FVM-based reservoir numerical simulation technology, such that there is not much technical difficulty to form a general-purpose EFVM-based reservoir numerical simulator, revealing the great application potential of EFVM. In this paper, based on the newly developed EFVM, we present the first meshless numerical modeling method for fractured reservoirs, which can avoid complicated gridding of the reservoir domain with complex geometry (e.g., complex fracture geometry and complex boundary conditions). So this novel method is named mesh-free discrete fracture model and abbreviated MFDFM. The proposed MFDFM based on EFVM provides the first meshless framework for numerical modeling of fractured reservoirs by completely replacing the mesh with the point cloud, which might be a bigger step from the perspective of the reservoir modeling compared with the development of EDFM based on DFM.

Furthermore, since the meshless methods have been confirmed by previous studies to be computationally more efficient in handling solid mechanics calculations compared with mesh-based FEM or extended finite element method, the developed meshless modeling framework may contribute to a better solution of the coupled geomechanics and multiphase flow problems in fractured reservoirs in the future. Therefore, we believe that the work done in this paper may create a new direction of meshless reservoir numerical simulation with multiphysics (e.g., geomechanics) considered and expected to induce a large number of new results in this direction in the future and even form a novel general purpose commercial reservoir numerical simulation software.

This paper is structured as follows: “Basic Governing Equations for Flow in Fractured Reservoirs” takes two-phase flow as an example to introduce studied governing equations. “Principles of GFDM” briefly reviews the meshless GFDM and the latest results on GFDM. “EFVM: Integral Form of GFDM” introduces the EFVM formulation. “The MFDFM Based on EFVM” gives the EFVM-based MFDFM, including developed matching and nonmatching point cloud generation algorithms to discretize the reservoir domain with fracture networks in “Meshless Point Cloud Discretization of the Fracture-Reservoir Domain,” the EFVM-based discrete schemes of the flow governing equations in “Discretization of Governing Equations,” the transmissibility of various connections between the matrix nodes and the fracture nodes in “Definition of Connection Types between Nodes and Transmissibility Calculation,” and the treatment of various boundary conditions by MFDFM in “Treatment of Boundary Conditions and Construction of Global Equations.” In “Numerical Examples,” four numerical examples are given to show the computational performances of the proposed MFDFM. The conclusions and future work are summarized in “Conclusions” and “Future Works,” respectively.

Basic Governing Equations and Numerical Methods

Basic Governing Equations for Flow in Fractured Reservoirs

This paper aims to develop a meshless numerical modeling method for flow simulation in fractured reservoirs that is independent of the currently widely used DFM and EDFM. Next, the paper demonstrates the technical details of the method using a common and representative two-phase flow model as an example and carries out the numerical cases. It should be noted that the meshless modeling method can be extended to three-phase black oil models and compositional models. Furthermore, for simplicity of illustration, this paper focuses on the 2D reservoir models, and the extension of this method to the 3D model is possible with significant work. Because the proposed meshless modeling framework for fractured reservoir MFDFM uses point clouds that are more flexible than grid division to discretize the computing domain, theoretically, the difficulty of MFDFM expanding from 2D to 3D may be less than that of EDFM based on a structured matrix grid and more likely than that of DFM based on a matching grid.

For a 2D two-phase flow model in conventional reservoirs, the conservative equation of phase a is generally expressed as:

,

\nabla \cdot (\frac{k k_{r a}}{B_{a} μ_{a}} \nabla p_{a}) + q_{a, w e l l} = \frac{\partial}{\partial t} (\frac{ϕ S_{a}}{B_{a}})

(1)

where $k$ is the permeability, a = o (i.e., oil phase) or w (i.e., water phase), $k_{r a} = k_{r a} (S_{w})$ is the relative permeability of phase a, $B_{a}$ is the volume coefficient of phase a, $μ_{a}$ is the viscosity of phase a; $p_{a}$ is the pressure of phase a; $t$ is the time; $S_{a}$ is the saturation of phase a; $ϕ = ϕ (p)$ is the reservoir porosity which is the function of pressure. $q_{a, w e l l}$ is the strength of the source or sink terms of phase a representing the production well or injection well.

For fractured reservoirs, where there is a fracture system in addition to the matrix system, the phase-a conservative equation is rewritten as

,

\nabla \cdot (\frac{k k_{r a}}{B_{a} μ_{a}} \nabla p_{a}) + q_{a, m f} + q_{a, w e l l} = \frac{\partial}{\partial t} (\frac{ϕ S_{a}}{B_{a}})

(2)

in which, similar to Eq. 1, a = o (i.e., oil phase) or w (i.e., water phase) and other physical parameters. Different from Eq. 16, the added $q_{a, m f}$ is the quantity of mass transfer between matrix and fracture of phase a.

Principles of GFDM

GFDM is a meshless method, which only needs to generate a point cloud (the point in a point cloud is also called the node in this paper) to discretize the reservoir computational domain. Compared with the mesh generation of the reservoir domain, node collocation in the point cloud has more flexibility. GFDM defines a node influence domain for each node. The influence domains of different nodes may intersect, that is as follows:

,

Ω = ⋃_{i} I_{i}, I_{i} | I_{j} \neq Φ

(3)

where $Ω$ is the calculation domain, and $I_{i}$ is the influence domain of each node.

As shown in Fig. 1, let node i be a node in the calculation domain, and the coordinate is $(x_{i}, y_{i})$ ⁠. There are $n_{i}$ other nodes in the influence domain of node i (i.e., $n_{i}$ is the count of the j’s), and the corresponding coordinates of these nodes are $(x_{j}, y_{j})$ (j = 1, 2, 3, ... ,n_i). The Taylor expansion of the unknown function at node i obtains:

u (x_{j}, y_{j}) = u (x_{i}, y_{i}) + Δ x_{j} {\frac{\partial u}{\partial x} |}_{i} + Δ y_{j} {\frac{\partial u}{\partial y} |}_{i} + \frac{1}{2} ({(Δ x_{j})}^{2} {\frac{\partial^{2} u}{\partial x^{2}} |}_{i} + 2 Δ x_{j} Δ y_{j} {\frac{\partial^{2} u}{\partial x \partial y} |}_{i} + {(Δ y_{j})}^{2} {\frac{\partial^{2} u}{\partial y^{2}} |}_{i}) + O (r^{3}),

(4)

Fig. 1

View large Download slide

Sketch of nodes and the node influence domain, in which the gray circular area is the influence domain of node i, and there are eight other nodes in the influence domain, that is, $n_{i}$ = 8.

where $Δ x_{j} = x_{j} - x_{i}$ ⁠, $Δ y_{j} = y_{j} - y_{i}$ ⁠.

Suppose $u_{i} = (x_{i}, y_{i})$ ⁠, $u_{x, i} = {\frac{\partial u}{\partial x} |}_{i}$ ⁠, $u_{y, i} = {\frac{\partial u}{\partial y} |}_{i}$ ⁠, $u_{x x, i} = {\frac{\partial^{2} u}{\partial x^{2}} |}_{i}$ ⁠, $u_{x y, i} = {\frac{\partial^{2} u}{\partial x \partial y} |}_{i}$ ⁠, $u_{y y, i} = {\frac{\partial^{2} u}{\partial y^{2}} |}_{i}$ ⁠.

Define the weighted error function $B (D_{u})$ ⁠:

,

B (D_{u}) = \sum_{j = 1}^{n} {[(u_{i} - u_{j} + Δ x_{j} u_{x, i} + Δ y_{j} u_{y, i} + \frac{1}{2} {(Δ x_{j})}^{2} u_{x x, i} + \frac{1}{2} {(Δ y_{j})}^{2} u_{y y, i} + Δ x_{j} Δ y_{j} u_{x y, i}) ω_{j}]}^{2}

(5)

where $D_{u} = {(u_{x, i}, u_{y, i}, u_{x x, i}, u_{y y, i}, u_{x y, i})}^{T}$ ⁠, $ω (Δ x_{j}, Δ y_{j})$ is the weight function. The purpose of adding the weight function in $B (D_{u})$ is to make the nodes closer to the central node have a lower truncation error in the Taylor expansion after minimizing $B (D_{u})$ ⁠, which requires the nodes closer to the central node to have larger weights, that is, the weight function is required to be monotonically decreasing with respect to the distance to the central node (Liu and Gu 2005). Commonly used weight functions include exponential functions, potential functions, cubic splines, quartic splines, etc. (Liu and Gu 2005). Benito et al. (2001), Gavete et al. (2003), and Chen et al. (2021) demonstrated that there is no significant difference in the computational accuracy of these weight functions when GFDM is applied to general second-order equations. Since the quadratic spline function has third-order continuity, and as can be seen from Eq. 4, the highest-order derivative term of the unknown function in the Taylor expansion in GFDM is only the second-order spatial derivatives. Many problems we study also require only the existence of the second-order derivative of the unknown function, so to avoid the added weight functions reducing the continuity order of the original problem, quadratic spline functions in Eq. 6 with third-order continuity are widely used in the study of GFDM (Gavete et al. 2003).

,

ω_{j} = {\begin{array}{cc} 1 - 6 {(\frac{r_{j}}{r_{m}})}^{2} + 8 {(\frac{r_{j}}{r_{m}})}^{3} - 3 {(\frac{r_{j}}{r_{m}})}^{4} & r_{j} \leq r_{m} \\ 0 & r_{j} > r_{m} \end{array}

(6)

where $r_{j}$ is the Euclidean distance from node j to central node i, and $r_{m}$ is the radius of the influence domain of node i.

To minimize the error function, the first-order derivative of $B (D_{u})$ to each component of $D_{u}$ is equal to zero:

.

\frac{\partial B (D_{u})}{\partial u_{x, i}} = 0, \frac{\partial B (D_{u})}{\partial u_{y, i}} = 0, \frac{\partial B (D_{u})}{\partial u_{x x, i}} = 0, \frac{\partial B (D_{u})}{\partial u_{y y, i}} = 0, \frac{\partial B (D_{u})}{\partial u_{x y, i}} = 0

(7)

The above equations are sorted into linear equations as follows:

,

A = L^{T} ω L, b = L^{T} ω U, L = {({L_{1}}^{T}, {L_{2}}^{T}, L, {L_{n}}^{T})}^{T}, L_{i} = (Δ x_{i}, Δ y_{i}, \frac{Δ {x_{i}}^{2}}{2}, \frac{Δ {y_{i}}^{2}}{2}, Δ x_{i} Δ y_{i}), ω = d i a g ({ω_{1}}^{2}, {ω_{2}}^{2}, L, {ω_{n}}^{2}), U = {(u_{1} - u_{i}, u_{2} - u_{i}, L, u_{n} - u_{i})}^{T} .

(8)

where $A = L^{T} ω L, b = L^{T} ω U, L = {({L_{1}}^{T}, {L_{2}}^{T}, L, {L_{n}}^{T})}^{T}, L_{i} = (Δ x_{i}, Δ y_{i}, \frac{Δ {x_{i}}^{2}}{2}, \frac{Δ {y_{i}}^{2}}{2}, Δ x_{i} Δ y_{i}),$

$ω = d i a g ({ω_{1}}^{2}, {ω_{2}}^{2}, L, {ω_{n}}^{2}), U = {(u_{1} - u_{i}, u_{2} - u_{i}, L, u_{n} - u_{i})}^{T}$ ⁠.

From Eq. 8, it can be obtained that:

,

D_{u} = {(u_{x, i}, u_{y, i}, u_{x x, i}, u_{y y, i}, u_{x y, i})}^{T} = A^{- 1} b = A^{- 1} L^{T} ω U = M U

(9)

where $M = A^{- 1} L^{T} ω$ ⁠. The matrix M has five rows, and the number of columns is equal to the total number of nodes in the domain of influence of node i except node i, that is, $n_{i}$ ⁠. Each row of the matrix M denotes the difference coefficients in the generalized difference expression for each spatial derivative of the unknown function at node i in $D_{u}$ ⁠, respectively. For clarity of description, the elements of the matrix $M$ are denoted as $m_{k j}^{i}$ ⁠, where $k = 1, 2, \dots, 5$ ⁠, $j = 1, 2, \dots, n_{i}$ ⁠, and the superscript i represents that the matrix $M$ is obtained when considering the influence domain of node i.

From the matrix $M$ ⁠, it is not difficult to obtain the difference expressions of the spatial derivatives at node i:

\begin{matrix} u_{x, i} = \sum_{j = 1}^{n_{i}} m_{1 j}^{i} (u_{j} - u_{i}), u_{y, i} = \sum_{j = 1}^{n} m_{2 j}^{i} (u_{j} - u_{i}), u_{x x, i} = \sum_{j = 1}^{n} m_{3 j}^{i} (u_{j} - u_{i}), u_{y y, i} = \sum_{j = 1}^{n} m_{4 j}^{i} (u_{j} - u_{i}), \\ u_{x y, i} = \sum_{j = 1}^{n} m_{5 j}^{i} (u_{j} - u_{i}) \end{matrix} .

(10)

It can be seen that GFDM only relies on coordinates of nodes in the influence domain of the considered node to approximate the spatial derivatives at the node, so it can realize the discretization of the differential equations without mesh division, which is a significant advantage of meshless GFDM compared with mesh-based FDM, FEM and FVM.

EFVM: Integral Form of GFDM

In reservoir numerical simulation, it is difficult to avoid the existence of a singular point source or point sink term (i.e., Dirac term) representing injection or production wells. However, classical GFDM does not use any integration treatment and does not have the concept of the control volume, so it is difficult to integrate the governing equation like FEM and FVM, to nonsingularize the singular Dirac terms. It is also difficult to approximate the strength of the singular source or sink to the strength of unit grid volume as FDM. Therefore, Rao (2022b) introduced the idea of node control volumes and calculated them for GFDM. Taking oil/water two-phase flow as an example, the extended finite volume discrete scheme satisfying local mass conservation is derived. EFVM can directly use the nonlinear solver in the widely used FVM-based reservoir numerical simulation technology, so it can easily form a general purpose EFVM-based reservoir numerical simulator.

Definition of the Node Control Volume

In addition to its influence domain, it is defined that each node has its control domain, and the control domains of distinct nodes do not intersect. Thus,

,

Ω = ⋃_{i} Ω_{i}, Ω_{i} | Ω_{j} = \emptyset

(11)

,

\sum_{i = 1}^{n_{p}} {\bar{V}}_{i} = V_{Ω}

(12)

where $Ω_{i}$ is the control domain of node i, ${\bar{V}}_{i}$ is the volume of $Ω_{i}$ ⁠, named the control volume of node i, ${\bar{V}}_{Ω}$ is the volume of the calculation domain $Ω$ ⁠, and $n_{p}$ is the total number of nodes in the calculation domain. Under the meshless framework, this study does not concern with determining the precise geometric shape of the node control domain $Ω_{i}$ due to possible complex node collocation in the computational domain. Instead, it focuses solely on calculating the appropriate node control volume.

Then, a connectable point cloud is defined as follows: for any two nodes in the computational domain (may as well be denoted as node i and node j), there must also be node i in the influence domain of node j if there is node j in the influence domain of node i. Thus, in case of a connectable point cloud, if node j is in the influence domain of node i, we denote that node i is neighboring to j, and j is also neighboring to i, and (i, j) is denoted as a neighboring pair. In practice, satisfying a connectable point cloud of the computational domain is a straightforward task. For example, simply setting the radius of each node’s influence domain (when the shape of the influence domain is selected as circular, of course, other shapes are also allowed) to be the same can be sufficient to get the desired connectable point cloud. However, there is still a valuable future work mentioned in “Future Works” on the construction of the connectable point cloud that has strongly heterogeneous node density on the calculation domain, for example, some of the matching point clouds in “Meshless Point Cloud Discretization of the Fracture-Reservoir Domain.”

It is worth noting that, there can be no concept of the node influence domain, because it is only necessary to determine which nodes participate in the construction of the generalized finite difference expressions of spatial derivatives at the considered node, and the introduction of the node influence domain just makes it easy and direct to determine the selection of these nodes. The set of these nodes is also named the local point cloud of the considered node.

Calculation of Node Control Volumes by Adding Virtual Nodes

For a boundary node, having nodes in its influence domain only on the inner side of the boundary node causes the center of gravity of the local point cloud of the boundary node to deviate significantly from the boundary node, resulting in the reduction of the accuracy of the generalized difference expressions of spatial derivatives in Eq. 10 obtained at the boundary node. The addition of a virtual node outside the boundary node can bring the center of gravity of the local point cloud in the domain of influence of the boundary node closer to the boundary node, thus improving the accuracy of the generalized difference approximation (Milewski 2012,Rao et al. 2022a). This is why we point out the need to add virtual nodes outside the boundary when calculating the nodal control volume in this section and the treatment of boundary conditions in “Treatment of Boundary Conditions and Construction of Global Equations.” This is similar to the treatments in the traditional FDM, where at a boundary node in the node-center FDM if that boundary point satisfies the boundary condition expressed by the first-order spatial derivative, there is only first-order accuracy if the boundary condition is approximated by a one-side finite difference scheme without adding a virtual node, but second-order accuracy if the boundary condition is approximated by a central difference with the addition of a virtual node. Readers can refer to Rao (2022b) for relevant details. After adding virtual nodes, Rao (2022b) derived the following equation for node control volumes:

,

V_{i} (m_{3 j}^{i} + m_{4 j}^{i}) = V_{j} (m_{3 i}^{j} + m_{4 i}^{j})

(13)

where $V_{i}$ is the control volume of each node after adding virtual nodes on the boundary, rather than the node control volume ${\bar{V}}_{i}$ only inside the calculation domain. The concept of the characteristic angle at each node is introduced to establish the relationship between $V_{i}$ and ${\bar{V}}_{i}$ ⁠. The characteristic angle at the node is identical to that defined in the BEM. The characteristic angle of a smooth node on the boundary is $π$ ⁠, the corner node on the boundary is the included angle $θ$ of tangents on both sides, and the internal node in the computation domain is $2 π$ ⁠. As a result, if denote the characteristic angle of the node i as $θ_{i}$ ⁠, then:

.

{\bar{V}}_{i} = \frac{θ_{i}}{2 π} V_{i}

(14)

Thus, Eq. 12 is rewritten as follows:

.

\sum_{i = 1}^{n_{p}} \frac{θ_{i}}{2 π} V_{i} = V_{Ω}

(15)

Finally, if there is $n_{c}$ pair of i-j, then each pair of i-j has an Eq. 13. Combined with Eq. 15, $n_{c}$ + 1 linear equations can be obtained. There are $n_{p}$ unknown node control volumes, and $n_{c} > n_{p}$ is generally true. Therefore, Eq. 13 and Eq. 15 form an overdetermined linear equation system, which can be solved by the least squares method, and then use Eq. 14 to calculate ${\bar{V}}_{i}$ ⁠. It is important to note that there is a premise in the above derivation process, that is, the point cloud of the computing domain is connectable. In addition, there is another detail to explain. Since the overdetermined system of equations is solved by least squares method in two-norm, Eq. 15 may not be strictly satisfied. However, the sum of the node control volumes is equal to the total volume of the computational domain, which is the basic measure of the mass conservation of the computational domain. Therefore, the “large number method” is used to multiply both sides of Eq. 15 by the larger number G and then obtained Eq. 16 to replace Eq. 15 in the overdetermined linear system of equations, to make Eq. 15 to be satisfied as strictly as possible.

.

G \sum_{i = 1}^{n_{p}} \frac{θ_{i}}{2 π} V_{i} = G V_{Ω}

(16)

Rao (2022b) showed that the convergence order of GFDM is less than or equal to that of the traditional low-order FDM (i.e., the convergence order is less than or equal to 2), so GFDM is also a low-order method. When ignoring the influence of the first-order single-point upstream scheme on the convergence order of the discrete scheme, the FVM commonly used in reservoir numerical simulation has second-order accuracy for the diffusion term expressed by the second-order spatial derivative. Thus it is a low-order FVM, which can be regarded as the integral form of the low-order FDM to a certain extent. Of course, many high-order FVMs have been developed in computational fluid dynamics, but they are rarely applied to reservoir numerical simulation. EFVM, the integral form of GFDM, has the same approximation accuracy of diffusion term as GFDM, which will be less than or equal to the second order, so EFVM is a low-order method. When the meshless collocation in EFVM is just composed of the grid center points in FVM, with selecting a small-radius node influence domain EFVM will reduce to the traditional low-order FVM. Since the complex fracture network in fractured reservoirs generally has complex geometries, the meshless EFVM allows for the simpler and more flexible discretization of the fractured reservoir computational domain using point clouds than the mesh-based methods.

It should be noted that when calculating the node control volume, the radius of the node influence domain cannot be too large, because when the node influence domain is too large, there will be many nodes in the influence domain, which will bring three disadvantages to the calculation of the node control volume: (i) the dimension of the overdetermined equations is large, which reduces the calculation efficiency; (ii) the connection between nodes is more complex, resulting in a large bandwidth of the Jacobian matrix in the nonlinear solver, which will significantly reduce the efficiency of constructing and calculating equations; and (iii) in recent work on upwind GFDM-based heat and mass transfer coupling calculation (Rao 2022a) and two-phase flow simulation (Rao et al. 2022a), it is found that for hyperbolic-like porous flow problems with relatively strong convection (such as multiphase flow with low capillary, heat and mass transfer, etc.), the increase of the radius of the node influence domain will increase the GFDM calculation error for the convection-controlled physical quantities (such as water saturation, temperature, etc.). As the integral form of GFDM, when the radius of the node influence domain is large, the EFVM-based calculation accuracy for these physical quantities may also be reduced.

Discretization of Governing Equations and Definition of Transmissibility

Taking the two-phase flow equation in Eq. 1 as an example, Rao (2022b) integrated the GFDM discrete scheme of the two-phase flow equation at node i with the node control volume $V_{i}$ ⁠, to obtain the EFVM discrete scheme in Eq. 17. For the specific derivation process, readers can refer to Rao (2022b).

,

\sum_{j = 1}^{n_{i}} [\frac{k_{r a, i j}^{t + Δ t}}{B_{a, i j}^{t + Δ t} μ_{a, i j}^{t + Δ t}} T_{i j}^{t + Δ t} (p_{a, (i, j)}^{t + Δ t} - p_{a, i}^{t + Δ t})] + Q_{a, w e l l} = \frac{V_{i}}{Δ t} [(ϕ_{i}^{t + Δ t} S_{a, i}^{t + Δ t}) - (ϕ_{i}^{t} S_{a, i}^{t})]

(17)

where $T_{i j}^{t + Δ t}$ is the transmissibility between node i and node j and is calculated by:

,

T_{i j}^{t + Δ t} = k_{i j}^{t + Δ t} V_{i} (m_{3, j}^{i} + m_{4, j}^{i}) = k_{j i}^{t + Δ t} V_{j} (m_{3, i}^{j} + m_{4, i}^{j})

(18)

in which, $m 3, j i$ and $m 4, j i$ are the elements of matrix M calculated according to Eq. 7 when node i is the central node. $m 3, i j$ and $m 4, i j$ are the elements of matrix M calculated according to Eq. 7 when node j is the central node.

,

\begin{array}{l} k_{r a, i j}^{t + Δ t} = {\begin{cases} k_{r a, i} (S_{w, i}^{t + Δ t}) i f p_{a, (i, j)}^{t + Δ t} \geq p_{a, i}^{t + Δ t} \\ k_{r a, i} (S_{w, i}^{t + Δ t}) i f p_{a, (i, j)}^{t + Δ t} < p_{a, i}^{t + Δ t} \end{cases} \end{array}

(19)

{k_{i j}}^{t + Δ t} = \frac{2}{1 / {k_{i}}^{t + Δ t} + 1 / {k_{j}}^{t + Δ t}}, {B_{a, i j}}^{t + Δ t} = \frac{{B_{a, i}}^{t + Δ t} + {B_{a, j}}^{t + Δ t}}{2}, {μ_{a, i j}}^{t + Δ t} = \frac{{μ_{a, i}}^{t + Δ t} + {μ_{a, j}}^{t + Δ t}}{2} .

(20)

It should be noted that EFVM can obtain the transmissibility satisfying local mass conservation between node i and node j without mesh division in the computational domain. In addition, EFVM can deal with various boundary conditions, especially the mixed boundary conditions that are difficult to deal with by traditional FVM. This part about the treatment of boundary conditions will be illustrated in detail in “Treatment of Boundary Conditions and Construction of Global Equations.”

The MFDFM Based on EFVM

Meshless Point Cloud Discretization of the Fracture-Reservoir Domain

MFDFM is based on EFVM, so it is essentially a meshless method. The meshless method requires only a point cloud to discretize the computational domain, and the point cloud contains much less topological information than the topological information in the mesh (e.g., based on the point cloud the mesh needs to determine which two nodes join to form an edge, which nodes form a mesh, the order of the nodes when forming the mesh, and the fact that the lengths of the edges of a mesh in the computational domain generally do not differ much, resulting in the small size meshes required in narrow areas of the computational domain needing to pass through more neighboring meshes before they become normal-sized, making the matching mesh will have a large number of small size meshes in narrow areas of the computational domain). The generation of a point cloud to discretize the computational domain is simpler, more flexible, and more diverse than mesh division, as the topological constraints on the generation of point clouds are much less than those on mesh division (Milewski 2012; Michel et al. 2021), which gives the meshless method an advantage over mesh-based methods in terms of computational domain discretization and adaptive analysis, which is one of the original and most important motivations for the development of meshless methods.

When applying the meshless method to the field of reservoir numerical simulation, fractured reservoirs or fracture-vuggy reservoirs with complex geometry are good application scenarios. For complex fracture networks in the computational domain, small segments/tips of fractures caused by fracture intersections and narrow matrix areas between almost contacting fractures lead to difficulties in generating matching meshes, resulting in many small size meshes at the fracture intersections, small segments/tips, and the narrow matrix regions between the fractures. To discretize the rectangular reservoir domain with a complex fracture network shown in Fig. 2a, when the Delaunay triangulation technique is used to obtain the matching mesh in Fig. 2b, it can be seen that there are many local high-density small size meshes. For such a small reservoir model, the total number of mesh nodes has reached 19,480. Of course, the EDFM can overcome this difficulty by simply generating a structured mesh division of the matrix domain of the computational domain and embedding the fracture network to the background matrix mesh, which can be called a nonmatching mesh. Similarly, when we use point clouds to discretize the reservoir computational domain, the point clouds can also be classified as matching or nonmatching depending on how the fracture network is treated. A matching point cloud is one in which there is a matrix node at the same location as the fracture node, that is, there are two systems at the fracture node, the matrix, and the fracture, and in the case of a two-phase flow simulation, there are four degrees of freedom at the fracture node, that is, the fracture system (node) pressure, the fracture (node) water saturation, the matrix system (node) pressure, and the matrix (node) water saturation. Therefore, the mass transfer between the matrix system and the fracture system based on a matching point cloud is achieved through the mass transfer between the fracture node and the matrix node sharing the same location with the fracture node. The nonmatching point cloud characterizes the mass transfer between the matrix system and the fracture system by determining in the node control domain of which matrix node the fracture node should be located and constructing a connection between the fracture node and that matrix node.

Fig. 2

$Sketches of the reservoir domain and matching Delaunay triangulation. (a) Reservoir domain with a complex fracture network. (b) The matching Delaunay triangulation.$

View large Download slide

Sketches of the reservoir domain and matching Delaunay triangulation. (a) Reservoir domain with a complex fracture network. (b) The matching Delaunay triangulation.

As mentioned earlier, a point cloud contains much less topological information compared with mesh division to discretize the computational domain, leading to the fact that discretizing the computational domain with a point cloud is simpler than the mesh division, and the corresponding point cloud generation methods are more diverse than mesh division methods. In Appendix A and Appendix B, we provide a self-written algorithm for the generation of matching and nonmatching point clouds for the reservoir domain with complex fracture networks, respectively, which will visually demonstrate the simplicity and flexibility of point cloud generation. Source codes for the algorithms are available on GitHub to facilitate the readers to reproduce the generated point clouds (https://github.com/Rao-Xiang/meshless-reservoir-numerical-simulation). From the algorithms in Appendices A and B, it can be seen that:

Similarities between the two algorithms
The matching point cloud generation algorithm is the same as the nonmatching point cloud generation algorithm in Step 1 of fracture intersection calculation and Step 2 of fracture node generation. The resolution of fracture nodes is controlled by the fracture segment characteristic length (denoted as $l_{f, c h a r}$ ⁠).
Differences between the two algorithms:
In the nonmatching point cloud, the nodes on the fracture are only fracture nodes, and the generation of matrix nodes is completely independent of the fracture nodes. Therefore, the nonmatching point cloud generation technique does not need to match the fracture morphology when generating the matrix nodes of the computational domain, so a variety of efficient and simple computational domain point cloud generation techniques can be used to obtain the matrix nodes in the matrix part of the computational domain. At this time, the matrix node distribution is generally uniform, which also ensures the accuracy of the meshless method for the local generalized difference estimation in the node influence domain. Therefore, the nonmatching point cloud generation algorithm directly ignores the fracture network morphology in Step 3 and generates the matrix node distribution in the computational domain. As mentioned earlier, because the point cloud is rarely limited by the topology, it is more flexible than mesh generation when generating matrix nodes in a reservoir domain with complex geometry. In the algorithms of Appendix A and Appendix B, the node characteristic spacing (denoted as $l_{c h a r}$ ⁠) is generally used to control the resolution of matrix nodes. Once the matrix nodes have been generated, Step 4 in the nonmatching algorithm establishes a non-neighboring connection (NNC) between each fracture node and some matrix node, which is determined based on the control volume of the matrix node and information on the distance from the fracture node to each matrix node.
In the matching point cloud, since there are also matrix nodes at the same locations as the fracture nodes (referred to as Type-I matrix nodes in the subsequent description for narrative clarity), the generation of matrix nodes in the matching point cloud as compared with the nonmatching point cloud needs to be limited by the distribution of the fracture nodes generated in Step 2 (i.e., the first class of matrix nodes), that is, the subsequentially generated matrix nodes in Step 3 must not be too close to or overlap with the fracture nodes generated in Step 2. To fully exploit the advantages of point cloud discretization of the computational domain over mesh division, in Step 3 of the matching point cloud generation algorithm, an initial matrix node distribution (referred to as Type-II matrix nodes for clarity) is generated in the same way as the matrix nodes generated in Step 3 of the nonmatching point cloud generation algorithm, and then in Step 4, Type-II matrix nodes that are closer to or overlap with the fracture nodes (i.e., Type-I matrix nodes) are removed in Step 4 to generate the final matrix node distribution.
Processing methods of well nodes
Note that the treatment of well nodes (including vertical well nodes, perforation points of fractured wells, etc.) is not specially pointed out in the point cloud generation algorithm in Appendix A and Appendix B. Here, we point out that the treatment of well nodes can also be divided into matching and nonmatching. The nonmatching treatment is to take the vertical well node as the matrix node closest to the vertical well node according to the principle of nearest distance or take the fracturing perforation point as the fracture node closest to the fracturing perforation point. The simple nonmatching well node treatment can be used in the nonmatching point cloud or the matching point cloud. The matching well node treatment is to take the vertical well node with known coordinates as the preset matrix node or the fracturing perforation point with known coordinates as the preset fracture node. This treatment can only be used in the matching point cloud, accordingly, before generating fracture nodes in Step 2 of the matching point cloud generation algorithm, the fracturing perforation points with known coordinates should be regarded as existing fracture points (i.e., they are included in the set of the fracture intersections calculated in Step 1), and in Step 4, the vertical well nodes with known coordinates should be included in the Type-I matrix nodes as preset matrix nodes, that is, remove the node close to the vertical well node of the known coordinate in the Type-II matrix nodes generated in Step 3. In general, although the matching well-node treatment can make the position of well nodes more in line with the actual position than the nonmatching well node treatment, the nonmatching well node treatment is simpler and more practical than the matching well node treatment.

Overall, whether it is the matching or the nonmatching point cloud generation algorithms, in the first three steps, the steps of generating fracture nodes and matrix nodes are nearly the same, and the generation of fracture nodes is completely independent of the generation of matrix nodes, but the processing of the relationship between matrix nodes and fracture nodes is different in Step 4. The whole process of these two-point cloud generation algorithms is simple. Nonmatching point cloud only needs $l_{c h a r}$ and $l_{f, c h a r}$ (later denoted as control parameters) in Appendix A to control the resolution of a generated point cloud. The nonmatching point cloud generation algorithm can use $l_{c h a r}$ ⁠, $l_{f, c h a r}$ ⁠, and $d i s_{t v}$ in Appendix B to control the resolution of the generated point cloud. Of course, the current algorithms we give may be relatively rough, readers can make relevant improvements to these algorithms or use more efficient or novel point cloud generation technology. After all, the point cloud generation technologies limited by little topology are diverse.

Next, we design a rectangular reservoir model with a complex fracture network and a cardioid reservoir model with an irregular boundary and a complex fracture network, by showing the nonmatching and matching point clouds given by the above two algorithms, to demonstrate that the simplicity and flexibility of generating a point cloud to discretize the computational domain alleviate the gridding challenges for fractured reservoirs. In the numerical tests in “Numerical Examples,” we will illustrate that reliable calculation results can be obtained based on these generated point clouds.

According to the matching point cloud generation algorithm in Appendix B, Figs. 3a through 3c show all nodes, matrix nodes, and fracture nodes in the matching point cloud of the computational domain in Fig. 2a, respectively, when making $l_{c h a r}$ = 10 m, $l_{f, c h a r}$ = 20 m, $d i s_{t v}$ = 5 m. It can be seen that this matching point cloud does not produce a large number of nodes at the fracture intersections and the narrow matrix region between the fractures compared with the matching triangular mesh in Fig. 2b. When different values of the control parameters are selected, for example, Figs. 4 and 5 show the point clouds generated for l_char = 10 m, l_f,char = 10 m, $d i s_{t v} = 3 m$ ⁠, dis_tv = 3 m, and $l_{c h a r}$ = 20 m, $l_{f, c h a r}$ = 20 m, $d i s_{t v}$ = 8 m, respectively; it can be seen that the point clouds in Figs. 4 and 5 are coarser and finer than those in Figs. 3a through 3c respectively, illustrating that the point cloud is flexible in terms of adjusting the point cloud resolution and is simpler to conduct adaptive analysis compared with the matching mesh.

Fig. 3

$Sketches of the matching and nonmatching point clouds to discretize the reservoir domain when setting lchar = 10 m, lf,char = 20 m, and distv = 5 m in the point cloud generation algorithms. (a) All nodes in the matching point cloud, in which the matrix node is marked in blue and the fracture node is marked in black. (b) Matrix nodes in the matching point cloud, in which the Type-I matrix nodes are marked in red and the Type-II matrix nodes are marked in blue. (c) Fracture nodes corresponding to fracture network in the matching point cloud. (d) All nodes in the non-matching point cloud, in which the matrix node is marked in blue and the fracture node is marked in black. (e) Matrix nodes in the non-matching point cloud, in which there is only one type of matrix node. (f) Fracture nodes corresponding to the fracture network in the non-matching point cloud. According to the point cloud generation algorithms, the fracture nodes are the same as those in the matching point cloud.$

View large Download slide

Sketches of the matching and nonmatching point clouds to discretize the reservoir domain when setting l_char = 10 m, $l f, c h a r$ = 20 m, and $d i s_{t v}$ = 5 m in the point cloud generation algorithms. (a) All nodes in the matching point cloud, in which the matrix node is marked in blue and the fracture node is marked in black. (b) Matrix nodes in the matching point cloud, in which the Type-I matrix nodes are marked in red and the Type-II matrix nodes are marked in blue. (c) Fracture nodes corresponding to fracture network in the matching point cloud. (d) All nodes in the non-matching point cloud, in which the matrix node is marked in blue and the fracture node is marked in black. (e) Matrix nodes in the non-matching point cloud, in which there is only one type of matrix node. (f) Fracture nodes corresponding to the fracture network in the non-matching point cloud. According to the point cloud generation algorithms, the fracture nodes are the same as those in the matching point cloud.

Fig. 4

View large Download slide

The matching point cloud when l_char = 10 m, l_f,char = 10 m, dis_tv = 3 m. (a) All nodes. (b) Matrix nodes. (c) Fracture nodes.

Fig. 5

View large Download slide

The matching point cloud when l_char = 20 m, $l_{f, c h a r}$ = 20 m, $d i s_{t v}$ = 8 m. (a) All nodes. (b) Matrix nodes. (c) Fracture nodes.

According to the nonmatching point cloud generation algorithm in Appendix A, Figs. 3d through 3f show all nodes, matrix nodes, and fracture nodes in the nonmatching point cloud in the computational domain of Fig. 1a when supposing $l_{c h a r}$ = 10 m, $l_{f, c h a r}$ = 20 m, respectively. It can be seen that the distribution of matrix and fracture nodes in the nonmatching point cloud is completely independent, and for the rectangular computational domain, a Cartesian matrix point cloud can be used directly.

For the case of irregular reservoir boundary, Figs. 6a and 6b shows the matching triangular mesh of the reservoir calculation domain with a complex fracture network in Fig. 6a. The matching mesh has a large number of small grids at intersections of fractures, small segments/tips, and narrow slivers of matrix between nearly touching fractures, and the total number of mesh nodes reaches 20,458. Figs. 7a through 7c and Figs. 7d through 7f, respectively, show the matching point cloud and nonmatching point cloud given by the algorithms. As mentioned earlier, since the matching point cloud algorithm and nonmatching point cloud algorithm given in the appendixes are independent when generating matrix nodes and fractures nodes in the first three steps; therefore, the mesh nodes of Delaunay triangular mesh in the reservoir calculation domain without considering the complex fracture network (i.e., all nodes in Fig. 7c) can be used as the Type-II matrix nodes in the matching point cloud and the matrix nodes in nonmatching point cloud, respectively. Among them, the number of Type-II matrix nodes in the matching point cloud is 1,814, the number of fracture nodes is 877, and the total number of nodes is 1,814 + 877 × 2 = 3,568. In the nonmatching point cloud, the number of matrix nodes is 2,258, the number of fracture nodes is 877, and the total number of nodes is 2,258 + 877 = 3,135. It can be seen that the total number of nodes in the generated point cloud is far less than the number of mesh nodes in the matching mesh, and the gridding challenge of the matching mesh to the fracture intersection, small fracture tips, and narrow matrix zone in Fig. 6b is avoided.

Fig. 6

$The matching triangulation of the cardioid reservoir domain with complex fracture network. (a) The heart-shaped reservoir domain with complex fracture network. (b) The matching Delaunay triangulation.$

View large Download slide

The matching triangulation of the cardioid reservoir domain with complex fracture network. (a) The heart-shaped reservoir domain with complex fracture network. (b) The matching Delaunay triangulation.

Fig. 7

$The matching and nonmatching point clouds of the cardioid reservoir domain with a complex fracture network. (a) All nodes in the matching point cloud. (b) Matrix nodes in the matching point cloud. (c) Fracture nodes in the matching point cloud. (d) All nodes in the non-matching point cloud. (e) Matrix nodes in the non-matching point cloud. (f) Fracture nodes in the$

View large Download slide

The matching and nonmatching point clouds of the cardioid reservoir domain with a complex fracture network. (a) All nodes in the matching point cloud. (b) Matrix nodes in the matching point cloud. (c) Fracture nodes in the matching point cloud. (d) All nodes in the non-matching point cloud. (e) Matrix nodes in the non-matching point cloud. (f) Fracture nodes in the

Taking the matching point cloud in Fig. 7 as an example, it can be seen that, the point cloud in Fig. 7a of the reservoir domain can be divided into the matrix node layer in Fig. 7b and the fracture node layer in Fig. 7c.

In the matrix node layer, EFVM needs to add virtual nodes outside the calculation domain to solve the node control volume of each real node and deal with the possible derivative boundary conditions. The method of adding virtual nodes is as follows: suppose a boundary node as node i and the external normal vector of node i at the boundary as $n$ ⁠, and the added corresponding virtual node as node j, then the coordinates of virtual node j are as follows:

.

(x_{j}, y_{j}) = (x_{i}, y_{i}) + Δ l n

(21)

These virtual nodes are added to obtain the node control volume of each real node in the calculation domain. Readers can refer to Rao (2022b) for specific details. Denote the calculated node control volume as $V_{i}$ (⁠ $i = 1, 2, 3, \dots, n_{m}$ ⁠), i is the number of real matrix nodes in the matrix node layer, and the control volume of each real node after adding virtual nodes meets Eq. 22, which is a copy of Eq. (13) for the convenience of sequential illutration.

.

V_{i} (m_{3 j}^{i} + m_{4 j}^{i}) = V_{j} (m_{3 i}^{j} + m_{4 i}^{j})

(22)

For the fracture node layer, the fracture is a 1D geometry lower than the reservoir domain. For the 2D reservoir domain shown in Fig. 7c, the fracture is composed of 1D line segments. Therefore, for a fracture node f, the control volume and corresponding surface area of the fracture node can be calculated as:

,

V_{f} = \sum_{g = 1}^{n_{f}} \frac{l_{f g} h_{f g} w_{f g}}{2} A_{f} = \sum_{g = 1}^{n_{f}} \frac{l_{f g} h_{f g}}{2}

(23)

where $n_{f}$ is the number of fracture nodes connected to the fracture node f, $h_{f g}$ and $w_{f g}$ are the average height (in the 2D model, if the fracture penetrates vertically through the reservoir, the average height is the reservoir thickness) and average aperture of the fracture segment, which are calculated in the arithmetic average scheme, that is, $h_{f g} = (h_{f} + h_{g}) / 2$ ⁠, $w_{f g} = (w_{f} + w_{g}) / 2$ ⁠. $l_{f g}$ is the distance between node f and node g. If it is a 3D reservoir calculation domain, the mutually affected fracture nodes are distributed on a 2D surface, and the control volume of each fracture node can be calculated on the fracture surface by using Eq. 10 in 2D GFDM.

Discretization of Governing Equations

For a matrix node i, integrate the node control volume of matrix node i on both sides of the phase-a conservative equation in fractures reservoir shown in Eq. 2 and take the fully implicit scheme to obtain Eq. 24, and detailed derivation process is shown in Appendix C.

.

\begin{array}{l} \sum_{j = 1}^{n_{i}} [\frac{{k_{i j}}^{t + Δ t} {k_{r a, i j}}^{t + Δ t}}{{B_{a, i j}}^{t + Δ t} {μ_{a, i j}}^{t + Δ t}} V_{i} (e_{3, j}^{i} + e_{4, j}^{i}) ({p_{a, (i, j)}}^{t + Δ t} - {p_{a, i}}^{t + Δ t})] + f l a g (i) \frac{{k_{r a, i f}}^{t + Δ t} {T_{i f}}^{t + Δ t}}{{B_{a, i f}}^{t + Δ t} {μ_{a, i f}}^{t + Δ t}} ({p_{a, f}}^{t + Δ t} - {p_{a, i}}^{t + Δ t}) \\ + Q_{a, w e l l} = \frac{V_{i}}{Δ t} [(\frac{{ϕ_{i}}^{t + Δ t} {S_{a, i}}^{t + Δ t}}{{B_{a, i}}^{t + Δ t}}) - (\frac{{ϕ_{i}}^{t} {S_{a, i}}^{t}}{{B_{a, i}}^{t}})] \end{array}

(24)

The specific meaning of each symbol in Eq. 24 can be seen in Appendix. It should be emphasized here that, $n_{i}$ is the number of nodes in the influence domain of matrix node i, $(i, j)$ refers to the sequence number of the j-th node in the global node except node i itself in the influence domain of matrix node i. Virtual nodes may exist in the influence domain of matrix node i at the boundary, but the role of virtual nodes in MFDFM is to assist in calculating the control volume of real nodes and dealing with derivative boundary conditions.

The transmissibility between matrix nodes can be extracted from the first items on the left hand side of Eq. 24, and this part is similar to that in the EFVM. When matrix node i is taken as the central node, the transmissibility between matrix node i and matrix node j (in fact, it should be the serial number of the j-th node in the influence domain of node i in the global node, which is denoted as j here for simplicity of expression) is as follows:

,

T_{i, i j}^{t + Δ t} = [k_{i j}^{t + Δ t} V_{i} (e_{3, j}^{i} + e_{4, j}^{i})]

(25)

where i in the subscript $i, i j$ denotes that the transmissibility between node i and node j is obtained from the perspective of the control volume of node i.

Similarly, from the perspective of the control volume of the node $(i, j)$ ⁠, the following is obtained as:

.

T_{j, i j}^{t + Δ t} = [k_{j i}^{t + Δ t} V_{j} (e_{3, i}^{j} + e_{4, i}^{j})]

(26)

The permeability $k_{j i} t + Δ t$ adopts the harmonic average scheme. Thus,

.

T_{i, i j}^{t + Δ t} = T_{j, i j}^{t + Δ t}

(27)

Similar to Eq. 17, the transmissibilities $T_{i, i j} t + Δ t$ and $T_{j, i j} t + Δ t$ are equal, so they can be denoted together as $T_{i j} t + Δ t$ ⁠, indicating that the local mass conservation is satisfied, that is written as:

.

T_{i j}^{t + Δ t} = T_{i, i j}^{t + Δ t} = T_{j, i j}^{t + Δ t} = k_{i j}^{t + Δ t} V_{i} (e_{3, j}^{i} + e_{4, j}^{i}) = k_{j i}^{t + Δ t} V_{j} (e_{3, i}^{j} + e_{4, i}^{j})

(28)

Then Eq. 24 can be rewritten as:

.

\begin{array}{l} \sum_{j = 1}^{n_{i}} [\frac{k_{r a, i j}^{t + Δ t}}{B_{a, i j}^{t + Δ t} μ_{a, i j}^{t + Δ t}} T_{i j}^{t + Δ t} (p_{a, (i, j)}^{t + Δ t} - p_{a, i}^{t + Δ t})] + f l a g (i) \frac{k_{r a, i f}^{t + Δ t} T_{i f}^{t + Δ t}}{B_{a, i f}^{t + Δ t} μ_{a, i f}^{t + Δ t}} (p_{a, f}^{t + Δ t} - p_{a, i}^{t + Δ t}) \\ + Q_{a, w e l l} = \frac{V_{i}}{Δ t} [(\frac{ϕ_{i}^{t + Δ t} S_{a, i}^{t + Δ t}}{B_{a, i}^{t + Δ t}}) - (\frac{ϕ_{i}^{t} S_{a, i}^{t}}{B_{a, i}^{t}})] \end{array}

(29)

Similarly, for the flow between fracture nodes in the fracture node layer, in the case of the 2D reservoir model, the flow along fractures is a 1D flow on 1D fracture line segments. Therefore, for a fracture node f, suppose the matrix node corresponding to the fracture node f is the matrix node i. Then the discrete scheme similar to Eq. 29 of oil/water two-phase flow in the discrete fracture system (i.e., the fracture node layer) can be obtained as follows:

.

\begin{array}{l} \sum_{g = 1}^{n_{i}} [\frac{k_{r a, f g}^{t + Δ t}}{B_{a, f g}^{t + Δ t} μ_{o, f g}^{t + Δ t}} T_{f g} (p_{a, (f, g)}^{t + Δ t} - p_{a, f}^{t + Δ t})] + \frac{k_{r a, i f}^{t + Δ t}}{B_{a, i f}^{t + Δ t} μ_{a, i f}^{t + Δ t}} T_{i f}^{t + Δ t} (p_{a, i}^{t + Δ t} - p_{a, f}^{t + Δ t}) \\ + Q_{a, w e l l} = \frac{V_{f}}{Δ t} [ϕ_{f}^{t + Δ t} S_{a, f}^{t + Δ t} - ϕ_{f}^{t} S_{a, f}^{t}] \end{array}

(30)

The second terms on the left hand side of Eq. 30 are the mass transfers between the matrix node and the corresponding fracture node, wherein, $f l a g (i)$ is the function to judge whether the matrix node i has a corresponding fracture node, that is, to judge whether a fracture pass across the node. If so, $f l a g (i) = 1$ ⁠. If not, $f l a g (i) = 0$ ⁠. $T_{i f}$ is the transmissibility to measure the mass transfer between the matrix node i and the fracture node f. In the case of the 2D reservoir domain, the fracture nodes are connected by 1D line segments, and the transmissibility between fracture nodes in Eq. 30 can be simply calculated by:

.

T_{f g}^{t + Δ t} = k_{f g}^{t + Δ t} h_{f g} w_{f g} / l_{f g}

(31)

It should be pointed out that, if the reservoir calculation domain is 3D, the fracture nodes are distributed on the 2D surface, then it is required to use Eq. 28 on the 2D surface based on EFVM to obtain the transmissibility between fracture nodes.

However, it should be noted that the source or sink term on the fracture node layer is generally the perforation point of the fractured horizontal well. At this time, the well index formula is no longer the Eq. C-3 in Appendix C, but is modified as follows:

.

W I_{f} = \frac{2 π k_{f} w_{f}}{\ln (r_{e, f} / r_{w}) + s}, r_{e, f} = 0.14 \sqrt{2 V_{f} / w_{f}}

(32)

Definition of Connection Types between Nodes and Transmissibility Calculation

It can be seen that the discrete schemes of the flow governing equations at the matrix node and the fracture node can be written into the same form in Eq. 33, but $n_{i}$ in the equation is replaced by $N_{i}$ which denotes the number of nodes (including matrix nodes and fracture nodes) neighboring to node i, that is, the first term on the left hand side of Eq. 33 contains all mass transfer related to the control volume of node i, which is similar to EDFM based on mesh division and FVM. According to Eq. 33, by determining the connection (neighboring) relationship between nodes (including matrix nodes and fracture nodes) and then calculating the corresponding transmissibility according to the type of connection (neighboring) relationship, it can be seen that there are three types of connections in MFDFM: matrix node-matrix node, fracture node-fracture node, and matrix node-fracture node.

.

\sum_{j = 1}^{N_{i}} [\frac{λ_{a, i j}^{t + Δ t}}{B_{a, i j}^{t + Δ t}} T_{i j} (p_{a, (i, j)}^{t + Δ t} - p_{a, i}^{t + Δ t})] + Q_{a, w e l l} = \frac{V_{i}}{Δ t} [ϕ_{i}^{t + Δ t} S_{a, i}^{t + Δ t} - ϕ_{f}^{t} S_{a, i}^{t}]

(33)

For the first-type connection between the matrix node and the matrix node, the corresponding transmissibility is calculated by Eq. 28.
For the second-type connection between the fracture node and the matrix node, transmissibility is calculated by Eq. 31.
For the third-type connection, although EFVM defines and calculates the control volume of each node, the actual shape of the control volume of each node is related to the node distribution in the influence domain of node i, so it is difficult to accurately describe the actual shape of the node control domain. Because a node is nearly located in the center of the node control domain, as shown in Fig. 8, the node control domain is simplified to a square area centered at the node on the x-y plane, and let the fracture segment connected to the fracture node is located on the symmetry axis. At this time, the side length of the square is $\sqrt{V_{i} / h_{i}}$ ⁠, and the average distance from the point in the control domain to the fracture line segment is $\sqrt{V_{i} / h_{i}} / 4$ ⁠. Suppose that the fracture node f is neighboring to matrix node i, the transmissibility $T_{i f} t + Δ t$ in Eq. 29 and Eq. 30 is calculated by Eq. 34, which is similar to the transmissibility formula in EDFM (Li and Lee 2008):

.

T_{i f}^{t + Δ t} = 2 {({(\frac{k_{i}^{t + Δ t}}{\sqrt{V / h_{i}} / 4} A_{f})}^{- 1} + {(\frac{k_{f}^{t + Δ t}}{w_{f} / 2} A_{f})}^{- 1})}^{- 1}

(34)

Fig. 8

View large Download slide

Sketch of geometry approximation of node control volume for calculation of transmissibility of the third-type connection.

The fully implicit discrete schemes of flow equations at each real node can be obtained by combining Eqs. 28 , 31, 33, and 34. As pointed out earlier, EFVM is a meshless FVM based on meshless collocation. Compared with the commonly used FVM-based reservoir numerical simulator, EFVM can deal with more complex boundary conditions. Therefore, in the next section, we will introduce the treatment of various boundary conditions in the proposed MFDFM, point out the specific composition of global equations, and specify which nodes participate in the calculation of global equations and the corresponding equation type of each node.

Treatment of Boundary Conditions and Construction of Global Equations

As mentioned above, the EFVM on which MFDFM is based can handle more complex boundary conditions than the traditional FVM, and the EFVM handles the boundary conditions in the same way as GFDM. As seen from “Principles of GFDM,” GFDM uses the node coordinate information in the node influence domain in combination with Taylor expansion to obtain the generalized difference operator that minimizes the weighted truncation errors of Taylor expansion, that is, the weighted least square method commonly used to obtain the spatial derivatives of unknown function in various meshless methods. Therefore, the accuracy of the generalized finite difference approximation is related to the quality of the node distribution participating in the construction of the generalized finite difference expression of spatial derivatives at the considered node. Since there are no other nodes outside the boundary, if no virtual nodes are added, the node distribution uniformity in the boundary node influence domain will be low, that is, the position of the center of gravity of all nodes in the boundary node influence domain is much different from that of the considered boundary node. At this time, the accuracy of the generalized finite difference approximation in Eq. 10 will be reduced. Therefore, in GFDM, scholars (Milewski 2012; Rao 2022a) have recognized that a virtual node can be added outside the boundary where the boundary node locates to improve the node distribution quality, then improve the accuracy of generalized finite difference approximation at the boundary node to handle derivative boundary conditions. For the details about the specific data analysis of this issue, readers can refer to Rao (2022b) and Rao et al. (2022a).

Fig. 9 is a schematic representation of adding virtual nodes at derivative boundary nodes in Rao (2022a). It can be seen that when blue virtual nodes are not added in Fig. 9a, the black nodes in the influence domain of the red boundary node are on the same side of the tangent at the red boundary node. After adding a blue virtual node in the influence domain of the boundary node in Fig. 9b, the gravity center of the new node distribution is closer to the red boundary node, such that the node distribution quality and approximation accuracy of Eq. 10 are improved. In addition, the virtual node corresponding to the derivative boundary node pointed out in “Meshless Point Cloud Discretization of the Fracture-Reservoir Domain” can be directly used as the virtual node required to deal with the derivative boundary conditions, so the virtual node will also participate in the construction of global equations.

Fig. 9

View large Download slide

Sketch of adding the virtual node. After Rao (2022a). (a) The original case. (b) The case after adding the virtual node.

Taking the pressure function as an example, if a boundary node i satisfies the Dirichlet boundary conditions $p_{o} | Γ = f (x, y)$ ⁠, $S_{w} | Γ = g (x, y)$ ⁠, the discrete equations at node i are given as follows:

,

p_{o, i} = f (x_{i}, y_{i})

(35)

,

S_{w, i} = g (x_{i}, y_{i})

(36)

where $x_{i}$ and $y_{i}$ are the x coordinate and y coordinate of node i, respectively.

If a boundary node i satisfies the second type of boundary conditions $\frac{\partial p_{o}}{\partial n} |_{Γ} = f (x, y)$ ⁠, suppose that the added virtual node for handling derivative boundary conditions of node i is node j, in the global equations, the equations at node i are the discrete scheme of the flow governing equation in Eq. 33. The equations at node j are the discrete schemes of corresponding boundary conditions, these are as follows:

,

n_{x} \sum_{l = 1}^{n_{i}} m_{1 l}^{i} (p_{o, (i, l)} - p_{o, i}) + n_{y} \sum_{l = 1}^{n_{i}} m_{2 l}^{i} (p_{o, (i, l)} - p_{o, i}) = f (x_{i}, y_{i})

(37)

where $n_{x}$ and $n_{y}$ are the components of the external normal vector $n$ ⁠, $n_{i}$ is the number of nodes (including virtual nodes) in the influence domain of node i except node i itself. Since virtual node j is in the influence domain of node i, there is l satisfying $(i, l) = j$ ⁠.

For the third type of boundary conditions (i.e., mixed boundary conditions) $α \frac{\partial p_{o}}{\partial n} |_{Γ} + β p_{o} = f (x, y)$ that are also derivative boundary conditions, the treatment method is the same as that for the second type of boundary conditions. The discrete schemes of the boundary conditions are as follows:

.

α [n_{x} \sum_{l = 1}^{n_{i}} m_{1 l}^{i} (p_{o, (i, l)} - {p_{o}}_{, i}) + n_{y} \sum_{l = 1}^{n_{i}} m_{2 l}^{i} (p_{o},_{(i, l)} - p_{o, i})] + β p_{o, i} = f (x_{i}, y_{i})

(38)

Therefore, if there are $n_{m}$ matrix nodes and $n_{f}$ fracture nodes in the reservoir calculation domain, including $n_{1 B C}$ first-type (Dirichlet) boundary condition nodes and $n_{2 B C}$ derivative-type boundary condition nodes, there are also $n_{2 B C}$ of all virtual nodes participating in the construction of global equations, and the nodes participating in the construction of global equations are $n_{m}$ real matrix nodes, $n_{f}$ real matrix nodes and virtual nodes corresponding to derivative boundary condition nodes, respectively. Among them, the equations corresponding to $n_{m} - n_{1 B C} - n_{2 B C}$ inner matrix nodes, $n_{f}$ fracture nodes, and $n_{2 B C}$ derivative boundary condition nodes are the discrete schemes in Eq. 33 of flow governing equations, the equations corresponding to $n_{1 B C}$ first-type boundary condition nodes are Eq. 35 and Eq. 36, and the equations corresponding to $n_{2 B C}$ virtual nodes of derivative boundary condition nodes are Eq. 37 or Eq. 38.

Therefore, there are only two types of equations contained in the global equations, one is the discrete schemes of the flow governing equations, and the other is the discrete equation of boundary conditions. Therefore, to construct the global equations efficiently, the first step is to construct the discrete schemes of flow governing equations efficiently by using the connection relationship and transmissibility similar to that in the numerical simulator based on FVM and then fill the derivative boundary condition Eq. 37 or Eq. 38 at the virtual node corresponding to the derivative boundary condition node. Finally, the first-type boundary condition Eq. 35 and Eq. 36 are filled at the virtual node corresponding to the first-type boundary condition node. Of course, if the bottomhole pressure (BHP) or production rate needs to be calculated, the equation of well production rate and BHP needs to be added. For the formed global equations, the nonlinear solver based on Newton iteration and automatic differentiation (Zhou et al. 2011; Younis 2011; Voskov and Tchelepi 2011, 2012) is used to solve them. It can be seen that Eqs. 35 through 38 representing boundary conditions are generally linear. Therefore, compared with EDFM based on grid and FVM, MFDFM does not increase the nonlinearity of the equations. As a result, it can be expected that the computational efficiency of MFDFM is comparable with that of EDFM. The numerical examples in “Numerical Examples” will prove this assertation with the specific data that spent Newton iteration steps vs. production time.

Additionally, two notes need to be given here:

It should be noted that the above global equations can be solved by adding linear equations corresponding to boundary conditions on the basis of the existing nonlinear solver in the FVM-based reservoir numerical simulator, which reduces the difficulty of constructing an MFDFM-based general purpose simulator. In particular, when it is the most common closed boundary condition in reservoir simulation, the discrete equation of the boundary condition is not required, but only the connection table composed of the connection between the nodes and the corresponding transmissibility is constructed without considering the connection between the boundary nodes and the virtual nodes in the first item on the left of the discrete scheme in Eq. 33 of flow governing equations at the boundary nodes. Based on this connection table, the equations are then directly solved using an existing nonlinear solver without adding equations of the boundary condition. That is to say, construct a corresponding data file via NNCs and call a commercial numerical simulation software.
We thought that the fracture often would not touch the boundary of the reservoir calculation domain, so there is not much attention to consider the case that fracture nodes are on the reservoir boundary. However, it is possible to imagine that if a fracture node is on a boundary and the boundary condition is closed, then no additional treatment of the fracture node is required because the fracture node is not connected to a node on the outside of the boundary, so there is no mass transfer between the fracture node and the outside of the boundary such that the closed boundary condition is naturally satisfied. If the boundary is not closed, but is a first-type boundary condition, the pressure or saturation at the fracture node can be taken to the value required by the first-type boundary condition. These two cases are easy to handle; however, if the boundary condition is a nonzero Newman boundary condition or a complex third-type boundary condition, the boundary condition needs only be assigned to the nodes of the matrix system with real edges on the boundary and not to the fracture system in one dimension lower than matrix system which has only a much small contact with the boundary. In summary, we prefer not to assign the boundary condition of the reservoir computational domain to the fracture system, which is one dimension lower than the reservoir computational domain, but only to the matrix system. From an approximation perspective, when there is a fracture node on the boundary, ignoring the effect of boundary conditions on the fracture node is actually equivalent to slightly moving the fracture node to the interior of the reservoir domain, which is also acceptable. In FVM-based EDFM, the fracture system is also considered in terms of its mass transfer with the matrix system and rarely in terms of its direct transfer with the boundary.

Numerical Examples

This section will give three comprehensive numerical examples to test the computational performances of MFDFM in regular and irregular collocations, homogeneous and heterogeneous reservoirs, regular boundary geometry and irregular boundary geometry, closed boundary condition, and Dirichlet boundary condition, thereby prove that the proposed MFDFM can obtain the nearly same calculation results as the traditional methods, but has the advantages of flexibly discretizing the computational domain and dealing with various boundary conditions. Note that since the proposed MFDFM is based on a meshless point cloud, the calculated pressure and saturation values are the values at the nodes, and the node distribution may be irregular. Therefore, the proposed MFDFM uses the nodal values to plot the field diagram through interpolation. The codes of the proposed MFDFM are currently written based on MATLAB, and the scatter function in MATLAB is used to plot the field diagrams in this paper, so the pressure and saturation values in field diagrams appear to be continuously distributed. To make EDFM results have the same resolution as the MFDFM results for effective comparisons, we interpolate the grid pressure and saturation values in EDFM according to the grid center coordinates as well.

MFDFM Performances for Nonintersecting Fractures With a Matching Point Cloud

As shown in Fig. 10a, this section designs a rectangular reservoir calculation area of 500×200 m, with a thickness of 10 m. In the model, the blue circles denote four water injection wells at the four corners of the reservoir, the black line segments denote the five hydraulic fractures, and the blue line denotes the wellbore of the fractured well. The fractured well is controlled by a constant BHP of 10 MPa, and four water injection wells are controlled by a fixed BHP of 25 MPa. The relevant physical properties and relative permeability data are shown in Tables 1 and 2, respectively.

Fig. 10

View large Download slide

Sketch of discretization of the reservoir model in MFDFM. (a) The reservoir model. (b) The matching point cloud. (c) Matrix-node layer in the matching point cloud generated by using the mesh vertices in Delaunay triangulation. (d) Fracture-node layer in the matching point cloud generated by using the mesh vertices in Delaunay triangulation. (e) The calculated node control volume profile of the matrix nodes.

Table 1

Physical properties of the reservoir model.

Properties	Values	Properties	Values
Matrix porosity	0.1	Matrix permeability	3 md
Fracture porosity	0.3	Fracture aperture	0.01 m
Oil compressibility	3×10⁻³ MPa⁻¹	Water compressibility	5×10⁻⁴ MPa⁻¹
Rock compressibility	1×10⁻⁴ MPa⁻¹	Oil viscosity	2 mPa·s
Water viscosity	0.6 mPa·s	Initial reservoir pressure	15 MPa
Initial water saturation	0.15	Oil volume factor	1.0
Water volume factor	1.0	Reservoir thickness	10 m
Well radius	0.1 m

Properties	Values	Properties	Values
Matrix porosity	0.1	Matrix permeability	3 md
Fracture porosity	0.3	Fracture aperture	0.01 m
Oil compressibility	3×10⁻³ MPa⁻¹	Water compressibility	5×10⁻⁴ MPa⁻¹
Rock compressibility	1×10⁻⁴ MPa⁻¹	Oil viscosity	2 mPa·s
Water viscosity	0.6 mPa·s	Initial reservoir pressure	15 MPa
Initial water saturation	0.15	Oil volume factor	1.0
Water volume factor	1.0	Reservoir thickness	10 m
Well radius	0.1 m

Table 2

Relative permeability data.

Water Saturation	Oil Relative Permeability	Water Relative Permeability	Water Saturation	Oil Relative Permeability	Water Relative Permeability
0.15	0	0.8	0.6276	0.207	0.063
0.45	0.04	0.4	0.6489	0.247	0.046
0.5212	0.08	0.227	0.6701	0.292	0.031
0.5425	0.103	0.184	0.6914	0.335	0.021
0.5638	0.125	0.156	0.7127	0.378	0.011
0.5851	0.145	0.125	0.76	0.475	0.005
0.6064	0.174	0.094	0.85	0.552	0

Water Saturation	Oil Relative Permeability	Water Relative Permeability	Water Saturation	Oil Relative Permeability	Water Relative Permeability
0.15	0	0.8	0.6276	0.207	0.063
0.45	0.04	0.4	0.6489	0.247	0.046
0.5212	0.08	0.227	0.6701	0.292	0.031
0.5425	0.103	0.184	0.6914	0.335	0.021
0.5638	0.125	0.156	0.7127	0.378	0.011
0.5851	0.145	0.125	0.76	0.475	0.005
0.6064	0.174	0.094	0.85	0.552	0

A matching point cloud of the reservoir model is shown in Fig. 10b, in which the blue dotted line denotes the fractured well with the matching well treatment. Figs. 10c and 10d show the matrix nodes and fracture nodes in the point cloud. Because the point cloud is matching, each fracture node in Fig. 10d corresponds to a matrix node in the same position in Fig. 10c. Fig. 10e illustrates the calculated node volumes by Eq. 13 and Eq. 16. It can be seen that the control volume of nodes in a region with a high node density is small, which is physical.

Case 1:Fig. 11 compares reservoir pressure [color code range (10 MPa, 23.5 MPa)] and oil saturation [color code range (0.15, 0.85)] profiles at the 800th day calculated by MFDFM and EDFM, and the radius of node influence domain in MFDFM is set as 20 m. Fig. 12 compares the oil production rates and water cut calculated by the two methods. It can be seen that the pressure and saturation distributions calculated by MFDFM and EDFM are very close, especially the pressure distribution. In addition, the oil production rates calculated by the two methods almost coincide, and the water-cut curves are slightly different, but also very similar. The comparison results show that MFDFM with the matching point cloud can obtain very similar calculation results with widely used EDFM, verifying the good performance of the proposed MFDFM in the homogeneous reservoir model with nonintersecting fractures.
Case 2: This case changes the matrix permeability profile to meet Eq. 39, which is a heterogeneous reservoir model.

Fig. 11

View large Download slide

Comparisons of pressure and oil saturation profiles at the 800th day calculated by EDFM and MFDFM. (a) Pressure, 800 days, MFDFM. (b) Pressure, 800 days, EDFM. (c) Oil saturation, 800 days, MFDFM. (d) Oil saturation, 800 days, EDFM.

Fig. 12

View large Download slide

Comparisons of oil production rates and water cut vs. production time. (a) Oil production rate. (b) Water cut.

.

k (x, y) = 3 \times e x p (2 [{(x / 500)}^{2} + {(y / 200)}^{2}]) m d

(39)

Fig. 13 compares reservoir pressure [color code range (10 MPa, 24 mPa)] and oil saturation [color code range (0.15, 0.85)] profiles at the 400th and 800th days calculated by MFDFM and EDFM in this heterogeneous case. Fig. 14 compares the oil production rates and water cut calculated by the two methods. We can observe that the pressure and saturation profiles calculated by MFDFM and EDFM are close. The oil production rates and the water cut vs. production time calculated by the two methods are slightly different but also very similar. Fig. 15 shows the cumulative Newton iteration steps vs. production time, from which it can be seen that EDFM requires more Newton iteration steps than MFDFM when dealing with heterogeneous models. This suggests that the computational cost of MFDFM may be less than or equal to that of EDFM. The comparisons verify the good performance of the MFDFM in the heterogeneous reservoir model.

Fig. 13

View large Download slide

Comparisons of pressure and oil saturation profiles at the 400th and the 800th days calculated by EDFM and MFDFM in the heterogeneous case. (a) Pressure, 400 days, MFDFM. (b) Pressure, 400 days, EDFM. (c) Pressure, 800 days, MFDFM. (d) Pressure, 800 days, EDFM. (e) Oil saturation, 400 days, MFDFM. (f) Oil saturation, 400 days, EDFM. (g) Oil saturation, 800 days, MFDFM. (h) Oil saturation, 800 days, EDFM.

Fig. 14

View large Download slide

Comparisons of oil saturation profiles calculated by EDFM and MFDFM in the heterogeneous case. (a) Oil production rate. (b) Water cut.

Fig. 15

View large Download slide

Comparisons of Newton iteration steps vs. production time.

MFDFM Performances for a Complex Orthogonal Fracture Network With Matching and Nonmatching Point Clouds

As shown in Fig. 16a, this example is a rectangular reservoir model with an orthogonal complex fracture network consisting of 124 explicitly characterized fractures parallel only to the x- or y-axes, with the blue line representing a fractured horizontal well and the intersection of this blue line with fractures being the perforating points of the fractured well. Except for that matrix permeability and fracture permeability become 1 and 50,000 md, respectively, other physical parameters are the same as those in Table 1. Water injection wells are controlled by a constant injection rate 15 m³/d, and production wells are controlled by a constant BHP 10 MPa. This example uses the nonmatching point cloud shown in Figs. 16b through 16d and the matching point cloud shown in Figs. 16f through 16h, respectively. Let the radius of the node influence domain of the nonmatching point cloud be 15 m and the distance between the virtual node and the corresponding boundary node be 10 m. The radius of the node influence domain of the matching point cloud is set to 25 m, and the distance between the virtual node and the corresponding boundary node is 15 m. Figs. 16e and 16i show the distribution of the node control volumes (not multiplied by the reservoir thickness) of these two point clouds calculated by Eqs. 13 and 16, respectively. It can be seen that in the nonmatching point cloud, since the matrix nodes are Cartesian, the internal node control volumes in the point cloud are all the same 100 m², while the control volumes of the four corner nodes of the domain are one-quarter of the control volumes of the inner nodes, and the control volumes of the boundary nodes other than these four corner nodes are one-half of the control volumes of the inner nodes, which are physical. In the matching point cloud, the control volumes of the matrix nodes corresponding to the denser fracture nodes are also generally smaller, and the areas of the overall node control volume distribution where the node control volume is small almost reflect the geometry of the fracture network, implying that the node control volumes calculated based on the matching point cloud are also physical.

Fig. 16

$Sketch of the reservoir domain, the nonmatching point cloud, and the matching point cloud. (a) Reservoir domain, in which two blue circles denote the injection wells, and the blue line denotes a fractured horizontal well. (b) All nodes in the non-matching point cloud. (c) Matrix nodes in the non-matching point cloud. (d) Fracture nodes in the non-matching point cloud. (e) Node control volume profile of the matrix nodes in the non-matching point cloud. (f) All nodes in the matching point cloud. (g) Matrix nodes in the matching point cloud. (h) Fracture nodes in the matching point cloud. (i) Node control volume profile of the non-matching point cloud.$

View large Download slide

Sketch of the reservoir domain, the nonmatching point cloud, and the matching point cloud. (a) Reservoir domain, in which two blue circles denote the injection wells, and the blue line denotes a fractured horizontal well. (b) All nodes in the non-matching point cloud. (c) Matrix nodes in the non-matching point cloud. (d) Fracture nodes in the non-matching point cloud. (e) Node control volume profile of the matrix nodes in the non-matching point cloud. (f) All nodes in the matching point cloud. (g) Matrix nodes in the matching point cloud. (h) Fracture nodes in the matching point cloud. (i) Node control volume profile of the non-matching point cloud.

Figs. 17 through 19a and 19b compare the pressure and oil saturation distributions at 1,000 and 2,500 days for EDFM, MFDFM with the nonmatching point cloud, and MFDFM with the matching point cloud, respectively, as well as oil production rates and water cut. It can be seen that MFDFM can achieve results very similar to EDFM, whether based on unmatched or matched point clouds, and the good performance of the matching point cloud strongly suggests that MFDFM with matching point clouds can avoid gridding challenges in narrow local areas compared with DFM based on matching grids.

Fig. 17

View large Download slide

Comparisons of the pressure and oil saturation profiles calculated by EDFM and MFDFM. (a) EDFM, pressure, 1,000 days. (b) EDFM, oil saturation, 1,000 days. (c) MFDFM with the non-matching point cloud, pressure, 1,000 days. (d) MFDFM with the non-matching point cloud, oil saturation, 1,000 days. (e) MFDFM with the matching point cloud, pressure, 1,000 days. (f) MFDFM with the matching point cloud, oil saturation, 1,000 days.

Fig. 18

View large Download slide

Comparisons of the pressure and oil saturation profiles calculated by EDFM and MFDFM. (a) EDFM, pressure, 2,500 days. (b) EDFM, oil saturation, 2,500 days. (c) MFDFM with the non-matching point cloud, pressure, 2,500 days. (d) MFDFM with the non-matching point cloud, oil saturation, 2,500 days. (e) MFDFM with the matching point cloud, pressure, 2,500 days. (f) MFDFM with the matching point cloud, oil saturation, 2,500 days.

Fig. 19

View large Download slide

Comparisons of production data and required Newton iteration steps. (a) water cut. (b) cumulative Newton iteration steps versus production time.

Fig. 19c compares the cumulative Newton iteration steps required for EDFM, MFDFM with the nonmatching point cloud, and MFDFM with the matching point cloud, and it can be seen that the MFDFM with the nonmatching point cloud requires fewer Newton iteration steps than EDFM over time. Although the MFDFM with the matching point cloud requires significantly more Newton iteration steps, this is because many of the matrix nodes in the matching point cloud have smaller control volumes, which increase the difficulty of convergence in the nonlinear solution, which is natural and acceptable. This comparison result argues for the assertion that the proposed MFDFM as described in “Treatment of Boundary Conditions and Construction of Global Equations” does not increase the degree of nonlinearity in the global set of equations compared with grid-based methods, and that MFDFM with a nonmatching point cloud may be computationally more efficient than EDFM.

MFDFM Performances for a Nonorthogonal Complex Fracture Network and a Mixed Boundary Condition

As shown in Fig. 20a, this example is a rectangular reservoir with a nonorthogonal complex fracture network, with two injection wells marked with blue circles on the left side of the reservoir and two production wells marked with blue circles on the right side of the reservoir. Except for that matrix permeability and fracture permeability become 5 and 50,000 md, respectively, other physical parameters are the same as those in Table 1. Water injection wells are controlled by a constant injection rate 30 m³/d, and production wells are controlled by a constant BHP 5 MPa. Figs. 20b through 20d show a nonmatching point cloud for this computational domain. This example was constructed with two main objectives in mind, one is to check the computational performance of MFDFM for nonorthogonal complex fracture networks, and the other is to check the computational performance of MFDFM when dealing with mixed boundary conditions. Therefore, this example first calculates the pressure and oil saturation distributions by EDFM at 500 and 2,500 days for the closed boundary condition, which are shown in Figs. 21a and 21b,andFigs. 22a and 22b, respectively and then calculates the pressure and oil saturation distributions by MFDFM at 500 and 2,500 days for the closed boundary condition, which are shown in Figs. 21c and 21d and Figs. 22c and 22d. It can be seen that MFDFM results are very close to the EDFM results, verifying the good performance of MFDFM to handle a nonorthogonal fracture network. Fig. 23 shows the cumulative Newton iteration steps vs. production time, which suggests MFDFM spends significantly fewer Newton iteration steps than EDFM.

Fig. 20

View large Download slide

Sketches of the reservoir domain and the point cloud. (a) The reservoir domain. (b) The point cloud. (c) Matrix nodes in the point cloud. (d) Fracture nodes in the point cloud.

Fig. 21

View large Download slide

Comparisons of calculated pressure and oil saturation profiles at the 500th day. (a) EDFM, closed boundary condition, pressure. (b) EDFM, closed boundary condition, oil saturation. (c) MFDFM, closed boundary condition, pressure. (d) MFDFM, closed boundary condition, oil saturation. (e) MFDFM, mixed boundary condition, pressure. (f) MFDFM, mixed boundary condition, oil saturation. (g) MFDFM, mixed boundary condition, pressure. (h) MFDFM, mixed boundary condition, oil saturation.

Fig. 22

View large Download slide

Comparisons of calculated pressure and oil saturation profiles at the 2,500th day. (a) EDFM, closed boundary condition, pressure. (b) EDFM, closed boundary condition, oil saturation. (c) MFDFM, closed boundary condition, pressure. (d) MFDFM, closed boundary condition, oil saturation. (e) MFDFM, mixed boundary condition, pressure. (f) MFDFM, mixed boundary condition, oil saturation. (g) MFDFM, mixed boundary condition, pressure. (h) MFDFM, mixed boundary condition, oil saturation.

Fig. 23

View large Download slide

Comparisons of the cumulative Newton steps vs. production time.

Considering the difficulty in obtaining the reference solution of EDFM for the general third-type boundary condition, the third type of boundary condition on pressure shown in Eq. 40 is first constructed in this section. It can be seen that the boundary condition in Eq. 40 is very close to the closed boundary condition because the value of $β$ in this case is too small. By comparing the MFDFM results with the EDFM results shown in Figs. 21e and 21f and Figs. 22e and 22f, which are based on the boundary condition in Eq. 40, we can see that the results are much similar with EDFM results, which verify the accuracy of the MFDFM for the third type of boundary condition. Finally, we construct the third type of boundary condition in Eq. 41 by increasing the value of $β$ ⁠. Since the pressure on the left side of the reservoir is gradually increased by the injection wells, the third type of boundary condition in Eq. 41 theoretically causes the pressure increase at the left side of the reservoir boundary to be lost to the outside of the computational domain, thus making the pressure difference drive of the injection water from the left side to the right side weaker and slowing down the advance of the injection water. The calculated pressure profiles in Figs. 21g and 22g illustrate that the overall reservoir pressure level is significantly lower than the EDFM calculations. The oil saturation profiles in Figs. 21h and 22h illustrate a slowing of the rate of advance of the injection water toward the producing well, which is particularly evident in Fig. 22h.

,

\frac{\partial p_{o}}{\partial n} | Γ + β (p_{o} - 15) = 0, β = 10^{- 6}

(40)

.

\frac{\partial p_{o}}{\partial n} | Γ + β (p_{o} - 15) = 0, β = 10^{- 3}

(41)

MFDFM Performances With Complex Reservoir-Boundary Geometry

As shown in Fig. 24a, in a fractured reservoir with a complex heart-shaped boundary geometry, there are two fractured horizontal wells (marked by blue lines) with a complex fracture network (marked by intersecting black lines). The nonmatching point cloud for the reservoir domain is used, in which Figs. 24b through 24d show all nodes, matrix nodes, and fracture nodes in the point cloud, respectively. To obtain the mesh-based reference solution, we construct a larger rectangular reservoir model in Fig. 24e capable of containing the reservoir model discretized by a Cartesian grid by assigning the permeability of the grid outside the cardioid reservoir boundary to 0 in Fig. 24f. The difference between MLDFM and mesh-based EDFM in dealing with such a heart-shaped complex reservoir boundary shows the advantage of meshless MLDFM over the mesh-based method in discretizing computational domains with complex geometry.

Fig. 24

View large Download slide

Sketch of discretization of the reservoir model in MFDFM. (a) Reservoir domain. (b) All nodes. (c) Matrix nodes. (d) Fracture nodes. (e) The EDFM grid to handle this model. (f) The permeability distribution of the EDFM grid. (g) The radius of the node influence domain is 40 m. (g) The radius of the node influence domain is 50 m.

The physical parameters of the reservoir model in this example are basically the same as those in Table 1, except that the permeability is changed to 1 md, the fracturing horizontal well is controlled by a constant BHP of 10 MPa, and the initial water saturation is the irreducible water saturation 0.15. Therefore, it is essentially a single-phase flow problem. Because the relative permeability curve shown in Table 2 is still used, the effective permeability is the oil phase permeability corresponding to irreducible water saturation of 0.8 md.

Node influence domain radii 30, 40, and 50 m, respectively, are used for calculation, and the corresponding distances $Δ l$ in Eq. 21 of adding the virtual nodes are selected as 20, 24, and 25 m, respectively. Figs. 24g and 24h show the calculated node control volumes. Noted that as described in “The MFDFM Based on EFVM,” the radius of the node influence domain radius of MFDFM should not be too large to ensure good accuracy for convection-dominated quantity (e.g., saturation, concentration, and temperature) and calculation efficiency, thus these three radii are selected. Figs. 25 and 26 show the calculated pressure distribution at 10 and 50 days by MFDFM and EDFM, respectively. Fig. 27compares the oil production rate vs. production time under different radii of the node influence domain. It is found that the MLDFM calculated results under different radii of the node influence domain are all very similar to those of EDFM, which verifies that MFDFM can effectively handle complex-reservoir-boundary geometry and illustrates that MFDFM can obtain relatively robust simulation results within a certain range of the influence domain radius, which indicates the possibility of constructing a robust MLDFM-based numerical simulator.

Fig. 25

View large Download slide

Comparisons of calculated pressure profiles at the 10th day by EDFM and MFDFM. (a) EDFM. (b) MDFDM, 30 m. (c) MFDFM, 40 m. (d) MFDFM, 50 m.

Fig. 26

View large Download slide

Comparisons of calculated pressure profiles at the 50th day by EDFM and MFDFM. (a) EDFM. (b) MFDFM, 30 m. (c) MFDFM, 40 m. (d) MFDFM, 50 m.

Fig. 27

View large Download slide

Comparisons of oil production rates vs. production time from EDFM and MFDFM. (a) Well 1. (b) Well 2.

Conclusions

Based on the newly proposed EFVM, this paper constructs a meshless numerical simulation method for fractured reservoirs. In general, this method has the following three characteristics:

MFDFM explicitly characterizes discrete fractures and uses the meshless point cloud to flexibly discretize the fractured reservoirs with complex geometry. Two types of point cloud discretization for the fractured reservoir domain including the matching type and the nonmatching type are developed, which illustrate that MFDFM avoids the gridding challenge of fracture networks and complex reservoir boundaries.
Based on the extended finite volume discrete scheme of flow governing equations, MFDFM constructs various types of internode connections and corresponding transmissibility formulas, which are then combined with various types of GFDM-based discretization equations of boundary conditions to construct the global equations.
The existing nonlinear solver in the current FVM-based reservoir numerical simulator can be used to solve the global equations, especially for closed boundary conditions, which reduces the difficulty of forming a general purpose MFDFM-based numerical simulator for fractured reservoirs.
Four numerical examples are implemented to verify that MFDFM achieves good computational performance for both homogeneous and heterogeneous reservoir models, disjoint fractures, orthogonal and nonorthogonal complex fracture networks, complex boundary geometry, and complex boundary conditions, providing the first meshless framework for numerical modeling of fractured reservoirs.

Future Works

We believe that six valuable works can be carried out in the future:

When the same radius of influence domain is used for each node to determine a connectable point cloud, the presence of some relatively dense local point clouds may lead to errors in the control volumes of the nodes calculated using Eqs. 13 and 16, for example, some node control volumes may be negative because in this case, the poor distribution quality of the local point clouds of some nodes leads to the calculated $(m 3 j i + m 4 j i)$ or $(m 3 i j + m 4 i j)$ in Eq. 13 being negative. Therefore, an important work that would help to improve the robustness of the proposed MFDFM is to develop a method for constructing a connectable point cloud of the reservoir calculation domain that guarantees the quality of the local point clouds, rather than simply using the same radius of influence domain at each node to achieve a connectable point cloud.
Expanding MFDFM to 3D to investigate its computational performances of 3D actual reservoir models may be an important future work.
This paper primarily demonstrates MFDFM’s computational performance in the fundamental two-phase flow model. Studying its application in compositional models is also a valuable future work.
By coupling efficient meshless methods to handle geomechanics in MFDFM, a new method can be formed to solve the hot topic of the coupled multiphase flow and geomechanics problems in fractured reservoirs.
In current MFDFM, a fracture node will only be connected to one matrix node, which is the same as in EDFM where a fracture cell is only connected to the matrix cell that contains it. As a result, the current MFDFM is likely to have the same limitations as EDFM (i.e., the limitations that pEDFM is trying to solve). As a result, it may be a valuable future work to study how to solve the limitations in MFDFM or build more complex connections between matrix nodes and fracture nodes.
Using the flexibility of the point cloud in the adaptive analysis, dynamic rather than static node refinement methods can be developed to improve the calculation accuracy at the early time production stage of fractured wells and avoid the problem that the static unchanged refined nodes hardly contribute to the improvement of calculation accuracy at the middle and late stage but increase the calculation cost.

Author Contributions

Xiang Rao: Conceptualization; Methodology; Software; Validation; Writing-original draft; Writing-review & editing; Formal analysis; Funding acquisition; Investigation, Visualization.

Hui Zhao: Visualization.

Yina Liu: Visualization.

Acknowledgments

Dr. Rao thanks the supports from the National Natural Science Foundation of China (No. 52104017), the Open Fund of Cooperative Innovation Center of Unconventional Oil and Gas (Ministry of Education & Hubei Province) (No. UOG2022-14), and the Open Fund of Hubei Key Laboratory of Drilling and Production Engineering for Oil and Gas (Yangtze University) (Grant No. YQZC202201).

Appendix A: An Algorithm for Generating a Nonmatching Point Cloud

Step 1: Calculating the Intersections between Fractures

The intersections between fractures are the hub of mass transfer between fractures, so the first step is to find out the intersections between fractures. Suppose the two endpoints of a fracture be A and B with coordinates $(a_{1}, b_{1})$ and $(a_{2}, b_{2})$ ⁠, respectively, then the parametric equation of the fracture line $\bar{A B}$ is:

(x, y) = (a_{1}, b_{1}) + ξ_{1} (a_{2} - a_{1}, b_{2} - b_{1})

(A-1)

where $(x, y)$ is the coordinate of the point on $\bar{A B}$ ⁠.

Suppose another fracture has two endpoints C and D with coordinates $(c_{1}, d_{1})$ and $(c_{2}, d_{2})$ ⁠, then the parametric equation of the fracture line $\bar{C D}$ is:

(x, y) = (c_{1}, d_{1}) + ξ_{2} (c_{2} - c_{1}, d_{2} - d_{1})

(A-2)

Combining Eqs. 1 and 2, we can obtain the system of linear equations in Eq. C-8 with respect to $ξ_{1}$ and $ξ_{2}$ ⁠.

.

[\begin{array}{cc} a_{2} - a_{1} & c_{1} - c_{2} \\ b_{2} - b_{1} & d_{1} - d_{2} \end{array}] [\begin{matrix} ξ_{1} \\ ξ_{2} \end{matrix}] = [\begin{matrix} c_{1} - a_{1} \\ d_{1} - b_{1} \end{matrix}]

(A-3)

Then, first judge whether the coefficient matrix of Eq. A-3 is singular, and if it is, then $\bar{A B}$ and $\bar{C D}$ are parallel. At this point, then there are no intersections to be computed. If nonsingular, by solving Eq. C-8 to obtain $ξ_{1}$ and $ξ_{2}$ ⁠. It is important to note that the intersection of the two fractures needs to be on both fractures, that is on the line segments $\bar{A B}$ and $\bar{C D}$ ⁠, so it is required that,

.

0 \leq ξ_{1} \leq 1, 0 \leq ξ_{2} \leq 1

(A-4)

If the calculated $ξ_{1}$ and $ξ_{2}$ satisfy the requirement in Eq. A-4, then $ξ_{1}$ and $ξ_{2}$ represent the values of the parameters in Eqs. A-1 and A-2 corresponding to the intersection of the two fractures, respectively.

Step 2: Generation of Fracture Nodes

After obtaining the intersections between the fractures in the first step, each fracture will be divided into several parts by the resulting fracture intersections, which can be characterized by a vector of parameters in each fracture parameter equation. For example, for a particular fracture f, the parameter corresponding to its two endpoints takes the value 0 or 1, respectively, and suppose the parameters corresponding to the fracture intersections of this fracture with other fractures calculated by Eq. A-3 as $ξ_{f} i$ ⁠, $i = 1, \dots, η f$ ⁠, where $η f$ is the number of intersections. Of course, it is possible that the values are equal to 0 or 1, or there will exist $ξ_{f} i = ξ_{f} j$ when $i \neq j$ ⁠, so it is necessary to eliminate duplicate values in the parameter vector. After obtaining the parameter vector for each fracture, suppose the parameter vector for a fracture f be $ξ_{f}$ ⁠, and the two adjacent parameter values in this parameter vector are $ξ_{f} i$ and $ξ_{f} i + 1$ ⁠, and the total length of this fracture is $l f$ ⁠, then the number of the divided segments of the fracture part between the points corresponding to the two adjacent parameter values is calculated as:

,

n_{f}^{i} = f l o o r [\frac{l_{f} (ξ_{f}^{i + 1} - ξ_{f}^{i})}{l_{c h a r}}] + 1

(A-5)

where $f l o o r (x)$ is the function refers to the rounding function, that is, when $n \leq x < n + 1$ and n is an integer, then $f l o o r (x) = n$ ⁠.

Then the true length of each segment in the fracture part is:

.

l_{r e a l} = \frac{l_{f} (ξ_{f}^{i + 1} - ξ_{f}^{i})}{n_{f}^{i}}

(A-6)

Step 3: Generation of Matrix Nodes

As mentioned before, the nonmatching point cloud directly ignores the fracture network morphology to generate the matrix nodes of the computational domain. For a rectangular boundary computational domain, as shown in Fig. 2a, the generation of matrix nodes can be done using a Cartesian point cloud. For a nonrectangular computational domain with irregular boundaries, a coordinate line parallel to the x-direction can first be constructed at an equal spacing in the y-direction, and then the nodes can be configured at an equal spacing in the x-direction according to the intersection of the coordinate line with the boundary. The spacing between nodes can be controlled by the characteristic length in the x-direction (denoted as $l_{x, c h a r}$ ⁠) and the characteristic length in the y-direction (denoted as $l_{y, c h a r}$ ⁠). In general, $l_{x, c h a r}$ can be equal to $l_{y, c h a r}$ ⁠, so they are togetherly denoted as $l_{c h a r}$ ⁠. As can be seen, this generation method of matrix nodes is indeed simpler than grid division in the computational domain. Of course, since nonmatching point clouds do not need to match the fracture network geometry, the various mesh division techniques can also be used directly to generate more uniformly distributed matrix nodes for complex boundary computational domains, for example, in Fig. 8, the mesh nodes of the Delaunay triangulation of the cardioid computational domain are used as matrix nodes in the point cloud of the domain.

Step 4: Construction of the NNCs between Matrix Nodes and Fracture Nodes

After obtaining fracture nodes and matrix nodes, this algorithm uses a simple method to establish the NNCs between fracture nodes and matrix nodes. In the meshless EFVM, although the node control volume is calculated, it is generally difficult to portray the specific shape of the node control domain. By using simple distance information between nodes, like that in Fig. 9, the node control domain can be assumed to be a square region. Suppose a fracture node f with coordinates $(x_{f}, y_{f})$ and matrix node i with coordinates $(x_{i}, y_{i})$ and calculated node control volume V_i, then the serial number (denoted by $ε_{f}$ ⁠) of the matrix node which has a NNC with f, is determined by

,

ε_{f} = \underset{i}{m i n} (| \frac{x_{f} - x_{i}}{l_{i}} | + | \frac{y_{f} - y_{i}}{l_{i}} |)

(A-7)

where $l_{i} = \sqrt{V_{i}}$ is the side length of the assumed square area.

It should be noted that the method to determine which matrix node has a NNC with f by finding which i to minimum Eq. A-7, reflecting the fact that a fracture node will generally have the NNC with only a matrix node in the proposed nonmatching point cloud generation algorithm, similar to EDFM in which a fracture cell only connected to one matrix cell which containing the fracture cell. In pEDFM, which addresses many of the limitations of EDFM, a fracture cell will have NNCs with multiple matrix cells, and this is a direction that this nonmatching point cloud generation algorithm could be further improved in the future.

Appendix B: An Algorithm for Generating a Matching Point Cloud

Step 1: Calculating the Intersections between Fractures

This step is the same as Step 1 in the nonmatching point cloud generation algorithm shown in Appendix A.

Step 2: Generation of Fracture Nodes

The fracture nodes are also denoted as the Type-I matrix nodes in the matching point cloud. This step is the same as Step 2 in the nonmatching point cloud generation algorithm shown in Appendix A.

Step 3: Generation of the Type-II Matrix Nodes

This step is the same as Step 3 in the nonmatching point cloud generation algorithm shown in Appendix A.

Step 4: Integration of the First Type of Matrix Nodes and the Second Type of Matrix Nodes

As mentioned before, since the matrix nodes in the matching point cloud consist of two parts, and the transmissibility between two matrix nodes is obtained based on EFVM, locally too close nodes may lead to poor quality node distribution in the influence domain of some nodes, resulting in reduced accuracy of the generalized difference approximation of spatial derivatives in the influence domain of the nodes. Therefore, this situation needs to be avoided as much as possible. Since the two types of matrix nodes are generated independently of each other (which is one of the main reasons why the generation of matching point clouds for complex fracture networks is easier and more flexible than matching mesh division), some nodes in the first-type matrix nodes may be close to or coincide with some nodes in the second-type matrix nodes.

For this reason, a threshold distance parameter (denoted $d i s_{t v}$ ⁠) is given in the algorithm, and suppose a node in the first-type matrix nodes be i with coordinates $(x_{i}, y_{i})$ ⁠, and a node in the second-type matrix nodes be j with coordinates $(x_{j}, y_{j})$ ⁠. If

,

\sqrt{{(x_{j} - x_{i})}^{2} + {(y_{j} - y_{i})}^{2}} \leq d i s_{t v}

(B-1)

then node j is removed or understood as node j in the second-type matrix nodes is replaced by node i in the first-type matrix nodes.

Appendix C: Derivation of the EFVM Discrete Scheme

Assuming that the pressure and water saturation of each node at time t have been known, integrate the node control volume of node i on both sides of Eq. 2, and take the implicit scheme of each physical quantity, that is, the values at the time $t + Δ t$ ⁠. Using GFDM based Eq. 10 and the work of EFVM by Rao (2022b), the integration of the first term of the left side of Eq. 2 is obtained as:

.

\int_{V_{i}} \nabla \cdot (\frac{k k_{r a}}{B_{a} μ_{a}} \nabla p_{a}) d Ω = \int_{V_{i}} \sum_{j = 1}^{n_{i}} [\frac{k_{i j}^{t + Δ t} k_{r a, i j}^{t + Δ t}}{B_{a, i j}^{t + Δ t} μ_{a, i j}^{t + Δ t}} (e_{3, j}^{i} + e_{4, j}^{i}) (p_{a, (i, j)}^{t + Δ t} - p_{a, i}^{t + Δ t})] d Ω

(C-1)

The SPU scheme in Eq. C-2 is adopted for phase permeability, the harmonic average scheme in Eq. C-3 is adopted for permeability, and the arithmetic average scheme in Eq. C-4 is adopted for volume coefficient and viscosity, which are expressed as follows, respectively.

,

\begin{array}{l} k_{r a, i j}^{t + Δ t} = {\begin{matrix} k_{r a, i} (S_{w, i}^{t + Δ t}) i f p_{a, (i, j)}^{t + Δ t} \geq p_{a, i}^{t + Δ t} \\ k_{r a, i} (S_{w, i}^{t + Δ t}) i f p_{a, (i, j)}^{t + Δ t} < p_{a, i}^{t + Δ t} \end{matrix} \end{array}

(C-2)

.

{k_{i j}}^{t + Δ t} = \frac{2}{1 / {k_{i}}^{t + Δ t} + 1 / {k_{j}}^{t + Δ t}}, {B_{a, i j}}^{t + Δ t} = \frac{{B_{a, i}}^{t + Δ t} + {B_{a, j}}^{t + Δ t}}{2}, {μ_{a, i j}}^{t + Δ t} = \frac{{μ_{a, i}}^{t + Δ t} + {μ_{a, j}}^{t + Δ t}}{2}

(C-3)

Eq. C-1 can be further approximated as:

.

\int_{V_{i}} \nabla \cdot (\frac{k k_{r a}}{B_{a} μ_{a}} \nabla p_{a}) d Ω = \int_{V_{i}} \sum_{j = 1}^{n_{i}} [\frac{k_{r a, i j}^{t + Δ t}}{B_{a, i j}^{t + Δ t} μ_{a, i j}^{t + Δ t}} [k_{i j}^{t + Δ t} V_{i} (e_{3, j}^{i} + e_{4, j}^{i})] (p_{a, (i, j)}^{t + Δ t} - p_{a, i}^{t + Δ t})] d Ω

(C-4)

For node i, if node i has no corresponding fracture node, then

.

\int_{V_{i}} q_{a, m f} d Ω = 0, a = o, w

(C-5)

If the matrix node i has a corresponding fracture node, suppose the serial number of the fracture node as f, then

,

\int_{V_{i}} q_{a, m f} d Ω = Q_{a, m f} = \frac{{k_{r a, i f}}^{t + Δ t}}{{B_{a, i f}}^{t + Δ t} {μ_{a, i f}}^{t + Δ t}} T_{i f} ({p_{a, f}}^{t + Δ t} - {p_{a, i}}^{t + Δ t}), a = o, w

(C-6)

where $T_{i f}$ is the transmissibility between the matrix node i and the corresponding fracture f. The schemes of phase permeability, volume coefficient, and viscosity are the same as those in Eqs. C-2 and C-3.

For the third term at the left hand side of Eq. 2, if there is no production well or injection well at node i, then:

.

\int_{V_{i}} q_{a, w e l l} d Ω = Q_{a, w e l l} = 0, a = o, w

(C-7)

If there are production wells or injection wells at matrix node i, one of the original motives of the EFVM proposed by Rao (2022b) is to handle the singular point source term in the meshless framework of GFDM, while the node control volume defined in EFVM is similar to the grid volume in the traditional grid-based FVM. However, in meshless EFVM, there is no specific description of the geometric shape of the node control domain, but only the value of the node control volume. Therefore, combined with the Peaceman formula, in MFDFM, the well index is calculated as:

W I_{i} = \frac{2 π k_{i} h_{i}}{\ln (r_{e, i} / r_{w}) + s}, r_{e, i} = 0.14 \sqrt{2 V_{i} / h_{i}},

(C-8)

where $r_{e, i}$ is the equivalent radius, $r_{w}$ is the well radius, and $s$ is the skin factor.

If the well is a production well, then

.

\int_{V_{i}} q_{a, w e l l} d Ω = Q_{a, w e l l} = \frac{{k_{r a, i f}}^{t + Δ t}}{{B_{a, i f}}^{t + Δ t} {μ_{a, i f}}^{t + Δ t}} W I_{i} ({p_{w f}}^{t + Δ t} - {p_{a, i}}^{t + Δ t}), a = o, w

(C-9)

If the well is a water injection well , then

,

\int_{V_{i}} q_{o, w e l l} d Ω = Q_{o, w e l l} = 0

(C-10)

.

\int_{V_{i}} q_{w, w e l l} d Ω = Q_{w, w e l l} = \frac{1}{B_{w, i f}^{t + Δ t}} (\frac{k_{r o, i f}^{t + Δ t}}{μ_{o, i f}^{t + Δ t}} + \frac{k_{r w, i f}^{t + Δ t}}{μ_{w, i f}^{t + Δ t}}) W I_{i} (p_{w f}^{t + Δ t} - p_{w, i}^{t + Δ t})

(C-11)

For the right hand side of Eq. 2, after integration, it is obtained that:

.

\int_{V_{i}} \frac{\partial (ϕ S_{a})}{\partial t} d Ω = \frac{V_{i}}{Δ t} [(\frac{ϕ_{i}^{t + Δ t} S_{a, i}^{t + Δ t}}{B_{a, i}^{t + Δ t}}) - (\frac{ϕ_{i}^{t} S_{a, i}^{t}}{B_{a, i}^{t}})]

(C-12)

Based on the above results, the discrete scheme of Eq. 2 can be obtained as follows:

\begin{array}{l} \sum_{j = 1}^{n_{i}} [\frac{{k_{i j}}^{t + Δ t} {k_{r a, i j}}^{t + Δ t}}{{B_{a, i j}}^{t + Δ t} {μ_{a, i j}}^{t + Δ t}} V_{i} (e_{3, j}^{i} + e_{4, j}^{i}) ({p_{a, (i, j)}}^{t + Δ t} - {p_{a, i}}^{t + Δ t})] + f l a g (i) \frac{{k_{r a, i f}}^{t + Δ t} {T_{i f}}^{t + Δ t}}{{B_{a, i f}}^{t + Δ t} {μ_{a, i f}}^{t + Δ t}} ({p_{a, f}}^{t + Δ t} - {p_{a, i}}^{t + Δ t}) \\ + Q_{a, w e l l} = \frac{V_{i}}{Δ t} [(\frac{{ϕ_{i}}^{t + Δ t} {S_{a, i}}^{t + Δ t}}{{B_{a, i}}^{t + Δ t}}) - (\frac{{ϕ_{i}}^{t} {S_{a, i}}^{t}}{{B_{a, i}}^{t}})], \end{array}

(C-13)

where it should be emphasized that, $n_{i}$ is the number of real nodes in the influence domain of matrix node i, $(i, j)$ refers to the sequence number of the j-th real node in global nodes except node i itself in the influence domain of matrix node i. $f l a g (i)$ is the function to judge whether the matrix node i has corresponding fracture nodes, if yes, $f l a g (i) = 1$ ⁠, if not, $f l a g (i) = 0$ ⁠.

Article History

Original SPE manuscript received for review 16 December 2021. Revised manuscript received for review 20 May 2022. Paper (SPE 210581) peer approved 24 May 2022.

References

Awadalla

,

T.

and

Voskov

,

D

.

2018

.

Modeling of Gas Flow in Confined Formations at Different Scales

.

Fuel

234

:

1354

–

1366

. 10.1016/j.fuel.2018.08.008.

Google Scholar

Crossref

Benito

,

J. J.

,

Ureña

,

F.

, and

Gavete

,

L

.

2001

.

Influence of Several Factors in the Generalized Finite Difference Method

.

Appl Math Model

25

(

12

):

1039

–

1053

. 10.1016/S0307-904X(01)00029-4.

Google Scholar

Crossref

Benito

,

J. J.

,

Ureña

,

F.

,

Gavete

,

L

. et al. .

2003

.

An H-Adaptive Method in the Generalized Finite Differences

.

Comput Methods Appl Mech Eng

192

(

5–6

):

735

–

759

. 10.1016/S0045-7825(02)00594-7.

Google Scholar

Bogomolny

,

A

.

1985

.

Fundamental Solutions Method for Elliptic Boundary Value Problems

.

SIAM J Numer Anal

22

(

4

):

644

–

669

. 10.1137/0722040.

Google Scholar

Crossref

Chavent

,

G.

and

Roberts

,

J. E

.

1991

.

A Unified Physical Presentation of Mixed, Mixed-Hybrid Finite Elements and Standard Finite Difference Approximations for the Determination of Velocities in Waterflow Problems

.

Adv Water Resour

14

(

6

):

329

–

348

. 10.1016/0309-1708(91)90020-O.

Google Scholar

Crossref

Cheng

,

A. H. D.

and

Cheng

,

D. T

.

2005

.

Heritage and Early History of the Boundary Element Method

.

Eng Anal Bound Elem

29

(

3

):

268

–

302

. 10.1016/j.enganabound.2004.12.001.

Google Scholar

Crossref

Fang

,

S.

,

Cheng

,

L.

, and

Ayala

,

L. F

.

2017

.

A Coupled Boundary Element and Finite Element Method for the Analysis of Flow through Fractured Porous Media

.

J Pet Sci Eng

152

:

375

–

390

. 10.1016/j.petrol.2017.02.020.

Google Scholar

Crossref

Fu

,

Z. J.

,

Tang

,

Z. C.

,

Zhao

,

H. T

. et al. .

2019

.

Numerical Solutions of the Coupled Unsteady Nonlinear Convection-Diffusion Equations Based on Generalized Finite Difference Method

.

Eur Phys J Plus

134

(

6

):

1

–

20

. 10.1140/epjp/i2019-12786-7.

Google Scholar

Crossref

Fu

,

Z. J.

,

Xie

,

Z. Y.

,

Ji

,

S. Y

. et al. .

2020

b.

Meshless Generalized Finite Difference Method for Water Wave Interactions with Multiple-Bottom-Seated-Cylinder-Array Structures

.

Ocean Eng

195

:

106736

. 10.1016/j.oceaneng.2019.106736.

Google Scholar

Crossref

Fu

,

Z. J.

,

Li

,

A. L.

,

Zhang

,

C

. et al. .

2020

a.

A Localized Meshless Collocation Method for Bandgap Calculation of Anti-Plane Waves in 2D Solid Phononic Crystals

.

Eng Anal Bound Elem

119

:

162

–

182

. 10.1016/j.enganabound.2020.07.014.

Google Scholar

Crossref

Chen

,

S. Y.

,

Hsu

,

K. C.

, and

Fan

,

C. M

.

2021

.

Improvement of Generalized Finite Difference Method for Stochastic Subsurface Flow Modeling

.

J Comput Phys

429

:

110002

. 10.1016/j.jcp.2020.110002.

Google Scholar

Crossref

Cheng

,

M.

and

Liu

,

G. R

.

2002

.

A Novel Finite Point Method for Flow Simulation

.

Int J Numer Meth Fluids

39

(

12

):

1161

–

1178

. 10.1002/fld.365.

Google Scholar

Crossref

Edwards

,

M. G

.

2011

.

Multi-Dimensional Wave-Oriented Upwind Schemes with Reduced Cross-Wind Diffusion for Flow in Porous Media

.

Int. J. Numer. Meth. Fluids

67

(

1

):

33

–

57

. 10.1002/fld.2329.

Google Scholar

Crossref

Dachanuwattana

,

S.

,

Jin

,

J.

,

Zuloaga-Molero

,

P

. et al. .

2018

.

Application of Proxy-Based MCMC and EDFM to History Match a Vaca Muerta Shale Oil Well

.

Fuel

220

:

490

–

502

. 10.1016/j.fuel.2018.02.018.

Google Scholar

Crossref

Gavete

,

L.

,

Gavete

,

M. L.

, and

Benito

,

J. J

.

2003

.

Improvements of Generalized Finite Difference Method and Comparison with Other Meshless Method

.

Appl Math Model

27

(

10

):

831

–

847

. 10.1016/S0307-904X(03)00091-X.

Google Scholar

Crossref

Gu

,

Y.

,

Qu

,

W.

,

Chen

,

W

. et al. .

2019

.

The Generalized Finite Difference Method for Long-Time Dynamic Modeling of Three-Dimensional Coupled Thermoelasticity Problems

.

J Comput Phys

384

:

42

–

59

. 10.1016/j.jcp.2019.01.027.

Google Scholar

Crossref

Gu

,

Y.

and

Sun

,

H. G

.

2020

.

A Meshless Method for Solving Three-Dimensional Time Fractional Diffusion Equation with Variable-Order Derivatives

.

Appl Math Model

78

:

539

–

549

. 10.1016/j.apm.2019.09.055.

Google Scholar

Crossref

Hajibeygi

,

H.

,

Karvounis

,

D.

, and

Jenny

,

P

.

2011

.

A Hierarchical Fracture Model for the Iterative Multiscale Finite Volume Method

.

J Comput Phys

230

(

24

):

8729

–

8743

. 10.1016/j.jcp.2011.08.021.

Google Scholar

Crossref

HosseiniMehr

,

M.

,

Tomala

,

J. P.

,

Vuik

,

C

. et al. .

2022

.

Projection-Based Embedded Discrete Fracture Model (PEDFM) for Flow and Heat Transfer in Real-Field Geological Formations with Hexahedral Corner-Point Grids

.

Adv Water Resour

159

:

104091

. 10.1016/j.advwatres.2021.104091.

Google Scholar

Crossref

Jiang

,

J.

and

Younis

,

R. M

.

2017

.

An Improved Projection-Based Embedded Discrete Fracture Model (PEDFM) for Multiphase Flow in Fractured Reservoirs

.

Adv Water Resour

109

:

267

–

289

. 10.1016/j.advwatres.2017.09.017.

Google Scholar

Crossref

Karimi-Fard

,

M.

and

Firoozabadi

,

A

.

2001

.

Numerical Simulation of Water Injection in 2D Fractured Media Using Discrete-Fracture Model

.

Paper presented at the

SPE Annual Technical Conference and Exhibition

,

New Orleans, Louisiana, USA

, 30 September–3 October. SPE-71615-MS. 10.2118/71615-MS.

Google Scholar

Crossref

Karimi-Fard

,

M.

,

Durlofsky

,

L. J.

, and

Aziz

,

K

.

2004

.

An Efficient Discrete-Fracture Model Applicable for General-Purpose Reservoir Simulators

.

SPE J.

9

(

2

):

227

–

236

. SPE-88812-PA. 10.2118/88812-PA.

Google Scholar

Crossref

Karimi-Fard

,

M.

,

Gong

,

B.

, and

Durlofsky

,

L. J

.

2006

.

Generation of Coarse-Scale Continuum Flow Models from Detailed Fracture Characterizations

.

Water Resour Res

42

(

10

). 10.1029/2006WR005015.

Google Scholar

Karimi-Fard

,

M.

and

Durlofsky

,

L. J

.

2016

.

A General Gridding, Discretization, and Coarsening Methodology for Modeling Flow in Porous Formations with Discrete Geological Features

.

Adv Water Resour

96

:

354

–

372

. 10.1016/j.advwatres.2016.07.019.

Google Scholar

Crossref

Li

,

L.

and

Lee

,

S. H

.

2008

.

Efficient Field-Scale Simulation of Black Oil in a Naturally Fractured Reservoir Through Discrete Fracture Networks and Homogenized Media

.

SPE Res Eval & Eng

11

(

4

):

750

–

758

. SPE-103901-PA. 10.2118/103901-PA.

Google Scholar

Crossref

Li

,

L.

and

Voskov

,

D

.

2021

.

A Novel Hybrid Model for Multiphase Flow in Complex Multi-Scale Fractured Systems

.

J Pet Sci Eng

203

:

108657

. 10.1016/j.petrol.2021.108657.

Google Scholar

Crossref

Li

,

P. W

.

2021

.

Space–Time Generalized Finite Difference Nonlinear Model for Solving Unsteady Burgers’ Equations

.

Appl Math Lett

114

:

106896

. 10.1016/j.aml.2020.106896.

Google Scholar

Crossref

Li

,

P. W.

and

Fan

,

C. M

.

2017

.

Generalized Finite Difference Method for Two-Dimensional Shallow Water Equations

.

Eng Anal Bound Elem

80

:

58

–

71

. 10.1016/j.enganabound.2017.03.012.

Google Scholar

Crossref

Liu

,

G. R.

and

Gu

,

Y. T

.

2005

.

An Introduction to Meshfree Methods and Their Programming

.

Berlin/Heidelberg, Germany

:

Springer Science & Business Media

.

Google Scholar

Li

,

B.

and

Tchelepi

,

H. A

.

2015

.

Nonlinear Analysis of Multiphase Transport in Porous Media in the Presence of Viscous, Buoyancy, and Capillary Forces

.

J Comput Phys

297

:

104

–

131

. 10.1016/j.jcp.2015.04.057.

Google Scholar

Crossref

Li

,

J.

,

Yu

,

W.

,

Guerra

,

D

. et al. .

2018

.

Modeling Wettability Alteration Effect on Well Performance in Permian Basin with Complex Fracture Networks

.

Fuel

224

:

740

–

751

. 10.1016/j.fuel.2018.03.059.

Google Scholar

Crossref

Matthäi

,

S. K.

,

Mezentsev

,

A.

, and

Belayneh

,

M

.

2005

.

Control-Volume Finite-Element Two-Phase Flow Experiments with Fractured Rock Represented by Unstructured 3D Hybrid Meshes

.

Paper presented at the

SPE Reservoir Simulation Symposium

,

The Woodlands, Texas, USA

, 31 January–2 February. SPE-93341-MS. 10.2118/93341-MS.

Google Scholar

Crossref

Michel

,

I.

,

Seifarth

,

T.

,

Kuhnert

,

J

. et al. .

2021

.

A Meshfree Generalized Finite Difference Method for Solution Mining Processes

.

Comp Part Mech

8

(

3

):

561

–

574

. 10.1007/s40571-020-00353-2.

Google Scholar

Crossref

Milewski

,

S

.

2012

.

Meshless Finite Difference Method with Higher Order Approximation—Applications in Mechanics

.

Arch Computat Methods Eng

19

(

1

):

1

–

49

. 10.1007/s11831-012-9068-y.

Google Scholar

Crossref

Moinfar

,

A.

,

Varavei

,

A.

,

Sepehrnoori

,

K

. et al. .

2013

.

Development of an Efficient Embedded Discrete Fracture Model for 3D Compositional Reservoir Simulation in Fractured Reservoirs

.

SPE J.

19

(

2

):

289

–

303

. SPE-154246-PA. 10.2118/154246-PA.

Google Scholar

Crossref

Monteagudo

,

J. E. P.

and

Firoozabadi

,

A

.

2004

.

Control-Volume Method for Numerical Simulation of Two-Phase Immiscible Flow in Two- and Three-Dimensional Discrete-Fractured Media

.

Water Resour Res

40

(

7

). 10.1029/2003WR002996.

Google Scholar

Noorishad

,

J.

and

Mehran

,

M

.

1982

.

An Upstream Finite Element Method for Solution of Transient Transport Equation in Fractured Porous Media

.

Water Resour Res

18

(

3

):

588

–

596

. 10.1029/WR018i003p00588.

Google Scholar

Crossref

Olorode

,

O.

,

Wang

,

B.

, and

Rashid

,

H. U

.

2020

.

Three-Dimensional Projection-Based Embedded Discrete-Fracture Model for Compositional Simulation of Fractured Reservoirs

.

SPE J.

25

(

4

):

2143

–

2161

. SPE-201243-PA. 10.2118/201243-PA.

Google Scholar

Crossref

Pruess

,

K.

and

Narasimhan

,

T. N

.

1985

.

A Practical Method for Modeling Fluid and Heat Flow in Fractured Porous Media

.

SPE J.

25

(

1

):

14

–

26

. SPE-10509-PA. 10.2118/10509-PA.

Google Scholar

Qu

,

W.

and

He

,

H

.

2020

.

A Spatial–Temporal GFDM with an Additional Condition for Transient Heat Conduction Analysis of FGMs

.

Appl Math Lett

110

:

106579

. 10.1016/j.aml.2020.106579.

Google Scholar

Crossref

Rao

,

B. N.

and

Rahman

,

S

.

2000

.

An Efficient Meshless Method for Fracture Analysis of Cracks

.

Comput Mech

26

(

4

):

398

–

408

. 10.1007/s004660000189.

Google Scholar

Crossref

Rao

,

X.

,

Cheng

,

L.

,

Cao

,

R

. et al. .

2019

a.

A Mimetic Green Element Method

.

Eng Anal Bound Elem

99

:

206

–

221

. 10.1016/j.enganabound.2018.12.004.

Google Scholar

Crossref

Rao

,

X.

,

Cheng

,

L.

,

Cao

,

R

. et al. .

2019

b.

A Modified Embedded Discrete Fracture Model to Improve the Simulation Accuracy During Early-Time Production of Multi-Stage Fractured Horizontal Well

.

Paper presented at the

SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition

,

Bali, Indonesia

, 29–31 October. SPE-196263-MS. 10.2118/196263-MS.

Google Scholar

Crossref

Rao

,

X.

,

Cheng

,

L.

,

Cao

,

R

. et al. .

2020

.

A Modified Projection-Based Embedded Discrete Fracture Model (PEDFM) for Practical and Accurate Numerical Simulation of Fractured Reservoir

.

J Pet Sci Eng

187

:

106852

. 10.1016/j.petrol.2019.106852.

Google Scholar

Crossref

Rao

,

X

.

2022

a.

An Upwind Generalized Finite Difference Method (GFDM) for Meshless Analysis of Heat and Mass Transfer in Porous Media

(accepted, in press).

Comput Part Mech

.

Google Scholar

Rao

,

X

.

2022

b.

A Novel Meshless Method Based on the Virtual Construction of Node Control Domains for Porous Flow Problems

. 10.48550/arXiv.2206.05531.

Google Scholar

Rao

,

X.

and

Liu

,

Y

.

2022

.

A Numerical Modelling Method of Fractured Reservoirs with Embedded Meshes and Topological Fracture Projection Configurations

.

Comput Model Eng Sci

131

(

3

):

1403

–

1429

. 10.32604/cmes.2022.018879.

Google Scholar

Rao

,

X.

,

Xin

,

L.

,

He

,

Y

. et al. .

2022

b.

Numerical Simulation of Two-Phase Heat and Mass Transfer in Fractured Reservoirs Based on Projection-Based Embedded Discrete Fracture Model (PEDFM)

.

J Pet Sci Eng

208

:

109323

. 10.1016/j.petrol.2021.109323.

Google Scholar

Crossref

Rao

,

X.

,

Liu

,

Y.

, and

Zhao

,

H

.

2022

a.

An Upwind Generalized Finite Difference Method for Meshless Solution of Two-Phase Porous Flow Equations

.

Eng Anal Bound Elem

137

:

105

–

118

. 10.1016/j.enganabound.2022.01.013.

Google Scholar

Crossref

Ren

,

G.

,

Jiang

,

J.

, and

Younis

,

R. M

.

2018

.

A Model for Coupled Geomechanics and Multiphase Flow in Fractured Porous Media Using Embedded Meshes

.

Adv Water Resour

122

(

DEC

):

113

–

130

. 10.1016/j.advwatres.2018.09.017.

Google Scholar

Ren

,

G.

and

Younis

,

R. M

.

2021

.

An Integrated Numerical Model for Coupled Poro-Hydro-Mechanics and Fracture Propagation Using Embedded Meshes

.

Comput Methods Appl Mech Eng

376

(

1–4

):

113606

. 10.1016/j.cma.2020.113606.

Google Scholar

Snow

,

D. T

.

1968

.

Rock Fracture Spacings, Openings, and Porosities

.

J Soil Mech and Found Div

94

(

1

):

73

–

91

. 10.1061/JSFEAQ.0001097.

Google Scholar

Crossref

Salam

,

N.

,

Haddade

,

A.

,

Clements

,

D. L

. et al. .

2017

.

A Boundary Element Method for A Class of Elliptic Boundary Value Problems of Functionally Graded Media

.

Eng Anal Bound Elem

84

:

186

–

190

. 10.1016/j.enganabound.2017.08.017.

Google Scholar

Crossref

Slough

,

K. J.

,

Sudicky

,

E. A.

, and

Forsyth

,

P. A

.

1999

.

Grid Refinement for Modeling Multiphase Flow in Discretely Fractured Porous Media

.

Adv Water Resour

23

(

3

):

261

–

269

. 10.1016/S0309-1708(99)00009-3.

Google Scholar

Crossref

Shakiba

,

M.

,

Cavalcante Filho

,

J. S. de A.

, and

Sepehrnoori

,

K

.

2018

.

Using Embedded Discrete Fracture Model (EDFM) in Numerical Simulation of Complex Hydraulic Fracture Networks Calibrated by Microseismic Monitoring Data

.

J Nat Gas Sci Eng

55

:

495

–

507

. 10.1016/j.jngse.2018.04.019.

Google Scholar

Crossref

Tene

,

M.

,

Al Kobaisi

,

M. S.

, and

Hajibeygi

,

H

.

2016

.

Algebraic Multiscale Method for Flow in Heterogeneous Porous Media with Embedded Discrete Fractures (F-AMS)

.

J Comput Phys

321

:

819

–

845

. 10.1016/j.jcp.2016.06.012.

Google Scholar

Crossref

Tene

,

M.

,

Bosma

,

S. B. M.

,

Al Kobaisi

,

M. S.

et al. .

2017

.

Projection-Based Embedded Discrete Fracture Model (PEDFM)

.

Adv Water Resour

105

:

205

–

216

. 10.1016/j.advwatres.2017.05.009.

Google Scholar

Crossref

Ureña

,

F.

,

Salete

,

E.

,

Benito

,

J. J

. et al. .

2012

.

Solving Third- and Fourth-Order Partial Differential Equations Using GFDM: Application to Solve Problems of Plates

.

Int J Comput Math

89

(

3

):

366

–

376

. 10.1080/00207160.2011.587871.

Google Scholar

Crossref

Voskov

,

D. V.

and

Tchelepi

,

H. A

.

2011

.

Compositional Nonlinear Solver Based on Trust Regions of the Flux Function Along Key Tie-Lines

.

Paper presented at the

SPE Reservoir Simulation Symposium

,

The Woodlands, Texas, USA

, 21–23 February. SPE-141743-MS. 10.2118/141743-MS.

Google Scholar

Crossref

Voskov

,

D. V.

and

Tchelepi

,

H. A

.

2012

.

Comparison of Nonlinear Formulations for Two-Phase Multi-Component EoS Based Simulation

.

J Pet Sci Eng

82–83

:

101

–

111

. 10.1016/j.petrol.2011.10.012.

Google Scholar

Wang

,

C.

,

Winterfeld

,

P.

,

Johnston

,

B

. et al. .

2020

.

An Embedded 3D Fracture Modeling Approach for Simulating Fracture-Dominated Fluid Flow and Heat Transfer in Geothermal Reservoirs

.

Geothermics

86

:

101831

. 10.1016/j.geothermics.2020.101831.

Google Scholar

Crossref

Wang

,

Y.

,

Gu

,

Y.

,

Fan

,

C. M

. et al. .

2018

.

Domain-Decomposition Generalized Finite Difference Method for Stress Analysis in Multi-Layered Elastic Materials

.

Eng Anal Bound Elem

94

:

94

–

102

. 10.1016/j.enganabound.2018.06.006.

Google Scholar

Crossref

Warren

,

J. E.

and

Root

,

P. J

.

1963

.

The Behavior of Naturally Fractured Reservoirs

.

SPE J.

3

(

3

):

245

–

255

. SPE-426-PA. 10.2118/426-PA.

Google Scholar

Wu

,

Y.

,

Cheng

,

L.

,

Fang

,

S

. et al. .

2020

.

A Green Element Method-Based Discrete Fracture Model for Simulation of the Transient Flow in Heterogeneous Fractured Porous Media

.

Adv Water Resour

136

:

103489

. 10.1016/j.advwatres.2019.103489.

Google Scholar

Crossref

Wu

,

Y. S.

,

Liu

,

H. H.

, and

Bodvarsson

,

G. S

.

2004

.

A Triple-Continuum Approach for Modeling Flow and Transport Processes in Fractured Rock

.

J Contam Hydrol

73

(

1–4

):

145

–

179

. 10.1016/j.jconhyd.2004.01.002.

Google Scholar

PubMed

Xia

,

H.

and

Gu

,

Y

.

2021

.

Short Communication: The Generalized Finite Difference Method for Electroelastic Analysis of 2D Piezoelectric Structures

.

Eng Anal Bound Elem

124

:

82

–

86

. 10.1016/j.enganabound.2020.12.012.

Google Scholar

Crossref

Xu

,

Y.

,

Fernandes

,

B. R. B.

,

Marcondes

,

F

. et al. .

2019

.

Embedded Discrete Fracture Modeling for Compositional Reservoir Simulation Using Corner-Point Grids

.

J Pet Sci Eng

177

:

41

–

52

. 10.1016/j.petrol.2019.02.024.

Google Scholar

Crossref

Yan

,

X.

,

Huang

,

Z.

,

Yao

,

J

. et al. .

2016

.

An Efficient Embedded Discrete Fracture Model Based on Mimetic Finite Difference Method

.

J Pet Sci Eng

145

:

11

–

21

. 10.1016/j.petrol.2016.03.013.

Google Scholar

Crossref

Yan

,

X.

,

Huang

,

Z.

,

Yao

,

J

. et al. .

2018

.

An Efficient Numerical Hybrid Model for Multiphase Flow in Deformable Fractured-Shale Reservoirs

.

SPE J.

23

(

4

):

1412

–

1437

. SPE-191122-PA. 10.2118/191122-PA.

Google Scholar

Crossref

Younis

,

R

.

2011

.

Modern Advances in Software and Solution Algorithms for Reservoir Simulation

.

PhD dissertation

,

Stanford University

,

Stanford, California

.

Yu

,

W.

,

Xu

,

Y.

,

Liu

,

M

. et al. .

2018

.

Simulation of Shale Gas Transport and Production with Complex Fractures Using Embedded Discrete Fracture Model

.

AIChE J

64

(

6

):

2251

–

2264

. 10.1002/aic.16060.

Google Scholar

Crossref

Zhang

,

N.

,

Yan

,

B.

,

Sun

,

Q

. et al. .

2017

c.

Improving Multiscale Mixed Finite Element Method for Flow Simulation in Highly Heterogeneous Reservoir Using Adaptivity

.

J Pet Sci Eng

154

:

382

–

388

. 10.1016/j.petrol.2017.04.012.

Google Scholar

Crossref

Zhang

,

Q.

,

Huang

,

Z.

,

Yao

,

J

. et al. .

2017

a.

A Multiscale Mixed Finite Element Method with Oversampling for Modeling Flow in Fractured Reservoirs Using Discrete Fracture Model

.

J Comput Appl Math

323

:

95

–

110

. 10.1016/j.cam.2017.04.005.

Google Scholar

Crossref

Zhang

,

Q.

,

Huang

,

Z.

,

Yao

,

J

. et al. .

2017

b.

Multiscale Mimetic Method for Two-Phase Flow in Fractured Media Using Embedded Discrete Fracture Model

.

Adv Water Resour

107

:

180

–

190

. 10.1016/j.advwatres.2017.06.020.

Google Scholar

Crossref

Zhou

,

Y.

,

Tchelepi

,

H. A.

, and

Mallison

,

B. T

.

2011

.

Automatic Differentiation Framework for Compositional Simulation on Unstructured Grids with Multi-Point Discretization Schemes

.

Paper presented at the

SPE Reservoir Simulation Symposium

,

The Woodlands, Texas, USA

, 21–23 February. SPE-141592-MS. 10.2118/141592-MS.

Google Scholar

Crossref

Zhuang

,

X.

,

Augarde

,

C. E.

, and

Mathisen

,

K. M

.

2012

.

Fracture Modeling Using Meshless Methods and Level Sets in 3D: Framework and Modeling

.

Int J Numer Meth Engng

92

(

11

):

969

–

998

. 10.1002/nme.4365.

Google Scholar

Crossref

A Meshless Numerical Modeling Method for Fractured Reservoirs Based on Extended Finite Volume Method

Summary

Introduction

Basic Governing Equations and Numerical Methods

Basic Governing Equations for Flow in Fractured Reservoirs

Principles of GFDM

EFVM: Integral Form of GFDM

Definition of the Node Control Volume

Calculation of Node Control Volumes by Adding Virtual Nodes

Discretization of Governing Equations and Definition of Transmissibility

The MFDFM Based on EFVM

Meshless Point Cloud Discretization of the Fracture-Reservoir Domain

Discretization of Governing Equations

Definition of Connection Types between Nodes and Transmissibility Calculation

Treatment of Boundary Conditions and Construction of Global Equations

Numerical Examples

MFDFM Performances for Nonintersecting Fractures With a Matching Point Cloud

MFDFM Performances for a Complex Orthogonal Fracture Network With Matching and Nonmatching Point Clouds

MFDFM Performances for a Nonorthogonal Complex Fracture Network and a Mixed Boundary Condition

MFDFM Performances With Complex Reservoir-Boundary Geometry

Conclusions

Future Works

Author Contributions

Acknowledgments

Appendix A: An Algorithm for Generating a Nonmatching Point Cloud

Step 1: Calculating the Intersections between Fractures

Step 2: Generation of Fracture Nodes

Step 3: Generation of Matrix Nodes

Step 4: Construction of the NNCs between Matrix Nodes and Fracture Nodes

Appendix B: An Algorithm for Generating a Matching Point Cloud

Step 1: Calculating the Intersections between Fractures

Step 2: Generation of Fracture Nodes

Step 3: Generation of the Type-II Matrix Nodes

Step 4: Integration of the First Type of Matrix Nodes and the Second Type of Matrix Nodes

Appendix C: Derivation of the EFVM Discrete Scheme

Article History

References

Cited By

Email Alerts

Suggested Reading

Explore

Connect

Resources

Engage

Sign In or Register to Continue