A modified projection-based embedded discrete fracture model (pEDFM) for practical and accurate numerical simulation of fractured reservoir

Journal of Petroleum Science and Engineering 187 (2020) 106852

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: http://www.elsevier.com/locate/petrol

A modified projection-based embedded discrete fracture model (pEDFM) for practical and accurate numerical simulation of fractured reservoir Xiang Rao a, *, Linsong Cheng a, Renyi Cao a, Pin Jia a, Hao Liu b, Xulin Du a a b

China University of Petroleum, Beijing, 102249, China Hohai University, Nanjing, 210098, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Embedded discrete fracture model Projection-based embedded discrete fracture model Complex fracture network Numerical simulation Flow barriers

pEDFM is a new numerical model which is proposed in the last three years to handle the flow simulation with complex fracture networks. Previous works show that classical EDFM has a good performance only for pro­ duction or injection simulation of multistage fractured well, but pEDFM can effectively handle more general cases than EDFM. Although pEDFM has been validated by by some works, this paper points out that there are still some limitations in previous pEDFM need to be resolved, lying on lack of a practical method to select projected faces, the proficiency in the transmissibility formula of f-m connections, and poor performances in some cases. Then, corresponding modifications are presented to resolve these limitations, and these modifications include A practical algorithm called “micro-translation method” to select projected-faces combinations, the physical transmissibility formula of additional fracture-matrix (f-m) connections and an extended pEDFM framework by considering additional fracture-fracture (f-f) connections. Some test cases are implemented to illustrate the ne­ cessity and validity of these modifications. Besides, two examples are implemented to illustrate the robustness of the modified pEDFM workflow for practical application of a fractured reservoir.

1. Introduction Accurate flow modeling and simulation in fractured reservoirs is quite challenging and important. To accurately and explicitly charac­ terize the geometric characteristics of fractures, discrete fracture model (DFM) was proposed to obtain high-quality simulation results (Snow, 1969). In DFM, fractures are handled by reducing a dimension and un­ structured matrix mesh is generated to conform to the fractures’ ge­ ometry such that fractures lie at the interfaces between matrix cells. DFM had been applied to single-phase, multiphase and compositional numerical simulation in the fracture media by combining various nu­ merical discretization methods, including finite difference method (FDM) (Slough et al., 1999), Galerkin finite element method (Noorishad and Mehran, 1982), finite volume method (FVM) (Karimi-Fard et al., 2004), control volume finite element method (Monteagudo and Fir­ oozabadi, 2004; Marcondes et al., 2010) and mixed finite element method (Alboin et al., 2002; Martin et al., 2005). However, it is much difficult to create a high-quality unstructured mesh to match arbitrary complex fracture networks. Besides, if adding or deleting some fractures or considering fractures’ dynamic propagation, the unstructured mesh needs to be re-created such that the computational efficiency is

significantly reduced. Above problems of DFM led to the emergence of embedded discrete fracture model (EDFM). Classical EDFM was proposed by Li and Lee (2008), and it was an extension of the algorithm proposed by Lee et al. (2001), in which fractures are treated in a similar way of the horizontal well in the reservoir. In EDFM, a structured grid is used to represent the matrix, and fractures are embedded in the structured grid by figuring out the connections between fracture cells and matrix cells. EDFM avoids the generation of complex unstructured conforming grids, and is convenient to be coupled with the existing numerical simulators which are based on FVM or FDM. EDFM had been successfully applied to nu­ merical simulations of various flows in fractured media (Norbeck et al., 2016; Wang and Shahvali, 2016). Panfili and Cominelli (2014) applied EDFM to simulate miscible gas injection in a fractured carbonate reservoir. Yu (2017) used EDFM to numerically simulate the shale gas transport and production with complex fractures. Dachanuwattana et al. (2018) applied EDFM and proxy-based MCMC to history match a Vaca Muerta shale oil well. Shakiba et al. (2018) used EDFM to numerically simulate flow problems in complex fracture networks calibrated by microseismic monitoring data. There are also some modifications of the numerical algorithm of EDFM, Hajibeygi et al. (2011) introduced a hi­ erarchical embedded fracture model by using iterative multiscale FVM

* Corresponding author. E-mail addresses: [email protected] (X. Rao), [email protected] (L. Cheng). https://doi.org/10.1016/j.petrol.2019.106852 Received 14 August 2019; Received in revised form 14 November 2019; Accepted 20 December 2019 Available online 23 December 2019 0920-4105/© 2019 Elsevier B.V. All rights reserved.

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Nomenclature f mi mj mk Γij Γjk Afmi Aij Ajk Apx fmi

Fracture cell The ith matrix cell The jth matrix cell The kth matrix cell Interface between mi and mj Interface between mj and mk Area of f within mi Area of Γij Area of Γjk Projection area of f on Γij along x direction

Afmi

Area of f

py Afmi

Vmi wf 〈d〉 dfmj dfmk dix djx T k

Volume of mi Aperture of f Average normal distance between a matrix cell and its containing fracture cell Distance from the center of mj to the center of f Distance from the center of mk to the center of f The half length of mi in x direction The half length of mj in x direction Transmissibility Permeability

Subscript K permeability T transmissibility

Projection area of f on Γjk along y direction

Fig. 1. Schematics of pEDFM for a fracture cell within a matrix cell, in which the blue segment denotes the fracture cell f, and the two red segments denote the xdirection and y-direction projections of f respectively, and connections are denoted by black arrows. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3. Schematics of the smaller distance criterion, in which the blue segment and the black square represent the fracture cell and a matrix cell respectively, and the red point F and O are respectively the center of the fracture cell and the center of the matrix cell. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2. The segments of a color on Γij is correspondingly the projection of the fracture cell of the same color, and the total projection of these fractures is the union of these segments on Γij which is plotted in black.

transfer flux which has a higher accuracy at early time, and extended 2D model to 3D which can effective handle inclined fractures with complex shapes in 3D reservoirs. Rao et al. (2019b) presented a modified EDFM by duplicating fracture cells to represent the cells of fracturing fluid layers, and studied the water blockage effect on water huff-n-puff pro­ cess. Rao et al. (2019c) presented a fast EDFM based on proper orthogonal decomposition (POD) method to improve computational efficiency. Rao et al. (2019d) proposed a modified EDFM to significantly improve the accuracy of early-time production simulation of multi-stage fractured well. Although classical EDFM has been validated by various studies

to improve the computational efficiency with ensuring enough accuracy. Ali Moinfar et al. (2014) developed an efficient EDFM for three-dimensional (3D) compositional reservoir simulation. Tene et al. (2016) presented a EDFM based on algebraic multiscale methods for heterogeneous flow media. Cao et al. (2019) presented an EDFM model for two-dimensional (2D) reservoir simulation, in which boundary element method (BEM) was used to handle pseudo-steady flow equation to obtain a new static transmissibility between a matrix cell and its containing fracture cells. Rao et al. (2019a) applied boundary integral equations of transient multi-phase flow PDEs to approximate the 2

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

general cases than EDFM. However, this paper points out there are still three limitations in previous pEDFM. They are: (i) There is lack of a practical algorithm to select projected faces to hold the fracture-cell projections. (ii) The transmissibility formula of additional fracture-matrix (f-m) connections is not physical. (iii) Poor performances of previous pEDFM in some cases indicate that previous pEDFM framework is not complete.

Fig. 4. Schematics of the different side criterion, in which the blue segment and the black square represent the fracture cell and a matrix cell respectively, and the red point F and O are respectively the center of the fracture cell and the center of the matrix cell.

To resolve above limitations, this paper proposes a practical algo­ rithm namely “micro-translation method” to select projected faces of fracture cells for general 3D cases, the modified physical transmissibility formula of additional f-m connections and the extended pEDFM frame­ work by considering additional fracture-fracture (f-f) connections. Some test cases are implemented to illustrate the necessity and validity of these modifications, and two examples are implemented to illustrate some 3D applications of the modified pEDFM. The modified model can effectively handle more general cases, that is, the modified pEDFM framework is more complete and practical than original pEDFM. The paper is structured as follows. In Section 2, previous pEDFM is briefly reviewed to introduce finite-volume discrete schemes of flow equations and the main differences between EDFM and pEDFM. Section 3 analyzes the limitations of previous pEDFM and makes corresponding modifications which are validated by giving test cases. Two examples are implemented to illustrated some applications of the extended pEDFM framework in Section 4. Finally, conclusions come in Section 5.

Fig. 5. Schematics of a non-trivial case, in which the black square is the matrix cell, and the blue segment is the fracture cell, and O and F are the matrix-cell center and fracture-cell center respectively. (For interpretation of the refer­ ences to color in this figure legend, the reader is referred to the Web version of this article.)

2. A brief review of pEDFM 2.1. Discrete scheme of flow equations Take the immiscible multiphase flow as an example to illuminate the discrete scheme of flow equations. The mass-conservative equations for phases are: � ��� � � � � kkra D ∂ φsa rpa ρa;sc gr þ~ qa;well ¼ (1) r⋅ ∂t Ba Ba Ba μa where k is the absolute permeability; kra , Ba , μa , pa , ρa;sc , sa are relative permeability, volume factor, viscosity, pressure, standard density and ~a;well is the flow term of phase a for wells (source saturation of phase a; q or sink term); D is the depth; t is the time. Block-center finite volume method is used to discretize the above flow equations, for the ith cell, that is, �tþΔt � �t � �� n X Vi φsa;i φsa;i qtij þ qa;well;i ¼ (2) Δt Ba;i Ba;i j¼1

Fig. 6. Schematics of a non-trivial case, in which the black square is the matrix cell, and the blue segment is the fracture cell, G is the point which is obtained by translating slightly along n1 from F. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) ⇀

described above. Tene et al. (2017) found that classical EDFM could not handle the cases in which the fracture permeability lies below that of matrix. To resolve this limitation, they proposed a projection-based EDFM (pEDFM) to extend classical EDFM adaptive to a wide range of fracture permeabilities (i.e. from flow barriers to highly conductive fractures). Besides, Jiang and Younis (2017) pointed out the other lim­ itation of classical EDFM that large errors could be induced for multi­ phase displacement across a fracture, and they made several improvements upon the classical pEDFM and applied the improved pEDFM to effectively resolve the above limitation. The improvements include a principle to select physical projected faces and modification of transmissibility formula, but the selection principle is not practical for 3D cases and the formula is not accurate enough. Ren et al. (2018) combined pEDFM and extended finite element method (XFEM) to establish a numerical model for coupling geomechanics and multiphase flow in fractured porous media. Their works (Tene et al., 2017; Jiang and Younis, 2017a,b; Ren et al., 2018) indicate that classical EDFM has a good performance only for production or injection simulation of multistage fractured well, but pEDFM can effectively handle more

where n is the number of the cells which are connected with the ith cell; The subscripts I and j denotes the quantities of corresponding physical parameters in ith and jth cells, respectively; The superscript t þ Δt and t denotes the quantities at time t þ Δt and t, respectively; Vi is the volume of the ith cell; Δt is the time interval. A two-point linear flux approximation is used to approximate the flux across a cell interface, that is, for two connected cells, � � �� � Di Dj qij ¼ Tij λa;ij pa;i pa;j ρa;sc g (3) Ba;i Ba;j Tij and λa;ij are the transmissibility and phase mobility of the connection between the ith and jth cells. λa;ij is expressed as: λo;ij ¼

3

kro;ij

μo;ij Bo;ij

;

(4)

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 7. A flowchart of the “micro-translation method”.

Fig. 8. Sketch of the second 2D non-trivial case, and in this case matrix-cell center O and fracture-cell center F coincide.

Fig. 10. Sketch of the second 3D non-trivial case, and in this case matrix-cell center O and fracture-cell center F coincide.

Fig. 11. Schematics of non-neighboring connection f-mj, in which there are two neighboring matrix cells mi and mj, and a fracture cell f in the matrix cell. Δx and Δy are the.

Fig. 9. Sketch of the first 3D non-trivial case.

4

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Table 1 Some physical properties of the reservoir and fluid. Properties

Value

Properties

Value

Matrix porosity Fracture porosity Fracture length

0.1 0.4 40 m

Oil coefficient of compressibility Oil viscosity Initial reservoir pressure

3.02 � 10 MPa 1 2 mPa s 20 MPa

Matrix permeability Fracture aperture coefficient of compressibility Water coefficient of compressibility Water viscosity Initial water saturation in matrix

1 mD 0.01 m 1.07 � 10 MPa 1 4 � 10 4 MPa 1 1 mPa s 0.2

3

4

2.2. The main difference between pEDFM and EDFM pEDFM is proposed by Tene et al. (2017) to extend classical EDFM adaptive to a wide range of fracture permeabilities (i.e. from flow bar­ riers to highly conductive fractures). Jiang and Younis (2017a,b) made several improvements upon the original pEDFM and applied the improved pEDFM to resolve the limitation that classical EDFM could induce large errors for multiphase displacement processes. Above facts indicate that, pEDFM can be applied to more general cases, unlike classical EDFM, it may be only adaptive to the production or injection

Fig. 12. Schematic of the regions division.

the upstream scheme is used for the terms (relative permeability) subject to the saturation and the arithmetic average scheme is used for the terms (viscosity and volume factor) subject to pressure, that is,

μo;ij ¼

μo;i þ μo;j 2

� 8 > > ​ ​ ​ ​ ​ ​ if ​ po;i k > ro;i <

; Bo;ij ¼

Bo;i þ Bo;j ; kro;ij ¼ � > 2 > > : kro;j ​ ​ ​ ​ ​ ​ if ​ po;i

ρa;sc g

Di Ba;i

ρa;sc g

Di Ba;i

�

� �

�

Ti Tj Ti þ Tj

ρa;sc g

Dj Ba;j

po;j

ρa;sc g

Dj Ba;j

� <

Tij is generally the half of the harmonic mean of two halftransmissibility, that is Tij ¼

po;j

� (5)

�

simulation for multistage fractured horizontal well. As shown in Fig. 1(a), there is a fracture cell f in the matrix cell mi, and mi has two neighboring matrix cells mj and mk. The area of the fracture cell is denoted by Afmi ., in classical EDFM, there is only one f-m connection (denoted by the black arrow) for this case, i.e., f-mi, and the transmissibility of this connection is � � � 1 � km 2Afmi kf 2Afmi ; Tf ¼ Tfmi ¼ T mi1 þ T f 1 ; Tmi ¼ i (8) 〈d〉 wf

(6)

so for the ith cell, the discrete scheme of the flow equation is obtained as: � � ���t n � X � Di Dj ρa;sc g þ qa;well;i Tij λa;ij pa;i pa;j Ba;i Ba;j j¼1 �tþΔt � �t � �� Vi φsa;i φsa;i ¼ (7) Δt Ba;i Ba;i

However, in pEDFM, as shown in Fig. 1(b), the fracture cell f is projected to the matrix-cell interfaces Γij and Γjk along x and y directions py with the corresponding projection areas denoted by Apx fmi and Afmi

respectively, and additional f-m connections including f-mj, f-mk are constructed besides f-mi. For f- mi connection, the transmissibility is calculated as:

Fig. 13. Schematic of the uniformly distributed points which are denoted by gray dots.

Fig. 14. Schematics of the reservoir model. 5

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 15. Oil saturation profiles at 30 days of the 4th case.

� 1 � km Afm kf Afmi Tfmj ¼ T mj1 þ T f 1 ; Tmi ¼ i i ; Tf ¼ 〈d〉 wf

where dfmj is the distance from the center of mj to the center of f. Then Tfmj can be calculated as

(9)

where kmi and kf are the permeabilities of mi and f respectively, wf is the aperture of f. 〈d〉 is the average normal distance between f and the point within mi, which is approximated by assuming that pressure varies lin­ early in the direction normal to each fracture cell in a matrix cell (Li and Lee, 2008). R !! j n ⋅ r jdV 〈d〉 ¼ mi (10) V mi

� � Tfmj ¼ T mj1 þ T f 1

of only Tm1j and T f 1 . Therefore, they presented a new formula for calculating Tfmj by adding Tmi into account, that is, � � Tfmj ¼ T mj1 þ T mi1 þ T f 1

measure to determine the transmissibility of the additional f-m con­ nections, the formula of which is illuminated as follows by taking f-mj connection as an example. For f-mj connection, Tene et al. (2017) proposed that the ‘matrix half’ of Tfmj is: kmj Apx fmi dfmj

(12)

Jiang and Younis (2017) discussed that the effect of mi should be considered in Tfmj . They gave a case that assuming the permeability of mi is ultra-low, and the fluid flows from the left to the right. In fact, there will be very little fluid from mj to f because of the ultra-low permeability of mi, however, lots of fluid more than actual case will be transferred from mj to f if Tfmj is simply approximated by using the harmonic average

where ! n is the unit normal vector on the surface of f, ! r is the vector starting at the center of f and ending at a point in mi. py The areas of the fracture-cell projections (i.e., Apx fmi and Afmi ) are the

T mj ¼

1

1

(13)

Similarly, the transmissibility of f-mk connection can be given as: � 1 � kmk Apy fmi Tfmj ¼ T mk1 þ T mi1 þ T f 1 T mk ¼ dfmk

(14)

where dfmk is the distance from the center of mk to the center of f. Because adding additional f-m connections results in that some part

(11)

6

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 16. Oil saturation profiles at 60 days of this case when mesh sizes are 1 m and 2 m

7

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 19. Schematics of the specific example, in which two fracture cells f1 and f2 are within matrix cells m1 and m2 respectively, and Γ12 is the interface of m1 and m2.

Tix ¼

� kmi Aij dix

Apx f

� ; Tjx ¼

� kmj Aij

Apx f

�

djx

(15)

where Aij is the area of the interface Γ ij between mi and mj. It can be seen that the effective flow area of mi-mj connection is equal to the area of Γ ij minus the total projection area of f. Jiang and Younis (2017a,b) pointed that, if there are multiple fracture cells projected onto a same interface, the total projection area should be the area of the union of these pro­ jections. As shown in Fig. 2, there are ni and nj fracture cells within mi and mj projected onto Γij , and the total projection area Apij on Γ ij is given

Fig. 17. Average error plots for results from EDFM, previous pEDFM and current pEDFM.

as:

� nj py i Apij ¼ S Di [ Dj ; Di ¼ [nl¼1 Dpx fl mi ; Dj ¼ [l¼1 Dfl mj

(16)

where Dpx fl mi is the projection on Γ ij of the fracture cell fl which is within py

mi, Dfl mj is the projection on Γjk of the fracture cell fk which is within mj,

Sð⋅Þ is a function whose independent variable is a region and the dependent variable is the area of the region. 3. Analysis of existing limitations and presentation of corresponding solutions Although the validity of pEDFM for single-phase and multiphase flow simulation has been illustrated by Tene et al. (2017), Jiang and Younis (2017a,b), there are still three major limitations of previous pEDFM, which will be analyzed and resolved in detail as follows.

Fig. 18. Average error plots for oil saturation profiles of 60 days with different mesh sizes.

3.1. The first limitation and its solution 3.1.1. Lack of a practical method to select projected faces Tene et al. (2017) did not explicitly discuss about the selection cri­ terion, but they indicated to select the faces that have smaller distance from the center of the fracture cell along each dimension, and this se­ lection criterion could be simply denoted by “the smaller distance cri­ terion”. As shown Fig. 3, if the fracture-cell center F is closer to face 2 in x direction, then in x direction face 2 is selected to hold the projection. Otherwise, face 3 is selected. If fracture-cell center F is closer to face 1 in y direction, then in y direction fracture cell is projected on face 1. Otherwise, it is projected on face 4. According to the above criterion, the fracture cell in Fig. 3 is projected on face 2 in x direction and on face 4 in y direction. Later, Jiang and Younis (2017) proposed that, when the distances between the fracture-cell center and two neighboring faces along one direction are the same, the smaller distance criterion in this direction will be invalid. If the projected faces are selected randomly at this time, it may result in unphysical fracture-matrix connections. As shown in Fig. 4, the distances from fracture-cell center F to face 1 and face 4 in y direction are the same. If face 3 and face 4 are respectively selected as the projected faces in x and y direction, then it will not conform to the physical impact of fracture geometry on fluid flow. To resolve above similar cases, Jiang and Younis added a selection criterion for 2D case which can be simply denoted as “the different side criterion”, that is, the

Table 2 The tabular data which was used for Fig. 17 30 days EDFM Previous pEDFM Current pEDFM

2.6708*10 1.0113*10 6.9030*10

60 days 4

9.4588*10 3.4991*10 1.7945*10

5 6

4 5 5

Table 3 The tabular data which was used for Fig. 18. 1m EDFM Previous pEDFM Current pEDFM

3.2271*10 7.6625*10 3.1192*10

2m 4 6 6

9.4588*10 3.4991*10 1.7945*10

4m 4 5 5

2.0455*10 1.5991*10 8.6610*10

3 4 5

of the flow areas of m-m connections may be blocked, the trans­ missibility of corresponding m-m connections should be modified. Take mi-mj connection as an example, the half-transmissibility Tix and Tjx for Tmi mj are given as:

8

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 20. Schematics of the transmissibilities of the specific example.

direction. As shown in Fig. 5, take a 2D non-trivial case as an example, a fracture cell (denoted by f) is contained in a matrix cell (denoted by m).

n1 and n2 are two opposite normal vectors to f, and the vector starting at the matrix-cell center (denote as O) and ends at fracture cell center ⇀

⇀

⇀

(denote as F) is denoted as OF . Then, to select the normal vector from n1 ⇀

⇀

and n2 which meets n ⋅OF � 0, that is, the angle between the selected ⇀

⇀

⇀

normal vector and OF is less or equal to 90� . If both n1 and n2 meet the requirement, both vectors can be selected. In the example shown in ⇀

⇀

Fig. 5, it can be seen that n1 meets the requirement. Step 3: As shown in Fig. 6, translate slightly the fracture-cell center F ⇀

Fig. 21. Sketch of the case with multiple fracture cells in two neighboring matrix cells, in which the segments of a color on Γij is correspondingly the projection of the fracture cell of the same color.

to point G along the selected normal vector n1 . Assuming the co­ ⇀

⇀ n1

ordinates of F is ðx; yÞ and is denoted in components as ðnx ; ny Þ, then coordinates of G is ðx þ lnx ;y þ lny Þ. “slightly” means that l is a very tiny number which can be recognized by computer, e.g., 0.0001 (the reason of a tiny translation is to make same distances along some directions become not the same, without damaging the relative distance relation that the distances between fracture-cell center to two neighboring faces along other directions are already not equal). Step 4: Based on the position of G, apply “the smaller distance cri­ terion” to select the projected faces along each direction. If there re­ mains the situation that the distances from point G to the two neighboring faces along a direction are the same, then any faces in this direction can be selected. In the above example, it can be seen that, because face 1 and face 3 are nearer to G along x and y directions respectively, face 1 and 3 should be selected as the projected faces.

two projection lines are not allowed to intersect with the fracture line on the same side of the fracture-cell center. According to the added crite­ rion, Fig. 4 shows that, only when face 1 and 3 are selected, the two intersections A and B will be on the different sides of F such that unphysical connections are avoided. However, “the different side criterion” is not practical especially for 3-D cases. This is because, in general 3D cases, the fracture cell maybe arbitrary-shaped convex polygons. The calculation and comparison of the intersecting lines between the projection planes of fracture cells in x, y, and z direction and the fracture-cell plane are complex, and there will be lots of iterations of face-selections in computing process. All these greatly increase complexity of the algorithm. Therefore, to make pEDFM adaptive to 3D general cases, a practical and simple face-selection al­ gorithm needs to be developed. 3.1.2. A practical algorithm to select projected-faces combinations In this section, a simple and practical algorithm for both 2D and 3D cases is presented, which is simply denoted as ‘micro-translation method’, and the detailed procedures including four steps are as follows. Step 1: If the distances from the fracture-cell center to the two neighboring faces along each direction are not the same in all directions, then the projected face in each direction can be determined by “the smaller distance criterion” without any subsequent steps. This is a trivial case. Step 2: For the non-trivial cases, the distances from fracture-cell center to the two neighboring faces are the same along at least one

Fig. 23. Schematics of reservoir model in this case.

Fig. 22. Schematics of the connections and corresponding transmissibility of the specific example from modified model. 9

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 24. Oil saturation profiles at 60 days when fracture 2 is a flow barrier.

10

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 25. Oil saturation profiles at 60 days when the permeability of fracture 2 is 1mD.

11

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 26. Oil saturation profiles at 60 days when the permeability of fracture 2 is 20000mD.

Fig. 27. Average error plots of results from EDFM, previous pEDFM and cur­ rent pEDFM.

Fig. 28. Average error plots at different mesh size when the permeability of fracture 2 is 1mD.

12

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Table 4 The tabular data which was used for Fig. 20 0 mD EDFM Previous pEDFM Current pEDFM

1 mD

0.0087 0.0118 1.0789*10

20000 mD

0.0015 1.1352*10 1.3396*10

4

0.0016 5.0755*10 8.6167*10

4 5

5 6

Table 5 The tabular data which was used for Fig. 21 1m EDFM Previous pEDFM Current pEDFM

2m

9.5176*10 3.1683*10 3.3647*10

4

4m

0.0015 1.1352*10 1.3396*10

5 6

0.0030 0.0011 3.2918*10

4 5

4

Fig. 29. Schematics of the reservoir model in example 1.

A simple flowchart of the presented above “micro-translation method” is shown in Fig. 7. 3.1.3. Test cases A non-trivial 2D and two 3D cases are given to illustrate the validity of the above algorithm, and these cases are clearly shown as follows. 3.1.3.1. Case 1. As shown in Fig. 8, it can be seen that the distances from F to the two neighboring faces in x and y directions are both the ⇀

⇀

same, which is a much special case. Under this case OF ¼ 0 and ⇀

Fig. 30. Oil-phase pressure profiles.

⇀ both n1

from F to the two neighboring faces in x, y, z directions are the same

and n2 meet n ⋅OF � 0, thus both n1 and n2 can be selected. Assume n1 ⇀

⇀

⇀

⇀

⇀

⇀

is selected and translate F slightly to G along , then by “the smaller distance criterion”, face 3 and 1 can be selected as the projected faces in

⇀

selected, translate F slightly to G along . It can be seen that G is closer to face 4 in x direction, closer to face 6 in z direction, but still the same close as to face 2 and 5 in y direction, thus any one of face 2 and 5 can be selected. That is, the projected faces combination can be selected as faces 4-2-6 or faces 4-5-6, which are both physical configurations.

Similarly, if n2 is selected then the projected faces combination can be selected as faces 4-2-1 or faces 3-5-1, which are also physical configurations. From the results, we can see that the proposed algorithm ‘microtranslation method’ can indeed determine the physical projected-faces combinations, and this method requires only simple algebraic opera­ tions, which can be efficiently handled by computers. ⇀

smaller. Firstly calculate the normal vectors n ¼ � α1 � α2 , and ⇀

⇀

because n1 meets n1 ⋅OF � 0 it can be selected as the needed normal ⇀

⇀

⇀

3.1.3.2. Case 2. As shown in Fig. 9, the distances from fracture-cell center F to face 1 and 6 are equal in z direction, and compared to face 3 and 5 the distances from F to face 4 and 2 in x and y direction are ⇀

⇀

⇀ n1

⇀

⇀

⇀

thus any of these two opposite normal vectors can be selected. If n1 is

x and y direction. Similarly, if n2 is selected face 2 and 4 will be selected as projected faces. It can be found that both above two selections are physical configurations, and proposed algorithm avoids the unphysical configurations, such as face 2-1 or face 3–4.

⇀

⇀

respectively. Under this case OF ¼ 0 and both n1 n2 satisfy n ⋅OF � 0,

⇀ n1

vector. Then translate F slightly to G along n1 , and it can be seen that the distance from G to face 6 is smaller in z direction. According to “the smaller distance criterion”, face 6, 2, and 4 are selected as the projected faces in x, y and z direction respectively, and the selection is indeed a physical configuration. ⇀

3.2. The second limitation and its solution 3.2.1. The proficiency in the transmissibility formula of f-m connections As shown in Fig. 11, there is a non-neighboring f-mj connection be­ tween the fracture cell f and the matrix cell mj in pEDFM, and we denote the projection area of f on interface Γij as Apx fmi ¼ Afmi . From Eq. (13), the

3.1.3.3. Case 3. Fig. 10 shows a more special situation, where matriccell center O coincides with fracture-cell center F, and the distances 13

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

transmissibility of the f-mj connection is calculated as: � � 1 Tfmj ¼ T mj1 þ T mi1 þ T f 1 where Tmj ¼

kmj Afmi dfmj ,

Tmi ¼

kmi Afmi 〈d〉 ,

Tf ¼

(17)

kf Afmi wf .

When the location of f gradually approaches the interface Γij , the

value of Tm1j in Eq. (17) should be able to decrease continuously to 0.

However, from Eq. (17), when f andΓ ij are coincident there will be T m1j ¼ Δx 2kmi Afmi

6¼ 0, and it shows that the expression of Tfmj in Eq. (17) is not

physical. We give some specific calculations to clearly illuminate this issue as follows: As shown in Fig. 11, suppose the dimensions of matrix cell in all directions are 2 m, and the fracture aperture is 0.01 m, and the fracturecell and matrix-cell permeability are kf ¼ 20D and kmi ¼ kmj ¼ 1mD respectively. When f and Γ ij are coincident, Tmj ¼ Tmi ¼ 4mD⋅m can be obtained, and because the fracture width is small, Tf ≫Tmj ¼ Tmi . Thus, Tfmj ¼ 4

1

þ4

1

þ Tf

1

�

1

� 2 ​ mD⋅m

(18)

But when f andΓ ij are coincident, Tfmj is physically obtained as � 1 Tfmj ¼ 4 1 þ Tf 1 � 4 ​ mD⋅m (19) It can be seen that the calculated transmissibility is only half of the physical result, and this error may induce errors. So it is necessary to construct a physical transmissibility formula for the additional f-m connections in pEDFM. 3.2.2. The modification of the transmissibility of f-m connections Section 3.2.1 indicates that the previous studies on pEDFM were somewhat crude for the calculation formulas of transmissibility of the additional connections. Therefore, in this paper we redefines Tmi 、Tmj and Tmk in Eq. (17) to solve this problem. As shown in Fig. 12, taking a two-dimensional mesh as an example, suppose ! n ¼ ðnx ; ny Þ is the unit normal vector obtained by the algorithm ‘micro-translation method’ in Section 3.1.2, and r ¼ ðx; yÞ is the vector whose start point is the fracture-cell center F and end point is a point P in the matric cell mi. Divide matrix cell mi into two regions Ω1 and Ω2 which are defined as follows: ⇀

Fig. 31. Oil saturation profiles.

Ω ¼ Ω1 [ Ω2 ⇀ ⇀ where fΩ1 j! n ⋅ r � 0g, fΩ2 j! n ⋅ r < 0g.

Fig. 32. Schematics of the numerical model of a well group in Changqing oil field.

14

(20)

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 33. The relative permeability and capillary pressure curves in matrix and FFL. 1x connection f-mj from Eq. (13) should be corrected toT Ω mi :

Table 6 Some physical parameters of the reservoir and fluid. Properties

Value

Properties

Value

Matrix porosity

0.13

0.6 mD

Matrix permeability in y direction Fracture porosity Fracture aperture

0.2 mD

Matrix permeability in x direction Matrix permeability in z direction Fracture permeability Rock compressibility coefficient water compressibility coefficient

1.311

1.00001

Oil volume factor at bubble point pressure Reservoir size

2 mPa s 16 MPa

Water viscosity Initial water saturation

2060mⅹ1580 mⅹ30 m 0.6 mPa s 0.47

820 kg/m3

Standard water density

1000 kg/m3

Oil compressibility coefficient Bubble point pressure Water volume factor at 20 MPa Oil viscosity Initial reservoir pressure Standard oil density

0.40 10 mm 3.02 � 10 3 Mpa 1 7.0

T Ωmi1 x ¼

R

⇀⇀

Ω1

jn⋅ r jdV

(21)

VΩ1

From Eq. (21), when VΩ1 →0, 〈dΩ1 〉→0,TΩm1i → þ ∞, and

0.06 mD 20000 mD 1.07 � 10 4 Mpa 1 5 � 10 4 Mpa

kmi Apx f ; 〈dΩ1 〉 ¼ 〈dΩ1 〉

� Tfmj ¼ T mj1 þ T Ωmi1 x

� 1

þ Tf

1

1

� � → T mj1 þ T f 1

1

(22)

It can be seen that, when the fracture cell gradually approaches to the interface, the transmissibility of f-mj connection naturally degenerates to that in the case of discrete fracture models. When calculating 〈dΩ1 〉, set A can be used to convert the integral operation in Eq. (21) into simple algebraic operations, that is: P ⇀⇀ P ⇀⇀ R ⇀⇀ ½n⋅ r ðpÞðΔxΔyÞ=ðmnÞ� ½n⋅ r ðpÞ� n ⋅ r dV p2A p2A Ω 〈dΩ1 〉 ¼ 1 ¼ (23) � VΩ 1 jAjðΔxΔyÞ=ðmnÞ jAj

1

where p is a point in set A, n⋅ r ðpÞ is the value of n⋅ r at p. Similarly, the Tmi in the transmissibility formula of f-mk connection 1y should corrected to T Ω mi : ⇀⇀

The approximate characterization of regions Ω1 , Ω2 can be done simply, and the specific process is given as follows: as shown in Fig. 13, mn points (m points per row, and n points per column) are uniformly distributed in the x and y direction respectively such that the control area of each point isðΔxΔyÞ =ðmnÞ, then determine if each point does not ⇀ satisfy with the condition of ! n ⋅ r � 0. Suppose the points satisfying ! n⋅

T Ωmi1 y ¼

⇀⇀

kmi Apy f 〈dΩ1 〉

(24)

Then the transmissibility Tfmk is given as: � Tfmk ¼ T mk1 þ T Ωmi1 y

r � 0 construct the set A, and the number of elements in set A is denoted ⇀ asjAj. The points that are not satisfied ! n ⋅ r � 0 construct set B, and the

⇀

� 1

þ Tf

1

1

(25)

Region Ω2 can be regarded as the fluid flow region for the nonneighboring connection f-mi, thus Tmi in the transmissibility formula of connection f-mi is redefined as TΩm2i :

number of elements in set B is denoted asjBj. It is obvious that jAj þ jBj ¼ mn, and the area (volume) of Ω1 and Ω2 are jAjðΔxΔyÞ =ðmnÞ and jBjðΔxΔyÞ =ðmnÞ respectively. Thus, region Ω1 can be characterized by set A and the control area of the points in set A, the region Ω2 can be characterized by set B and the control area of the points in set B. As shown in Fig. 12, if fluid flows from mj to f, the regions’ order through which the fluid physically flows is mj Ω1 f. In this case pEDFM directly constructs the non-neighboring connection f-mj, and the transmissibility of the connection mi-mj is reduced to 0, that is, all the fluid is assumed to flows directly from mj to the fracture element f. Therefore, it can be known that the closer the fracture cell is to the projected interface, the flow model from pEDFM is more in line with the actual physical process, and the accuracy will be higher. From the above analysis, the f-mj connection should be the connection between a series of the mediummj , Ω1 and f, so Tmi in the transmissibility formula of

R T Ωmi2 ¼

k mi A f ; 〈dΩ2 〉 ¼ 〈dΩ2 〉

⇀⇀ Ω2

n⋅ r dV

(26)

VΩ2

〈dΩ2 〉 can also be approximated by simple algebraic operations: P ⇀ R ⇀⇀ ½n⋅ r ðpÞ� n ⋅ r dV p2B Ω 〈dΩ2 〉 ¼ 2 � jBj VΩ 2

(27)

where p is a point in set B, n⋅ r ðpÞ is the value of n⋅ r at p. From above, It can be seen that the main idea of resolving problem Ω1 y Ω2 1x (iii) is to refine initial Tmi to T Ω mi , T mi and T mi by dividing matrix cell mi into two regions Ω1 and Ω2 , region Ω1 can be regarded as the region in ⇀⇀

15

⇀⇀

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 34. NTG profiles of three layers.

which f interacts with mj and mk, and the remaining region Ω2 is the region where f interacts with mi.

profiles from different methods at 30 days is shown in Fig. 15, and to study the performances of these models with different mesh sizes, the oil saturation profiles at 60 days when mesh sizes are 1 m and 2 m are compared in Fig. 16. The average error plots of classical EDFM, previous pEDFM and current pEDFM at 30 days and 60 days are shown in Fig. 17. From the results we can see that the saturation profile from current pEDFM after the modifications of transmissibilities of additional f-m connections has a more satisfactory match to that from ECLIPSE and has a less average error. Besides, from Figs. 15 and 16, it can be seen that EDFM induces obvious errors in oil saturation profiles in which the injected water converges to the right end of the fracture, but the right profile should be that the injected water moves forward almost evenly along the fracture. The error plots of these models with different sizes which are shown in Fig. 18 illustrate that current pEDFM has a high accuracy and good convergence. Compared to EDFM and previous pEDFM, and these indicate that the modification is effective. The tabular data for Figs. 17 and 18 are listed in Table 2 and Table 3.

3.2.3. A test case The pEDFM after the above improvement is denoted as current pEDFM, and pEDFM presented by Tene et al. (2017), Jiang and Younis (2017a,b) is named as previous pEDFM. A test case is given in which embedded fracture is vertical and aligned with the coordinate axes, so local grid refinement (LGR) method can handle this case and the results from LGR can be regarded as the reference solution. In the test case, we illustrate the validity of the improvement by comparing the perfor­ mances of classical EDFM, previous pEDFM, current pEDFM and LGR. To measure the error, we define the average error as follows. Ae ¼

N 1 X ðeÞ So i N i¼1

ðnÞ �2

So i

(28)

where the superscript (e) denotes the calculated oil saturation value from classical EDFM, previous pEDFM, or current pEDFM, and the su­ perscript (n) denotes the reference solution from LGR, and N is the grid number. In this case, as shown in Fig. 14, there is a fracture of permeability 20000mD in the reservoir, and an injector and a producer penetrate the whole reservoir thickness and locate at the left lower corner and right upper corner respectively. The injector is at a constant injection rate 1 m3/d, and the producer is at a constant bottom hole pressure 10 MPa. Other basic physical parameters are shown in Table 1. Oil saturation

3.3. The third limitation and its solution 3.3.1. Poor performances of previous pEDFM in some cases Compared with classical EDFM, new f-m connections are added and original m-m connections are simultaneously weakened in previous pEDFM. However, only adding new f-m connections is not enough, because in some cases the results of previous pEDFM is not physical. A specific example is taken to illustrate this critical issue. As shown in Fig. 19, there are fracture cells f1 and f2 within matrix

16

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 35. Saturation profiles of the upper layer at 1 year, 10 years, 30 years and 50 years.

cells m1 and m2 respectively, and it can be known that f1 and f2 are projected on the common interface Γ12 in x direction, and the corre­ px sponding projection areas are denoted as Apx fi and Afj :

of important non-neighboring connections, that is, the connections be­ tween fracture cells within different matrix cells but projected on the same interface. As shown in Fig. 21. Assume there are ni fracture cells in i the matrix cell mi (denoted as f m r ; r ¼ 1; 2⋯ni ) projected onΓ ij , the corresponding projection and projection area are denoted asDpfmi ðr ¼ 1; 2⋯ni Þ and SðDpfmi Þðr ¼ 1; 2⋯ni Þ respectively. There are nj

(29)

px Apx fi ¼ Afj ¼ Aij

From Eq. (15), it can be known that Tjx ¼ 0, Tix ¼ 0, so Tmi mk ¼ 0. Suppose kfj ≫0, kfi ¼ 0(i.e. fracture cell fj is a flow barrier), kmi 6¼ 0, kmj 6¼ 0, so Tfj ≫0, Tfi ¼ 0. From Eq. (13) we obtain that: � � Tfi mj ¼ T mj1 þ T mi1 þ T fi 1

1

¼0

(30)

� 1 � Tfj mi ¼ T mi1 þ T mj1 þ T fj 1 ≫0

(31)

r

r

j fracture elements in the matrix cell mj (denoted as f m s ; s ¼ 1; 2⋯nj ) projected on Γij , the corresponding projection and projection area are denoted as Dpmj ðs ¼ 1; 2⋯ni Þ andSðDpmj Þðs ¼ 1; 2⋯ni Þ.

fs

fs

Compared with previous pEDFM, the additional non-neighboring mj i connections between f m r in mi and f s in mj are added, and the flow area between these two fracture cells is measured by the area of the p p i j intersection of the projections of f m and f m r s , that is, SðDf mi \D mj Þ.

In Fig. 20, the connection whose transmissibility is not 0 is marked with a blue arrow, and the connection whose transmissibility is 0 is marked with a red arrow. Suppose that there is fluid flowing from the left to the right, because f2 is a flow barrier, the fluid cannot flow from the left side to m3. However, due to Tf1 m2 ≫0, Tm2 m3 > 0, there is still fluid unphysically flowing from the left side to m3 along the route of m1-f1-m2m3 indicated by the blue arrow in Fig. 20. Therefore, it is necessary to made a modification of previous pEDFM to obtain physical results for general cases.

r

fs

Because these two fracture cells are respectively in the matrix cells mi and mj, the modification idea in Section 3.2.2 is combined to define the mj i transmissibility of f m r -f s connection as: � � 1 Tf mi f mj ¼ T f m1j þ T Ω1 xmi mj 1 þ T Ω1 xmi mj 1 þ T f m1i (32) r

s

s

where TΩ1 xmi

mj mi ;f r f s

mj ;f r f s

¼

mi ;f r f s

kmi SðDpmi \Dpmj Þ fr

Ω 〈di 1 〉

fs

, TΩ1 xmi

mj mj ;f r f s

r

¼

kmj SðDpmi \Dpmj Þ fr

〈dj 1 〉 Ω

fs

i fm r

.

For the non-neighboring connection between and mj, because a portion of its flow area is occupied by the flow between the connections mi i j between f m and f m r s ; s ¼ 1; 2⋯nj , the effective flow area of f r -mj

3.3.2. Extended pEDFM framework: the additional f-f connections Section 3.3.1 indicates that previous pEDFM neglected another class

17

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 36. Pressure profiles of the upper layer at 1 year, 10 years, 30 years and 50 years. n

between the added f-f connection and f-m connections, and the sum of these occupied flow areas equals to the area of the union of all fracture cells’ projections, thus the effective flow area of connection mi-mj is the nj i Dpfmi Þ [ð[s¼1 Dpmj ÞÞ. Thus, the area of interface SðΓ ij Þ minus Sðð[nr¼1

j connection equals to SðDpfmi Þ minus Sð[s¼1 ðDpfmi \Dpmj ÞÞ. Then we combine r

r

fs

the idea in Section 3.3 to define the transmissibility Gf mr i mj as: � � Tf mr i mj ¼ T f m1i þ T Ωm1;fxmi m 1 þ T mj1 r

T Ωm1;fxmi m i r j

T mj ¼

¼

i

1

(33)

j

r

h � � kmi S Dpfmi

� � ��i nj S [s¼1 Dpfmi \ Dpmj

r

〈dΩi 1 〉

h � � kmj S Dpfmi

(35)

fs

r

djx

where the permeability of

i fm r

Tix ¼

is kf mr i , the permeability of mi is kmi , the

permeability of mj is kmj . j m Similarly, for f m s -mi connection, the transmissibility Tf j m is defined s

as:

� Tf mj mi ¼ T f m1j þ T Ω1 xmj s

T mj ¼

1

mj ;f s mi

s

T px mi ¼

�

h � � kmi S Dpmj fs

1

� � ��i nj S [s¼1 Apfmi \ Apmj fs

r

〈dΩj 1 〉 h � � kmj S Dpmj fs

Tjx ¼

i

� � ��i nj S [s¼1 Apfmi \ Apmj fs

r

djx

(39)

� � kmi S Γij

S

� � ��� �� nj i [nr¼1 Dpfmi [ [s¼1 Dpmj

(40)

fs

r

dix � � kmj S Γij

S

� � ��� �� nj i [nr¼1 Dpfmi [ [s¼1 Dpmj

(41)

fs

r

djx

As can be seen from the above definitions, when ni ¼ 0 or nj ¼ 0, that is, only one of mi and mj contains fracture cells projected on the interface Γ ij , the above transmissibility formulas of various connections naturally degenerate to those of previous pEDFM. Taking nj ¼ 0 as an example to illuminate it: nj When nj ¼ 0,SðDpmj Þ ¼ 0ðr ¼ 1; 2⋯nj Þ,Sð[s¼1 Dpmj Þ ¼

(36)

þ T mi1

1

where

� � ��i nj S [s¼1 Dpfmi \ Dpmj

r

� � Tmi mj ¼ T ix1 þ T jx1

(34)

fs

r

fs

r

transmissibility Tmi mj is defined as:

(37)

fs

fs

n

j i 0SðDpfmi \Dpmj Þ ¼ 0Sð[nr¼1 ðDpfmi \Dpmj ÞÞ ¼ 0, Sð[s¼1 ðDpfmi \Dpmj ÞÞ ¼ 0, r

fs

Then we obtain that

(38)

Tf mj mi ¼ 0; Tf mi f mj ¼ 0 s

j m where the permeability of f m s is kf j . s

For mi-mj connection, since the flow area is occupied by the fluid flow 18

r

s

r

fs

r

fs

(42)

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 37. Saturation profiles of the middle layer at 1 year, 10 years, 30 years and 50 years.

Ti ¼

� � kmi S Γij

� �� nj S [s¼1 Dpmj fs

dix

; Tj ¼

� � kmj S Γ ij

� �� nj S [s¼1 Dpmj

(the permeability is zero), 1mD fracture (the permeability is same as matrix permeability) or 20000mD fracture (the permeability is same as fracture 1) respectively. The dimensions of every matrix cell are all 2 m in x, y and z directions. A water-injection well is at the lower-left quarter with a constant injection rate of 1 m3/d, and a production well is at the upper-right corner with a constant bottom hole pressure of 10 MPa. Other physical parameters are the same with those in Table 1. The comparisons of oil saturation profiles from different methods at 60 days with different permeabilities of fracture 2 are shown in Fig. 24, Fig. 25

(43)

fs

djx

Above formulas are consistent with Eqs. (14) and (15), thus previous pEDFM is a special case of the model presented here. Next investigate whether the above modified model can effectively solve the limitation stated in section 3.3.1. In that example, n1 ¼ 1, n2 ¼ 1, Apx ¼ Apxmj ¼ Aij , then, m f i f1

1

Aij Apmj fs

� � � � nj i [nr¼1 Apfmi [ [s¼1 Apmj ¼ 0; Ti ¼ Tj ¼ 0; Tij ¼ 0; Apfmi f r r � � s nj Apmj \ Apfmi ¼ 0; Tf mr i mj ¼ 0; Tf mj mi ¼ 0 [s¼1 fr

s

� � nj Apfmi \ Apmj ¼ 0; [s¼1 r

fs

and Fig. 26 respectively. Moreover, the average error plots of classical EDFM, previous pEDFM and current pEDFM are shown in Fig. 27. To study the performances with different mesh sizes, different mesh sizes (1 m, 2 m and 4 m) are adopted to run the case when the permeability of fracture 2 is 1mD, and oil saturation profiles at 60 days when mesh sizes are 1 m and 2 m are illustrated in Fig. 25, and the error plots with different mesh sizes are compared in Fig. 28. The tabular data for Figs. 27 and 28 are listed in Table 4 and Table 5. From the comparison results, we can see that, compared with clas­ sical EDFM and previous pEDFM, current pEDFM has a higher accuracy, especially when fracture 2 is a flow barrier. Fig. 24 shows that the flow barrier does not function well in the profiles from EDFM and previous

Because kf mj ¼ 0, Tf mr i ¼ 0, and it can be obtained that Tf mi f mj ¼ 0. As 1

(44)

s

r

s

shown in Fig. 22, the connection whose transmissibility is not 0 is marked with a blue arrow, and the connection whose transmissibility is 0 is marked with a red arrow. The above modified model converts Figs. 20–22. It can be seen that fluid cannot flow from the matrix cell m1 to the right side, and f2 can physically play the role of a flow barrier. 3.3.3. A test case As shown in Fig. 23, the permeability of matrix and fracture 1 are 1mD and 20000mD respectively, and fracture 2 is set as a flow barrier 19

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 38. Pressure profiles of the middle layer at 1 year, 10 years, 30 years and 50 years.

20

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 39. Saturation profiles of the lower layer at 1 year, 10 years, 30 years and 50.

pEDFM, that is the fluid still flows across the flow barrier, but a much satisfactory match of the oil saturation profiles between proposed improved pEDFM and LGR. From Figs. 24 and 25, the water displace­ ment edge from current pEDFM is wider along x direction and shorter along y direction than that from previous pEDFM, which is more consistent with the result from ECLIPSE. The average error plots also quantitatively and intuitively indicate current pEDFM makes a signifi­ cant step-up in computational accuracy. Besides, Fig. 21 vividly illus­ trate that the accuracy will increase rapidly when the mesh size decreases, that is, current pEDFM has a good convergence. Therefore, this case illustrates that additional f-f connections indeed should be added in pEDFM, that is the current pEDFM framework is the most complete up to now.

500 days are shown in Fig. 30 and Fig. 31 respectively. 4.2. Example 2 This example is to illustrate the robustness of the proposed EDFM workflow for practical application to an improved oil recovery (IOR) study for a fractured reservoir. As shown in Fig. 32, the studied reservoir model is a well group in Changqing oil field of China which includes three multi-stage fractured horizontal well (producer) and eight straight well (injectors). The heights and length of these hydraulic fractures are generally different, and asymmetry relative to horizontal wellbore. The initial reservoir pressure and initial water saturation are 10 MPa and 0.47. Producers are at a constant bottom hole pressure (BHP) 8 MPa, and the injectors are at a constant BHP 35 MPa. The relative permeability and capillary pressure curves are plotted in Fig. 33, and other physical parameters of the well group are summarized in Table 6. The ratios of effective layer thickness to layer thickness (i. e., the same conception as NTG in commercial software ECLIPSE) are plotted in Fig. 34. Then the oil saturation and oil-phase pressure profiles of three layers at 1 year, 10 years, 30 years and 50 years can be obtained in Figs. 35–40 via the proposed pEDFM model. It can be seen that, due to the low NTG in someplace, the oil there is hard to be displaced. And because natural fractures are not developed, there is no obvious fingering of injected water. The oil production rate and watercut curves of these three producers are plotted in Fig. 41. It can be seen that, the watercut is generally lower than 50% and grows slowly because no water fingering happens, but the oil production rates decrease a little rapidly and are generally lower than 5 m3/d after five years, this is because the injected water does not make

4. Two application examples Based on the extended pEDFM framework, two examples are implemented to illustrate the applications of pEDFM. 4.1. Example 1 This example is mainly used to illustrate the significance of the flow barrier, that is, flow barriers can be embedded to construct complex reservoir boundary, the faults or interlayer in reservoirs. These cases are not effectively handled by classical EDFM, but not pEDFM. As shown in Fig. 29, there are two flow barriers and a highlyconductive fracture (20000mD) in the reservoir. An injector is at a constant injection rate of 1 m3/day, and a producer is at a constant bottom hole pressure of 10 MPa. Initial reservoir pressure is 20 MPa. Other physical parameters used in this example are the same as those in Table 1. The oil-phase pressure and saturation profiles at 100, 250 and 21

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 40. Pressure profiles of the lower layer at 1 year, 10 years, 30 years and 50.

22

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

Fig. 41. The oil production rate and watercut curves of these three producers.

a good performance to maintain the average reservoir pressure due to extremely low permeability of matrix and lack of natural fractures.

Declaration of competing interest We would like to submit the enclosed manuscript entitled " A modified projection-based embedded discrete fracture model (pEDFM) for practical and accurate numerical simulation of fractured reservoir ", which we wish to be considered for publication in “Journal of Petroleum Science & Engineering”. No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publi­ cation. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved the manuscript that is enclosed.

5. Conclusions This paper firstly points out three limitations among previous pEDFM, lying on lack of a practical method to select projected faces, the proficiency in the transmissibility formula of f-m connections, and poor performances in some cases. To resolve above limitations, correspond­ ing modifications to previous pEDFM are made. A practical algorithm namely “micro-translation method” to select proper projected faces is proposed, and several test cases are illustrated to show the validity of the algorithm. Then, the modified physical transmissibilities of additional fm connections and the extended pEDFM framework by considering additional f-f connections are presented. Some test cases are imple­ mented to illustrate the necessity and validity of these modifications and indicate that proposed extended pEDFM framework maybe the most complete up to now. Two examples are implemented to illustrate the 3D applications of the extended pEDFM framework.

Acknowledgements This work was supported by National Natural Science Foundation of China (No. U1762210, No. 51674273), the Fundamental Research Funds for the Central Universities (2019B07514), China Postdoctoral Science Foundation (No. BX20180381, No.2018M640218) and Science Foundation of China University of Petroleum, Beijing (No. 2462018YJRC015), and National Science and Technology Major Project of China (No. 2017ZX05013004-003 and No. 2017ZX05013002-003).

Author contribution Rao Xiang: Idea, Writing, Establishing models; Cheng Linsong: su­ pervision; Cao Renyi: supervision; Jia Pin: Reviewing and Editing; Liu Hao: English Checking and editing; Du Xulin: Commercial software running.

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.petrol.2019.106852. 23

X. Rao et al.

Journal of Petroleum Science and Engineering 187 (2020) 106852

References

Norbeck, J.H., McClure, M.W., Lo, J.W., Horne, R.N., 2016. An embedded fracturemodeling framework for simulation of hydraulic fracturing and shear stimulation. Comput. Geosci. 20 (1), 1–18. Panfili, P., Cominelli, A., 2014. Simulation of miscible gas injection in a fractured carbonate reservoir using an embedded discrete fracture model. In: Abu Dhabi International Petroleum Exhibition and Conference, 10–13 November. UAE, Abu Dhabi (July). Rao, X., Cheng, L., Cao, R., et al., 2019. An efficient three-dimensional embedded discrete fracture model for production simulation of multi-stage fractured horizontal well[J]. Eng. Anal. Bound. Elem. 106, 473–492. Rao, X., Cheng, L., Cao, R., et al., 2019. A modified embedded discrete fracture model to study the water blockage effect on water huff-n-puff process of tight oil reservoirs[J]. J. Pet. Sci. Eng. 106232. Rao, X., Cheng, L., Cao, R., Jia, P., Du, X., 2019. A fast embedded discrete fracture model based on proper orthogonal decomposition POD method. In: SPE Kuwait Oil & Gas Show and Conference. Rao, X., Cheng, L., Cao, R., Jia, P., Dong, P., Du, X., 2019. A modified embedded discrete fracture model to improve the simulation accuracy during early-time production of multi-stage fractured horizontal well. In: SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition. Ren, G., Jiang, J., Younis, R.M., 2018. A Model for coupled geomechanics and multiphase flow in fractured porous media using embedded meshes. Adv. Water Resour. 122, 113–130. Shakiba, M., Cavalcante Filho, J.S. de A., 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. Slough, K.J., Sudicky, E.A., Forsyth, P.A., 1999. Grid refinement for modeling multiphase flow in discretely fractured porous media. Adv. Water Resour. 23 (3), 261–269. Snow, D.T., 1969. Rock fracture spacings, openings, and porosities closure. Proc. Am. Soc. Civ. Eng. 95 (SM3), 873–883. Tene, M., Al Kobaisi, M.S., Hajibeygi, H., 2016. Algebraic multiscale method for flow in heterogeneous porous media with embedded discrete fractures (f-AMS). J. Comput. Phys. 321, 819–845. Tene, M., Bosma, S.B.M., Al Kobaisi, M.S., Hajibeygi, H., 2017. Projection-based embedded discrete fracture model (pEDFM). Adv. Water Resour. 105, 205–216. https://doi.org/10.1016/j.advwatres.2017.05.009. Wang, Y., Shahvali, M., 2016. Discrete fracture modeling using Centroidal Voronoi grid for simulation of shale gas plays with coupled nonlinear physics. Fuel 163, 65–73. Yu, 2017. Simulation of shale gas transport and production with complex fractures using embedded discrete fracture model. AIChE J.

Alboin, C., Jaffr�e, J., Roberts, J., Serres, C., 2002. Modeling fractures as interfaces for flow and transport in porous media. In: Chen, Ewing (Ed.), Fluid Flow and Transport in Porous Media, vol. 295. American Mathematical Society, pp. 13–24. Cao, R., Fang, S., Jia, P., Cheng, L., Rao, X., 2019. An efficient embedded discretefracture model for 2D anisotropic reservoir simulation. J. Pet. Sci. Eng. 174, 115–130. Dachanuwattana, S., Jin, J., Zuloaga-Molero, P., Li, X., Xu, Y., Sepehrnoori, K., Miao, J., 2018. Application of proxy-based MCMC and EDFM to history match a Vaca Muerta shale oil well. Fuel 220, 490–502. Hajibeygi, H., Karvounis, D., Jenny, P., 2011. A hierarchical fracture model for the iterative multiscale finite volume method. J. Comput. Phys. 230 (24), 8729–8743. Jiang, J., Younis, R.M., 2017. An efficient fully-implicit multislope MUSCL method for multiphase flow with gravity in discrete fractured media. Adv. Water Resour. 104, 210–222. Jiang, J., 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. https://doi.org/10.1016/j.advwatres.2017.09.017. Karimi-Fard, M., Durlofsky, L.J., Aziz, K., 2004. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Li, L., Lee, S.H., 2008. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reserv. Eval. Eng. 11 (4), 750–758. Lee, S.H., Lough, M.F., Jensen, C.L., 2001. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resour. Res. 37 (3), 443–455. Marcondes, F., Varavei, A., Sepehrnoori, K., 2010. An Element-Based Finite-Volume Method Approach for Naturally Fractured Compositional Reservoir Simulation. Brazilian Thermal Sciences Meeting. Uberlandia, Brazil, December. Martin, V., Jaffr�e, J., Roberts, J., 2005. Modeling fractures and barriers as interfaces for flow in porous media. SIAM 26, 1667–1691. Moinfar, A., Varavei, A., Sepehrnoori, K., et al., 2014. Development of an efficient embedded discrete fracturemodelfor 3d compositional reservoir simulation infractured reservoirs. SPE J. 19 (2), 289–303. SPE-154246-PA. Monteagudo, J.E.P., 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). Noorishad, J., Mehran, M., 1982. An upstream finite element method for solution of transient transport equation in fractured porous media [soil transport of pollutants]. Water Resour. Res. 18 (3), 588–596.

24