# A Coupled Novel Green Element and Embedded Discrete Fracture Model for Simulation of Fluid Flow in Fractured Reservoir

Authors: Linsong Cheng; Xulin Du; Renyi Cao; Zhikai Wang; Junjie Shi
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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 Section2, 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 Figure1. 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.Figure 1
Schematic of EDFM concept (modified from Li and Lee (2008) [9]).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.(1)CI=Ad,(2)d=∫Vmxnx′dx′Vm,where CI is the geometric index between matrix grid and fracture control volume, A is the surface area of fracture element connected to matrix, d is the average distance between the matrix grid and the fracture grid, xn denotes the distance between fracture and matrix grid, and Vm 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:(3)∇⋅K∇p=c∂p∂t+f,where p represents the pressure over matrix domain, f indicates the distribution of internal source or sink strengths and generally equal to 0, K is the matrix permeability, and c denotes the medium attributes.The above flow equation also can be modified as follows:(4)∇2p=−∇ψ⋅∇p+σ∂p∂t+νf,where ψ=lnK, σ=cv, and ν=1/K.As shown in Figure2, 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 a, b, and c, respectively.Figure 2
One triangle cell for the novel GEM of double nodes of pressure and flux.Then, the boundary integral equation of Equation (4) is adopted to each triangular element as follows:(5)−λpi+∫Γep∇G⋅n+GqnKds+∬ΛeG−∇ψ⋅∇p+σ∂p∂t+νfdA=0.For each pressure point, the corresponding basis functions of the pressure nodes are expressed asϕ1,ϕ2,ϕ3, 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.(6)p=∑i=13piϕi.When a pressure nodei is chosen as a source term, the boundary integral equation of Equation (5) is then rewritten as follows:(7)∑j=13Rijpj+∑z=acLizqnKz−∑j=13∑m=13Uimjψmpj+∑j=13∑m=13Wimjσmdpjdt+υmfj=0,i,j=1,2,⋯,n,where(8)Rij=∫Γeϕj∇GM,Mi⋅nds−δijλ,Liz=∫ΓzGM,Mids,Γz=theedgeoffluxnode,Uimj=∬ΛeGM,Mi∂ϕm∂x∂ϕj∂xdA+GM,Mi∂ϕm∂y∂ϕj∂ydA,Wimj=∬ΛeGM,MiϕmϕjdA.Then, the system of linear equations of triangular block can be approximated as follows:(9)Biz3×3qn=−Eij3×3p−Cij3×3dpjdt−Fi3,where(10)Eij=∑m=13Rij−Uimjψm,Biz=Liz1Kz,Cij=∑m=13Wimjσm,Fi=∑j=13∑m=13Wimjυmfj.Inmatrixcellnotcontainingfractures,Fi=0,p=p1,p2,p3T,q=qA,qB,qCT.Based on Equation (9), the inverse matrix of coefficient matrix B can be solved, denoted B−1. The two sides of Equation (9) are multiplied by B−1; then, it is obtained that(11)qn=−B−1Ep−B−1Cdpdt−B−1F.The above expression can be rewritten in an implicit form as follows:(12)qn=−θB−1E+B−1CΔtpn+1+−1−θB−1E+B−1CΔtpn−B−1F.Further, make the above equation concise:(13)qn=Mpn+1+b,where(14)M=−θB−1E+B−1CΔt,b=−1−θB−1E+B−1CΔtpn−B−1F.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 Figure3, 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.Figure 3
Two-dimensional EDFM based on triangle elements and the novel GEM.Put the following cell as a case to explain the process of matrix cells which contain fracture cells. In the cell,f in Equation (3) is substituted by qomf; then, the formulas of the grid containing fracture segments are expressed as follows:(15)−λpi+∫Γep∇G⋅n+GqnKds+∬ΛeG−∇ψ⋅∇p+σ∂p∂t+νqomfdA=0,where qomf is the interflux between matrix cell and fracture cell.Since the apertures of fractures are really small, the source termqomf over a fracture segment f can be assumed identical. Then,(16)∬ΛeνGM;MiqomfdA=∑j=1nf∑m=13ωjqomf,j∬ΓjνmϕmGM;MidA,where ωj is the aperture of the jth fracture cell, qomf,j is the interflow between the matrix cell and the jth fracture element, Γj is the line of the jth fracture element, and nf is the number of fracture elements.Then, Equation (11) of the cell containing fracture segments is rewritten as follows:(17)qn=Biz3×3−1−Eij3×3p−Cij3×3dpjdt−Fi3,where(18)Eij=∑m=13Rij−Uimjψm,Biz=Liz1Kz,Cij=∑m=13Wimjσm,Fi=∑j=1nf∑m=13ωjqomf,j∬ΓjνmϕmGM;MidA,p=p1,p2,p3T,q=qA,qB,qCT.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 I1. 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 e1 and e2 are explained as an example to introduce the detail of the coupling process as follows.Figure 4
Coupling process for the novel GEM based on double nodes of pressure and flux.May as well assume that equations of elemente1 are as follows:(19)qa2qb2qc2=M11e1M12e1M13e1M21e1M22e1M23e1M31e1M32e1M33e1p12p22p32+b1e1b2e1b3e1.Since the edge 2-3 is shared by these two cells, the equations of elemente2 are as follows:(20)−qa2qd2qe2=M11e2M12e2M13e2M21e2M22e2M23e2M31e2M32e2M33e2p22p32p42+b1e2b2e2b3e2.Couple Equations (19) and (20) to the following equation:(21)0qb2qc2qd2qe2−b1e1+b1e2b2e1b3e1b2e2b3e2=M11e1M12e1+M12e2M13e1+M13e2M13e2M21e1M22e1M23e10M31e1M32e1M33e100M21e2M22e2M23e20M31e2M32e2M33e2p12p22p32p42.Assuming the numbers of nodes of matrix cells, edges of matrix cells, and fracture segments arenp, nedge, and nfe, respectively, as stated by Archer et al. [15], the number of the first part of global equations is nedge, which is always larger than np, so I1 is an overdetermined system.The second part of global equations denoted byI2 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, nfe equations can be obtained, and these expressions are as the second part of equations. Also, another needed nfe equations denoted by the third part I3 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 obtainI3. As shown in Figure 5, it is a fracture divided into nf segments. To be coupled with I1 and I2, using FDM in an implicit form shown in Equation (22), needed nf equations can be obtained, which is the same as the method proposed by Jia et al. (2015) [31].(22)Bi′pf,i−1n+1+Ci′−Fi′Δtpf,in+1+Di′pf,i+1n+1+Ei′qomf,i+Ei′qwell+Fi′Δtpf,in=0,where Bi′, Ci′, Di′, Ei′, and Fi′ are the coefficients term, qwell is the flux from the fracture cell to the well, qwell=0 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., pwf and pf=pwf.Figure 5
Sketch of divided fracture segments on a fracture.What is to be solved in a fracture cell with a well isqomf and qwell; however, the unknowns of other fracture segments without a well are qomf and pf. 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 arene+2nf equations in I1, I2, I3, and np+2nf unknowns including np nodal pressures pm of matrix cells and 2nf unknowns of fracture segments. Since ne is larger than np, 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:(23)bine+2nf=Mijne+2nf×np+2nfpin+1ne+2nf.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.(24)MTMpn+1=MTb.


## 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.Table 1
Input parameters of numerical case.
ParametersValueParametersValueMatrix porosity0.17Matrix permeability1 mDRock compressibility coefficient1.07×10−4MPa−1Fracture aperture10 mmFracture porosity0.4Fracture permeability20000 mDFluid viscosity2 mPa·sFluid compressibility coefficient3.02×10−3MPa−1Figure 6
Sketch of reservoir domain and a horizontal well with one fracture.Figure 7
Top view of locally refined fracture grids and their surrounding matrix grids.Figure 8
Comparison of pressure distribution maps over 1000 days of three various simulators.
(a)
The original EDFM(b)
The novel EDFM(c)
LGR by tNavigator®Figure 9
Comparison of wellbore flow rate over 1000 days of three various simulators.


### 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 Table1. 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 2×2 and 10×10meters, 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 (10×10m) 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 (2×2m and 10×10m) 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.(25)Ae=1N∑i=1Nqi−qirefqiref×100,where qi represents the wellbore flow rate calculated by EDFM, qiref represents the wellbore flow rate of tNavigator®, and N is the number of time steps.Figure 10
Comparison of pressure distribution maps over 10 days of three various simulators.
(a)
The original EDFM(b)
The novel EDFM(c)
LGR by tNavigator®Figure 11
Comparison of early-time results of wellbore flow rate.
(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 mFigure 12
Error comparison between the original EDFM and the modified EDFM.


## 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.Table 1
Input parameters of numerical case.
ParametersValueParametersValueMatrix porosity0.17Matrix permeability1 mDRock compressibility coefficient1.07×10−4MPa−1Fracture aperture10 mmFracture porosity0.4Fracture permeability20000 mDFluid viscosity2 mPa·sFluid compressibility coefficient3.02×10−3MPa−1Figure 6
Sketch of reservoir domain and a horizontal well with one fracture.Figure 7
Top view of locally refined fracture grids and their surrounding matrix grids.Figure 8
Comparison of pressure distribution maps over 1000 days of three various simulators.
(a)
The original EDFM(b)
The novel EDFM(c)
LGR by tNavigator®Figure 9
Comparison of wellbore flow rate over 1000 days of three various simulators.


## 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 Table1. 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 2×2 and 10×10meters, 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 (10×10m) 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 (2×2m and 10×10m) 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.(25)Ae=1N∑i=1Nqi−qirefqiref×100,where qi represents the wellbore flow rate calculated by EDFM, qiref represents the wellbore flow rate of tNavigator®, and N is the number of time steps.Figure 10
Comparison of pressure distribution maps over 10 days of three various simulators.
(a)
The original EDFM(b)
The novel EDFM(c)
LGR by tNavigator®Figure 11
Comparison of early-time results of wellbore flow rate.
(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 mFigure 12
Error comparison between the original EDFM and the modified EDFM.


## 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 Figure13, 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.Figure 13
Sketch of reservoir domain and a multistage fractured horizontal well.Figure 14
Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days.Figure 15
Prediction of well production performance over 1000 days.
(a)
Oil production performance(b)
Average formation pressure and reserve recovery


### 5.2. Heterogeneous Fractured Porous Media

The prior example in Figure13 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.Figure 16
Matrix permeability distribution map.Figure 17
Matrix porosity distribution map.Figure 18
Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days.Figure 19
Prediction of well production performance over 1000 days.
(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 Figure20, 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).Figure 20
Sketch of reservoir and complex fracture networks.Table 2
Input parameters of numerical case.
ParametersValueParametersValueMatrix porosity0.12Matrix permeability0.1 mDRock compressibility coefficient1.07×10−4MPa−1Fracture aperture10 mmFracture porosity0.4Fracture permeability5000 mDFluid viscosity2 mPa·sFluid compressibility coefficient3.02×10−3MPa−1Figure 21
Pressure distribution maps in the 1st, 10th, and 100th days and 1st, 3rd, and 10th years.Figure 22
Prediction of well production performance over 10 years.
(a)
Oil production performance(b)
Average formation pressure and reserve recovery


## 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 Figure13, 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.Figure 13
Sketch of reservoir domain and a multistage fractured horizontal well.Figure 14
Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days.Figure 15
Prediction of well production performance over 1000 days.
(a)
Oil production performance(b)
Average formation pressure and reserve recovery


## 5.2. Heterogeneous Fractured Porous Media

The prior example in Figure13 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.Figure 16
Matrix permeability distribution map.Figure 17
Matrix porosity distribution map.Figure 18
Pressure distribution maps in the 1st, 10th, 100th, 300th, 500th, and 1000th days.Figure 19
Prediction of well production performance over 1000 days.
(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 Figure20, 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).Figure 20
Sketch of reservoir and complex fracture networks.Table 2
Input parameters of numerical case.
ParametersValueParametersValueMatrix porosity0.12Matrix permeability0.1 mDRock compressibility coefficient1.07×10−4MPa−1Fracture aperture10 mmFracture porosity0.4Fracture permeability5000 mDFluid viscosity2 mPa·sFluid compressibility coefficient3.02×10−3MPa−1Figure 21
Pressure distribution maps in the 1st, 10th, and 100th days and 1st, 3rd, and 10th years.Figure 22
Prediction of well production performance over 10 years.
(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.



---

*Source: 9910424-2021-07-05.xml*