Abstract
In this study, we present a hybrid model coupled with two-set nodes Green element method (GEM) and embedded discrete fracture model (EDFM) for capturing the effect of transient flow in inhomogeneous fractured porous media. GEM is an excellent advanced algorithm, which can solve nonlinear problems in heterogeneous media. That is also an obvious advantage of GEM against the original boundary element method (BEM). The novel GEM has double nodes of pressure and flux and it is an improvement of classical GEM, which has a second-order precision and fits for triangle structured grids. In the place of adopting the linear flow approximation for original EDFM, the interflows between local triangle matrix grids and fracture elements are derived using the novel GEM, which has higher accuracy than those in previous EDFMs. Consequently, the modified hybrid model can indeed calculate the pressure and flux distribution of transient flow in multifracture porous media. Three numerical cases are presented to show the practicability of our novel model which include (i) multistage fractured horizontal well, (ii) heterogeneous fractured porous media, and (iii) complex fracture networks (CFNs) in an unconventional reservoir.
1. Introduction
Accurate modeling and simulation of pressure and flux distribution in naturally and hydraulically fractured reservoirs is a hot topic and quite challenging work for many oilfield engineering technicians. So far, there are lots of numerical approaches that have been employed to solve this physics problem, such as the equivalent continuum model [1–3], discrete fracture model (DFM or DFN) [4, 5], and embedded discrete fracture model (EDFM) [6–8]. Each approach has advantages and disadvantages, but the most widely used method is EDFM at present because of its high precision and high efficiency. The original EDFM was firstly presented by Li and Lee (2008) [9], in which discrete fractures are embedded in the structured matrix cell, and the mass transfer equations can be expressed by the nonneighboring connection (NNC) geometrical relationship of fractures and matrix. Compared with DFM, EDFM only adopts a set of fixed structured matrix blocks and avoids the complexity of grid division. Besides, it is easy to be coupled with the existing reservoir simulation codes based on the finite volume method (FVM). Therefore, EDFM is an important improvement of DFM and has attracted the attention of many scholars in recent years. But, for the original EDFM, computing the interflux conductivity between matrix grids and discrete fracture elements is only related to geometrical properties, which did not in accord with the actual physical fact. To overcome this problem, the boundary element method (BEM) is possibly a compelling treatment to calculate the mass transmission between matrix and fracture. However, the original BEM cannot deal with this problem effectively because the actual reservoir is always strongly heterogeneous. Consequently, it is essential to develop an effective numerical approach for coupled EDFM simulators to characterize the transient mass transfer in heterogeneous fractured porous medium.
Green function method (GEM), intrinsically a variation of BEM, can effectively deal with nonlinear problems in heterogeneous media. The original GEM was proposed by Taigbenu et al. (1995a, 1995b, 1997, 1999) [10–14]; in this method, the calculation domain was meshed by polygonal cells, and cell vertices were considered as unknown nodes. In previous GEMs, the normal flux at each internal node was obtained by distinguishing the pressure expressed by the weighted average of nodal values of the pressure and basis functions. However, this approximation approach reduced overall accuracy and errors will widen with the decrease of polygonal size. Then, some scholars put forward some improved methods for original GEM. Archer (1999, 2006) et al. [15, 16] pointed out the error of original GEM can be reduced by adopting the overhauser interpolation functions. Pecher et al. (2001) [17] and Lorinczi et al. (2006, 2008, 2010, 2011) [18–21] presented a modified GEM, in which they used the continuity of flux vector component to approximate the flux term and increased the number of degrees of freedom of each internal node. Lorinczi et al. do not use the traditional way of coupling the equations at the same source point into one equation. Instead, all the equations at the same source point are listed as independent equations to construct the overdetermined equations, and then, the unknowns are solved by the least squares method. Therefore, the pressure and flux terms are explicitly solved in this method, which improved to second-order accuracy. But, the system of equations may be singular in this method due to the excessive number of degrees of freedom of each internal node, which reduces the robustness of numerical simulation. Taigbenu (2008) [22] proposed a flux-correct GEM and thought that the sum of flux term on the cell edge in the clockwise direction (or both in the counterclockwise direction) around the same node is zero, and uses these additional equations to form a closed linear equation to solve the unknowns. Taigbenu (2012) [23] returned to the original GEM and improved the approximate accuracy of normal flux by using the quadratic interpolation basis function. But in essence, the flux term is not solved explicitly and the accuracy is not high enough. By referring to the calculation of the flux term on the edge of cell in the mixed finite element method (MFEM), Fang et al. (2017) [24] presented a novel mixed boundary element method (MBEM) to improve the accuracy by moving the pressure and flux nodes to the edge-center point to satisfy the local mass conservation. Rao et al. (2018a, 2018b) [25, 26] presented a novel GEM that incorporates double nodes of pressure and flux; in this method, the pressure unknowns included all cell vertices and flux unknowns contained all edge-center points. The advantage of this approach is that it satisfies the local mass conservation and has high accuracy, and it can be used for solving stably the second-order partial differential equations (PDEs), such as the mass transmission and heat transfer equations. To further improve the strong robustness of GEM, Rao et al. (2019) [27] proposed a mimetic Green element method (MGEM) which coupled the idea of the mimetic finite difference method (MFDM).
At present, no matter the Green element method (GEM) itself or the application in embedded discrete fracture model (EDFM) is still at the exploration stage. Because the GEM has the unique advantage of high accuracy in dealing with unsteady and heterogeneous problems, it is important to further investigate the novel GEM and its application in the engineering field. In this study, a modified EDFM based on novel GEM with two-set nodes triangular element was proposed to solve accurately the interflow between matrix and fractures, which can effectively capture the dynamic flow characteristics for heterogeneous fractured reservoirs. The paper is organized as follows. In Section 2, the principle of novel GEM and its coupling method in EDFM are illuminated. In Section 3, the model validation process is introduced and we compare the accuracy of early-time and long-time performance of novel model and original method. In Section 4, three cases are implemented to show the application of our novel model which include (i) multistage fractured horizontal well, (ii) heterogeneous fractured porous media, and (iii) complex fracture networks in unconventional reservoirs. Finally, the conclusions are presented in Section 5.
2. A Brief Review of Original EDFM
Firstly, we explain the concept of EDFM approach and illustrate the reason why the accuracy of transient simulation of original EDFM is not perfect for approximating the interflow between matrix cell and fracture elements. According to the pretreatment of EDFM, the fractured porous medium can be independently divided into a matrix cell system and some discrete fracture elements as plotted in Figure 1. For two-dimensional media, fracture elements are regarded as 1D pathways discretized independently from matrix cells. Besides, the matrix and fracture elements are illustrated on top of each other with overlapped matrix grids highlighted with the color blue.
In previous EDFMs, the mass transfer between matrix and fracture is approximated by the geometric index in the transmissibility term, which is based on a steady flow. Li and Lee [9] firstly assume that the pressure changes linearly along each fracture’s vertical direction in the grid block and adopt the following equations to calculate the interflux between fracture and matrix as shown in Equation (1) and Equation (2), which has been extended to many models and studies [28–30]. However, this linear flow assumption that the pressure linearly distributed along each fracture’s vertical direction in matrix cell is not reasonable in the production process, especially in the early stage of production, resulting in low precision in the early stage owing to the great difference between ultralow matrix permeability and ultrahigh fracture permeability.where is the geometric index between matrix grid and fracture control volume, is the surface area of fracture element connected to matrix, is the average distance between the matrix grid and the fracture grid, denotes the distance between fracture and matrix grid, and indicates the matrix grid volume.
3. Methodology of Novel EDFM
In this report, we explain the idea of GEM of double nodes of pressure and flux combined with the flow equation in porous medium. The single-phase fluid flow in matrix is generally described as follows:where represents the pressure over matrix domain, indicates the distribution of internal source or sink strengths and generally equal to 0, is the matrix permeability, and denotes the medium attributes.
The above flow equation also can be modified as follows:where , , and .
As shown in Figure 2, the novel GEM has double nodes of pressure and flux, in which pressure nodes are distributed at the three vertices of triangular matrix element and flux nodes are distributed at the edge-center points of triangular matrix block. The pressure nodes are marked as 1, 2, and 3, respectively. The flux nodes are marked as , , and , respectively.
Then, the boundary integral equation of Equation (4) is adopted to each triangular element as follows:
For each pressure point, the corresponding basis functions of the pressure nodes are expressed as , respectively. The parameters of pressure within triangle element are obtained by the weighted average of the nodal values and basis functions. For each flux node, owing to the piecewise continuous of normal flux, normal flux at one edge is considered as a constant, which is obviously different from the previous GEMs.
When a pressure node is chosen as a source term, the boundary integral equation of Equation (5) is then rewritten as follows:where
Then, the system of linear equations of triangular block can be approximated as follows:where
Based on Equation (9), the inverse matrix of coefficient matrix can be solved, denoted . The two sides of Equation (9) are multiplied by ; then, it is obtained that
The above expression can be rewritten in an implicit form as follows:
Further, make the above equation concise:where
There are physical quantities in the above formulas to be clearly defined. For an actual physical fact, the interflux between matrix and fracture is considered as the source or sink term in matrix cell, such as Equation (3). Consequently, it is necessary to clear out in which unit is a source term or sink term, which is related to the geometrical properties between fractures and matrix cells.
As shown in Figure 3, it is a sketch of two-dimensional EDFM, in which fractures are plotted in blue lines and divided into several cells by boundaries of matrix grids and intersecting fractures. There is one node in each fracture cell, and the node is positioned to the midpoint of the fracture cell.
Put the following cell as a case to explain the process of matrix cells which contain fracture cells. In the cell, in Equation (3) is substituted by ; then, the formulas of the grid containing fracture segments are expressed as follows:where is the interflux between matrix cell and fracture cell.
Since the apertures of fractures are really small, the source term over a fracture segment can be assumed identical. Then,where is the aperture of the th fracture cell, is the interflow between the matrix cell and the th fracture element, is the line of the th fracture element, and is the number of fracture elements.
Then, Equation (11) of the cell containing fracture segments is rewritten as follows:where
Equations (11) and (17) are equations of matrix cells with and without fracture segments contained, respectively, which are expressed as the first part of global equations . These equations are coupled by utilizing the physical reality that the sum of the normal fluxes at flux nodes on shared edge is zero, which is a main characteristic of GEM based on two-set nodes. As shown in Figure 4, the two adjacent triangular cells denoted by and are explained as an example to introduce the detail of the coupling process as follows.
May as well assume that equations of element are as follows:
Since the edge 2-3 is shared by these two cells, the equations of element are as follows:
Couple Equations (19) and (20) to the following equation:
Assuming the numbers of nodes of matrix cells, edges of matrix cells, and fracture segments are , , and , respectively, as stated by Archer et al. [15], the number of the first part of global equations is , which is always larger than , so is an overdetermined system.
The second part of global equations denoted by is obtained in the cells containing fracture segments. In each cell with fracture segments contained, source points are selected at midpoints containing fracture segments, respectively, equations can be obtained, and these expressions are as the second part of equations. Also, another needed equations denoted by the third part can be obtained by handling flow equations in fractures; specific details are as follows.
The finite difference method (FDM) can be applied to discretize the flow equation to obtain . As shown in Figure 5, it is a fracture divided into segments. To be coupled with and , using FDM in an implicit form shown in Equation (22), needed equations can be obtained, which is the same as the method proposed by Jia et al. (2015) [31].where , , , , and are the coefficients term, is the flux from the fracture cell to the well, represents that there is no well in the fracture cell. If there is a well in the fracture cell, it can be assumed that the well is producing at constant bottom hole pressure (BHP), i.e., and .
What is to be solved in a fracture cell with a well is and ; however, the unknowns of other fracture segments without a well are and . Consequently, there are two uncertain factors in each fracture segment. The problem of the constant production rate of production well can be solved similarly.
In all, there are equations in , , , and unknowns including nodal pressures of matrix cells and unknowns of fracture segments. Since is larger than , the whole global equations are overdetermined owing to the application of novel GEM. Adding the sealed outer boundary condition, the resultant equations may as well be assumed to be expressed as follows:
Calculating this overdetermined system is equivalent to the solution in Equation (24) by orthogonal projection theorem in functional analysis. Thus, all the unknowns can be obtained.
4. Model Validation and Discussion for Transient Effect
4.1. Model Validation
The quality of modified two-dimensional EDFM based on the novel GEM is compared with the original EDFM by Li and Lee (2008) [9] and the local grid refinement (LGR) by tNavigator®. The tNavigator® is a commercial high-performance reservoir simulation software developed by Rock Flow Dynamics (RFD) group, which results can be as the exact solution. The first case is to compare the accuracy of three simulators in capturing the long-time productivity of one horizontal well with one fracture. Values of physical properties of low-permeability oil reservoir are indicated in Table 1. Figure 6 shows the sketch of a rectangular reservoir domain (length: 1010 meters, width: 1010 meters). The reservoir thickness is 10 meters, and the initial reservoir pressure is 20 MPa. There is a 210-meter fracture and a production well in the center of formation. The horizontal wellbore crosses the midpoint (505, 505) of the fracture in a 2D drawing. Figure 7 illustrates the local refinement grids which can be used for representing one fracture in tNavigator®. The fracture can be explicitly described by EDFM but it must be along the direction of grid line in the LGR method. The horizontal well is simulated with a constant BHP of 10 MPa. Figure 8 shows the computational pressure distribution maps over 1000 days of three various simulators. Wellbore flow rate curves obtained from three various simulators are shown in Figure 9. It can be indicated that the results of three different simulators are in good agreement, which verifies the overall accuracy for long-time production performance of proposed method.
(a) The original EDFM
(b) The novel EDFM
(c) LGR by tNavigator®
4.2. Discussion for Transient Effect
This discussion is conducted to verify that modified EDFM based on novel GEM can reflect the early-time transient flow effect and to analyze the average relative error between the original EDFM and the modified EDFM. The basic properties of this case are the same as in Table 1. There are two important factors that are considered in the discussion: grid size and matrix permeability. Consequently, the dimension of a grid is set as and , and various matrix permeabilities include 0.1 mD, 1 mD, and 10 mD, respectively. The computing time is set as 10 days. To achieve an accurate response to the early transient effect, a smaller time step is used in the simulation process compared with the verification example. Comparison of pressure distribution maps over 10 days of three various simulators () is shown in Figure 10. It can be indicated from the results of pressure distribution maps that there was no significant difference between three different simulators. Wellbore flow rate curves obtained from three various simulators (original EDFM, modified EDFM, and tNavigator®) in a condition of various grid dimensions ( and ) and matrix permeability (0.1 mD, 1 mD, and 10 mD) are compared in Figure 11. Compared with the results of tNavigator®, it is apparent from the early-time results that the modified EDFM has higher precision and robustness than the original EDFM. This is because proposed novel GEM can effectively reflect the transient interflow between local triangle matrix grids and fracture elements replacing the linear flow approximation in the original EDFM. The degree of transient flow effect is closely related to grid dimension and matrix permeability. Average relative errors are obtained by Equation (25), and the corresponding results are indicated in Figure 12. It can be found that the average relative error of the results of modified EDFM is much lower than that of original EDFM. Also, with the decrease of matrix permeability and the increase of grid size, the effect of transient flow is magnifying such that the average relative errors are increasing. The proposed novel GEM indeed achieves a more clear characterization of transient flow effects and improves the early accuracy of simulation.where represents the wellbore flow rate calculated by EDFM, represents the wellbore flow rate of tNavigator®, and is the number of time steps.
(a) The original EDFM
(b) The novel EDFM
(c) LGR by tNavigator®
(a) Matrix permeability 0.1 mD and grid size 2 m
(b) Matrix permeability 0.1 mD and grid size 10 m
(c) Matrix permeability 1 mD and grid size 2 m
(d) Matrix permeability 1 mD and grid size 10 m
(e) Matrix permeability 10 mD and grid size 2 m
(f) Matrix permeability 10 mD and grid size 10 m
5. Numerical Examples
5.1. Multistage Fractured Horizontal Well
Technologies of multistaged fracturing and drilling horizontal wells are widely applied in the production of low-permeability oil reservoirs. As illustrated in Figure 13, it is a rectangular reservoir (length: 1000 meters, width: 600 meters). The thickness of reservoir is 10 meters. The outer boundary is sealed, and the initial reservoir pressure is 20 MPa. A horizontal well is evenly staged fractured into 11 fractures and it is producing at a constant BHP of 10 MPa. Other input parameters of this case are as same as in Table 1. Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days and prediction of well production performance over 1000 days are plotted in Figures 14 and 15, respectively. The simulation results indicated that the initial productivity of horizontal well is very high in the early time of production (the first 100 days) and the oil produced in the wellbore is mainly in fractures at this stage. After that, the formation pressure and daily oil production decrease rapidly and maintain at a low value. Therefore, the low-permeability reservoir development is dominated by the “fracture-controlled reserves.” When evaluating the productivity of fractured wells, the influence of half-length and number of hydraulic fractures are very significant.
(a) Oil production performance
(b) Average formation pressure and reserve recovery
5.2. Heterogeneous Fractured Porous Media
The prior example in Figure 13 is an ideal case, and the actual reservoir is heterogeneous. The proposed model has the obvious advantage against the original boundary element method (BEM) in solving nonlinear problems in heterogeneous media. Therefore, the above model is repeatedly simulated based on the actual reservoir properties. The new scheme retains the same parameters of hydraulic fracture as above and adds the random distribution of permeability and porosity, focusing on the evaluation of the reservoir heterogeneity influence on well performance. The matrix permeability distribution map is plotted in Figure 16. The matrix permeability median value is still 1 mD and its gradation ranges from 0.1 to 10 mD. The matrix porosity distribution map is illustrated in Figure 17, which shows that the matrix porosity median value is 0.17 and its gradation ranges from 0.1 to 0.25. Other physical parameters of this case are same as the example 1. Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days and prediction of well production performance over 1000 days are plotted in Figures 18 and 19, respectively. Compared with the prior example, it can be found from Figure 19 that the daily oil production and cumulative oil production of horizontal well in heterogeneous reservoir are less than that in homogeneous reservoir under the condition of keeping other factors unchanged. This is because the reservoir heterogeneity increases the fluids flow resistance. The application of novel EDFM based on two-set nodes GEM and triangle element in this example shows that the proposed model can effectively capture the production characteristic of fractured horizontal wells in actual oilfield.
(a) Oil production performance
(b) Average formation pressure and reserve recovery
5.3. Complex Fracture Networks in Unconventional Reservoir
Unconventional reservoirs, such as tight or shale oil reservoir, are characterized by extensive natural fractures in matrix and the reaction of complex fracture networks (CFNs) after hydraulic fracturing stimulation. As indicated in Figure 20, it is a rectangular oil reservoir (length: 2000 meters, width: 1200 meters). The reservoir is 10 meters thick and the outer boundary condition is sealed, and the initial reservoir pressure is 25 MPa. There are 30 fractures with different length and azimuth, which constitute the multiscale complex fracture networks in the formation. A production horizontal well is producing at a constant BHP of 10 MPa. Input parameters of tight oil reservoir are shown in Table 2. Pressure distribution maps in the 1st, 10th, and 100th days and 1st, 3rd, and 10th years are plotted in Figure 21. Prediction of well production performance over 10 years is illustrated in Figure 22. The computational results illustrated that the development of unconventional reservoir depends on the CFNs stimulated by hydraulic fracturing because the matrix permeability of unconventional reservoir is smaller than that of low-permeability reservoir. The scale and degree of CFNs is an important factor to affect the productivity of horizontal well. Overall, the inspiration from the perspective of engineering technicians is that both two aspects need to be considered to improve significantly the total production as much as possible. On the one hand, it is necessary to establish long-wide hydraulic fractures with high conductivity. On the other hand, it is important to increase the area and utilization rate of stimulated reservoir volume (SRV).
(a) Oil production performance
(b) Average formation pressure and reserve recovery
6. Conclusions
A modified EDFM by mixing the idea of two-set nodes GEM based on the triangle element is developed in this paper and it can effectively solve the transient effect of single phase in fractured heterogeneous reservoir. Owing to the novel GEM, the proposed approach can investigate the mass transfer between matrix and fracture more exactly. Besides, the novel model has the obvious advantage against the original BEM in solving nonlinear problems in heterogeneous medium. Some examples are illustrated to verify the precision of proposed model and show the application in fractured reservoirs.
Data Availability
The [DATA TYPE] data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We sincerely appreciate the support from the National Natural Science Foundation of China (No. U19B6003, No. U1762210, and No. 51774297), the National Science and Technology Major Project of China (No. 2017ZX05069), and the China Petroleum Science and Technology Program (Major Program, ZLZX2020-02-04).