- Email: [email protected]
Multiscale mimetic method for two-phase flow in fractured media using embedded discrete fracture model
Advances in Water Resources 107 (2017) 180–190
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Multiscale mimetic method for two-phase flow in fractured media using embedded discrete fracture model Qingfu Zhang, Zhaoqin Huang, Jun Yao∗, Yueying Wang, Yang Li China University of Petroleum (East China), No. 66, Changjiang West Road, Huangdao District, Qingdao 266580, China
a r t i c l e
i n f o
Article history: Received 21 February 2017 Revised 22 June 2017 Accepted 24 June 2017 Available online 27 June 2017 Keywords: Multiscale mimetic method Embedded discrete fracture model Multiphase flow Fractured reservoirs
a b s t r a c t A multiscale mimetic method is developed for the simulation of multiphase flow in fractured porous media in the context of an embedded discrete fracture model (EDFM). The EDFM constructs independent grids for matrix and fracture system. Therefore, it is an efficient and practical flow model as it avoids the complicated unstructured grid subdivision and computing process. In order to extend the EDFM to field-scale applications, we integrate EDFM into a multiscale mimetic method. In this work, we use the multiscale basis functions to capture the detailed interactions between the fractures and the background. The multiscale basis functions are calculated numerically by solving EDFM on the local fine-grid with mimetic finite difference (MFD) method. The MFD method is conservative and robust, which makes it possible to deal with highly complex grid systems. Through combination of multiscale mimetic method and EDFM, this formulation can generate accurate velocity field and pressure field on the fine-scale grid more efficiently than the traditional methods. Numerical results are presented for verification of this multiscale mimetic approach for embedded discrete fracture media, and demonstrate its computational efficiency. The results show that this method is an accurate and efficient method for flow simulation in real-field fractured heterogeneous reservoirs. © 2017 Published by Elsevier Ltd.
1. Introduction Lots of engineering problems, such as reservoir development, groundwater transport, geothermal exploitation, largely depend on the information provided by simulation of flow through porous media. Apart from their intrinsic heterogeneous properties, these geological formations commonly contain complex fractures, with multiple scales and different conductivity properties. Given their essential influence on flow patterns, the fractures should be represented accurately. Among the methods targeted at describing fractures explicitly, single-porosity model (Ghorayeb and Firoozabadi, 20 0 0) regards fractures as a narrow high-permeability region. Fully discrete single-porosity models may contain up to 100 million grid cells. Hence, the simulation of such models is deemed intractable even with the advent of supercomputers. Discrete fracture model (DFM) (Huang et al., 2011a, 2011b; Karimi-Fard and Firoozabadi, 2001; Hauge and Aarnes, 2009; Hoteit and Firoozabadi, 2008) is based on conforming unstructured grid generation techniques which treat fractures as inner boundaries of the matrix cells. However, this technique would cause difficulties in grid generation and comput-
∗
Corresponding author. E-mail address: [email protected] (J. Yao).
http://dx.doi.org/10.1016/j.advwatres.2017.06.020 0309-1708/© 2017 Published by Elsevier Ltd.
ing progress. Especially when the distances between fractures are small, the grid generation often has a poor quality which will lead to miscalculation. As an alternative, the embedded discrete fracture model (EDFM) (Lee et al., 2001; Li and Lee, 2008; Yan et al., 2016) constructs independent grids for matrix and fracture system. That is, the discrete fracture network is embedded into the matrix structured grid system directly. Matrix and fractures are coupled by transfer functions. Therefore, EDFM is appealing for its ability to bypass the challenges related to unstructured grid (Moinfar et al., 2012; Zhou et al., 2014). Even though EDFM is an efficient flow model due to its ability to avoid the complicated unstructured grid subdivision and computing process, its ability for field-scale simulation is limited by numerical methods used to solve the flow model. Even after homogenization of small-scale fractures, the remaining degrees of freedom still exceed the controllable level of traditional numerical methods. This issue motivated the development of multiscale methods for EDFM. The multiscale methods discussed herein originate from the seminal paper (Hou and Wu, 1997; Efendiev and Wu, 2002), in which a governing equation with zero right hand side (i.e., homogeneous problem) is solved in each element to construct basis functions. It is similar to the ideas proposed earlier by Babuska et al. (1992) and Babuska and Osborn (1981). Then this
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
idea has been expanded by many researchers and spawned a series of relevant methods, including mixed multiscale finite-element (MsMFEM) (Chen and Hou, 20 03; Aarnes, 20 06, 20 05), Generalize multiscale finite element method (GMsFEM) (Efendiev et al., 2013, 2015), multiscale finite-volume (MsFV) method (Jenny et al., 20 03, 20 04; Hajibeygi and Jenny, 2009). All these methods construct basis functions by solving local flow problems to incorporate fine-scale effects into coarse-scale flow equations. Early versions of multiscale methods are proposed as a robust alternative to upscaling methods. Shortly thereafter, these methods have been expanded to model multiphase flow in porous media (Juanes, 2005; Krogstad et al., 2009; Efendiev et al., 2006). During the past few years, multiscale methods have been extended to simulate fluid flow in fractured media. They were initially proposed to model fractured porous media which treat fractures as narrow high-permeability regions (Natvig et al., 2011; Gulbransen et al., 2009). Later, instead of regarding fractures as volumetric element, discrete fracture model was integrated into multiscale finite element method (Zhang et al., 2016) and GMsFEM (Akkutlu et al., 2015). Recently, a combination of iterative multiscale finite volume method and hierarchical fracture model was presented (Hajibeygi et al., 2011; Matei et al., 2016). These multiscale methods have achieved satisfying results and made great contributions to the exploration of multiscale methods toward fractured porous media. In this work, we develop a multiscale mimetic method in the context of an embedded discrete fracture model (EDFM). The final matrix formulation of the multiscale mimetic method is similar to the multiscale mixed finite method. However, the discretization of the multiscale mimetic method is similar to finite volume method. Mimetic finite difference (MFD) method builds computational formulation on separate gridcells by introducing pressures located at the middle of the cell faces. This makes MFD could deal with any complex grid. To summarize, different from existing multiscale formulations, by combing multiscale mimetic method and EDFM, this method could avoid the complicated unstructured grid generation and computing process. Moreover, as discussed in many works (Skaflestad and Krogstad, 2008; Aarnes et al., 2008, 2014), multiscale mimetic method can deal with highly complex grid systems because of its applicability for the complex unstructured grid. Besides, the excellent local conservation property of MFD method makes multiscale mimetic method produce conservative velocity fields which are necessary for flow simulation. This paper proceeds as follows: we start by introducing the embedded discrete fracture model and its discretization. Next, we describe the multiscale mimetic method briefly and then extend it to deal with the EDFM system. In the numerical experiment section, several 2D and 3D numerical test cases are presented to validate the correctness and effectiveness of the multiscale method. Finally the concluding remarks were given in the final section.
181
Fracture system:
−∇ · (Kf λf · ∇ pf ) = qf −
qmf q − ff δff Vf Vf
(2)
Here, Ki (i = m, f) stands for the permeability tensor; pi, Vi represent pressure, volume of grid cells, respectively; qi = qin + qiw is sink/source term (the subscript w denotes the wetting phase and n denotes the non-wetting phase); λi = λin + λiw is the total mobility where λiα = kriα /μα is the mobility of phase α , which depends on relative permeabilitykriα and viscosity μα ; qmf is transfer flow between matrix and fracture; qff is transfer flow between intersecting fractures; δ mf = 1 if matrix cell contains fracture cells, else δ mf = 0; δ ff = 1 if fracture cell intersects with another fracture cells, else δ ff = 0. The pressure in matrix is assumed to be continuous, therefore, the crossflow term can be written in this form
qmf = −Tmf ( pm − pf )
(3)
Here, Tmf = kmf λmf Amf /d, where Amf is interfacial area, d is equivalent distance between matrix cell and fracture cell, λmf is total mobility decided by upstream weight method, and kmf is harmonically averaged permeability. The crossflow term qff of intersected fractures (Fig. 1) is given by (Karimi-Fard et al., 2004)
qff = Tff pfi − pf j
(4)
, T = k d /d ; Here, Tmf = λff Tfi Tfj /(Tfi + Tfj ), where Tfi = kfi dfi /d i j fj fj fj λff is total mobility decided by upstream weight method; df and kf stand for fracture aperture and permeability respectively; d represents the average normal distance between the centre of fracture segments and the intersections of fractures. The saturation equations are written as: Matrix system
φm
∂ Smw q + ∇ · vmw = qmw + mfw δm f ∂t Vm
vmw = fmw [vm + Km λmn · ∇ pmc ]
(5) (6)
Fracture system
φf
∂Sfw q f fw q + ∇ · v f w = q f w − mfw δm f − δ ∂t Vf Vf f f
v f w = f f w [v f + K f λ f n · ∇ p f c ]
(7) (8)
Here, φ i denotes the porosity, fiw = λiw /λi denotes the fractional flow function, and pic is the capillary pressure. In this work, considering the computational efficiency, we use IMPES strategy to solve the coupled system. IMPES means the flow equations are solved implicitly to obtain the velocity field. Then the velocity field is employed to solve the saturation equations. 2.2. Discretization
2. Embedded discrete fracture modelling 2.1. Mathematical model The flow system for incompressible and isothermal two-phase flow without considering the influence of gravity in fractured media is governed by pressure and saturation equations. The pressure equations are written as Matrix system:
−∇ · (Km λm · ∇ pm ) = qm +
qmf δmf Vm
Matrix grid and fracture grid are constructed to solve the flow equations (Fig. 2). The fracture grid is embedded into matrix grid independently. Therefore, this method can alleviate the gridding complexities. To efficiently handle polyhedral cells, MFD is used as our discretization method. As shown in Fig. 3, Ak is the interface between the matrix cells, nk is the normal to the interface, xik is the vector pointing from cell centroid to the interface centroid. f Let vmi = [vmi1 , vmi2 , · · ·,vmin ]T be the fluxes of interfaces of i , f
(1)
pemi and pmk is the average pressures located at the gridcell centres and interfaces respectively. These quantities are related through a transmissibility matrixTmi , that is
182
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
li1 lj1
lj2
li2
Fig. 1. Schematic of the intersecting fractures in embedded discrete fracture model.
(a) matrix grid
(b) fracture grid
Fig. 2. Schematic of the independent grids in embedded discrete fracture model.
The matrix system can be coupled with the fracture system to form the resulting equation
⎡ nk
Ak Ði
x ik
pei
pfk
⎢CmT ⎢ ⎢ T ⎢Dm ⎢ ⎣0
pej
Ðj
vmf i = Tmi · ei pemi − pmf k
(9)
1]T .
Here, ei = [1, 1, ..., By using the divergence theorem, we can derive the following equation from Eq. (1)
k=1
i
qmi d + qm f δm f
(10)
In this work, the dimension of fractures is reduced to (d-1)dimension. By using implicit difference method, the flow equations in the fracture system are discretized as
Tξ i+ 1 pfi+1 − pfi − Tξ i− 1 pfi − pfi−1 = ffi + qmfi + qffi δffi 2
2
Here, Tξ i± 1 = Kf λf 0.5( ξ 2
0
0
Tmf1 + Tmf2
0
−Tmf1
−Tmf2
0
0
0
0
Tmf1
0
Tf1 − Tmf1 − Tff
Tff
Tmf2
0
Tff
Tf2 − Tmf2 − Tff
⎡ ⎤ ⎡ ⎤ vm 0 ⎢ pm ⎥ ⎢ f m ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ×⎢ ⎢πm ⎥ = ⎢0 ⎥ ⎣ pf1 ⎦ ⎣ ff1 ⎦
Fig. 3. Schematic of mimetic finite difference method.
vmf k =
Dm
0
face centroid cell centroid
n
−Cm
Bm
dfi
i±1 + ξ i )
(11)
, and ffi = Vfi qfi , where Vfi and dfi
represent the volume and aperture of i thfracture element, respectively.
pf2
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(12)
ff2
Where Bm and Cm are block diagonal for which the gridblocks are Tm−1 and an Ne × 1 vector with all entries equal one, where Tmi i denotes transmissibility between cell central and cell boundaries of gridcell i . This denotes that the mimetic method works directly on each gridcell. Ne denotes the number of grids. Tmfi is the matrix of transfer coefficient between fracture cell and matrix cell, Tff is the matrix of transfer coefficient between two intersected fractures. The sought unknowns are the face-velocity vm , the matrix gridcell pressures pm , matrix face-pressure π m , the fracture pressure pfi . Developing an efficient approach for Eq. (12) is a challenging problem in field applications. One reason is the system may contain up to several millions of sought unknowns in field-scale applications. Another reason is the great differences between fracture and matrix properties would deteriorate the value of condition numbers. Apparently traditional upscaling method cannot be used to solve this system because of the high resolution fractures, which have essential influence on whole flow patterns. This creates a mo-
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) matrix grid
183
(b) fracture grid
Fig. 4. Schematic of multiscale grids in multiscale mimetic method.
⎧ ⎨ ωi (x ), x ∈ i ∇ · ψi j = −ω j (x ), x ∈ j ⎩ 0, x∈ / i j
In ij with ψ ij • n = 0 ∀x ∈ ∂ ij . No-flow boundaries are imposed. We should notice that it is generally not convergent as the mesh size tends to zero when K is an anisotropic tensor. A homogenization-based mixed multiscale finite element could be applied to avoid this problem (Arbogast, 2011). The ωi (x) must satisfy ωi (x )dx = 1 so that there is a unit i average flow over the ij between two coarse gridcells. The most simplest option is ωi (x) = 1/|i |, in this work, in order to avoid the computational error resulting from homogenization, the ωi (x) is chosen as
Fig. 5. Schematic description of basis functions.
tivation for the multiscale method, which have an attractive advantage of solving coarse problems while honouring fine-scale information. The introduction of this method is presented in the following section.
ω (x )i =
In this part, we briefly describe the multiscale mimetic framework, an efficient multiscale procedure for fractured porous media. This method is formulated using a coarse-scale grid and a finescale grid. The multiscale mimetic method first superposes a finescale grid on both matrix and fractures as shown in Fig. 3. The petrophysical properties of matrix and fracture are given in the fine-scale grid. The coarse-scale grid is introduced on top of the fine-scale grid as seen in Fig. 4. Each block is composed of a connected collection of fine-scale gridcells. In this work, the fracture gridcells are also decomposed to coarse-scale grid and fine-scale grid. Over this coarse-scale grid, two approximation spaces are defined: each coarse gridcell corresponds to a pressure approximation space while each interface corresponds to a flux approximation space. Next, all pairs of adjacent blocks are detected for the construction of multiscale basis functions. This is a crucial step of multiscale method since the fine-scale features and the interaction between the fractures and the matrix are accounted for with multiscale functions. Consider two coarse-scale grid cells i and j (Fig. 5) and let ij consists of i and j , that is ij = i ∪ ij ∪j .
ij is the interface between i and j . The basis functions ϕ ij for pressure and ψ ij for velocity are constructed by solving following local flow equations (Aarnes et al., 2005)
(13)
⎧ ⎪ σ x / σ ξ d ξ , qdx = 0 ( ) ( ) ⎨
⎪ ⎩q(x )/
i
i
q(ξ )dξ
i
i
(15) qdx = 0
Where σ (x) = trace(K(x))/d. By using MFD, the hybrid system of Eqs. (13)–(15) can be written as (Aarnes et al., 2006)
⎡
3. Multiscale mimetic method for fractured porous media
ψi j = −λt K ∇ ϕi j
(14)
−C
B
⎣C T
0
DT
0
⎤ ⎤ ⎡ ψi j 0 0 ⎦⎣ϕi j ⎦ = ⎣ωi θi j ⎦ π 0 0
D
⎤⎡
(16)
Here B and C are block diagonal for which the gridblocks are Tm−1 and an ne × 1 vector with all entries equal one, Tm is transmissibility matrix, ne is the number of gridcells in the block ij , θ ij = 1 if x ∈ i , θ ij = −1 if x ∈ j . The velocity basis function ψ ij is represented as a vector of fluxes on the set of faces inside the coarse gridcells. Similarly, the pressure basis function ϕ ij is represented as a vector of fine-scale gridcells pressure. When the coarse blocks contain fractured regions, the fractures are represented by EDFM. This model is effective for constructions of basis functions, especially when the fractures are close to each other, which may lead to a poor unstructured grid generation and calculation error. The flowing system is solved for basis functions
−∇ · Km λm · ∇ϕimj − ψ f m δmf = ω (x )i
−∇ · Kf λf · ∇ϕifj + ψ m f δm f + ψ f f δ f f = ω (x )i
(17) (18)
The basis functions, including the fracture basis functions ψ fm and ψ ff , are calculated using a linear system like Eq. (12). The fracture velocity basis function ψ fm is represented as a vector of fluxes on the faces between fracture and matrix, the fracture pressure baf sis function ϕi j is represented as a vector of fine-scale fracture cells
184
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Fracture
(a) Coarse grid cell with partially overlapping fracture
(b)
(c)
(d)
(f)
Fig. 6. Illustration of all four pressure basis functions of a coarse cell. (b) is basis function of bottom edge, (c) is the basis function of right edge, (d) is the basis function of upper edge, (f) is the basis function of left edge.
important factors (such as matrix, fractures) in the construction of the global equation system (Fig. 6). If we construct multiscale basis functions in a larger domain, it means the multiscale basis functions are constructed using oversampling technique. Multiscale basis functions normally only have support in the corresponding coarse gridcells. Oversampling means that we extend the support of the multiscale basis functions (see Fig. 7), that is, we include more fine gridcells when computing the basis function. The use of oversampling is particularly attractive in cases of highly heterogeneous reservoirs. In fractured blocks, the use of oversampling technique is more important than its application in the matrix blocks. By using oversampling technique, we can extend the support of the basis functions and lower the impact of boundary. Specifically, the multiscale basis functions ψ Ej (also called oversampling functions) are constructed in the larger domain KE ⊃K (see Fig. 7), that is, the support of the multiscale basis functions is expanded from K to KE . We then form the actual ψ i by linear combination of ψ Ej ,
Fig. 7. Illustration of oversampled region.
pressures, similar to the matrix basis functions. Here, basis function ψ ij is constructed for matrix-matrix impacts, ψ fm is for the matrix-fracture coupling, ψ ff is for fracture-fracture interactions. Multiscale basis functions are applied to capture the impacts of all
ψi =
e
ci j ψ jE
(19)
j=1
where cij are determined by the corresponding boundary conditions of target block.
2-D problem
3-D problem
e5 I
e5 e4
e1
I e3
e2
e4
e1
e3
Fig. 8. Schematic with intersecting fractures at I.
e2
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
185
No-flow
1. 0
( 0. 2, 0. 52)
( 0. 8, 0. 52)
P=0Mpa
P=1Mpa
( 0. 52, 0. 8)
( 0. 52, 0. 2) 0
1. 0
No-flow
(a) Fractures model
(c) EDFM pressure
(b) Multiscale grid
(d) Multiscale coarse-scale results
(e) Multiscale fine-scale results
Fig. 9. Reference pressure solution and multiscale pressure solutions.
Fig. 10. Velocity and pressure distribution along y = 0.52 L of EDFM and multiscale solutions.
Once the multiscale basis functions are constructed, we then assemble the coarse-scale system. To form the coarse system, we split the basis functions into two parts
ψi j = ψiHj − ψ jiH
(20)
Where
organized as ψ ff and ψ fm . Then the global coarse system reads
⎡
c −Cm
Dcm
0
0
Tmc f 1 + Tmc f 2
0
−Tmc f 1
−Tmc f 2
0
0
Bcm
⎢C c T ⎢ m ⎢ cT ⎢D m ⎢ ⎢ 0 ⎣
Tmc f 1
0
−ψi j (E ), E ∈ j ψ i j ( E ) , E ∈ i j \ j H H ψi j ( E ) = · ψ ji (E ) = 0, E∈ / i 0, E∈ / j (21) Next, we assemble all the velocity basis functions ψiHj as
columns in ψ and assembled the pressure basis function as columns in matrix . Likewise, the fractures basis functions are
⎡
vc
⎤
0
Tmc f 2
⎡
0
⎤
⎢ pc ⎥ ⎢ f c ⎥ ⎢ c⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ×⎢ ⎢π ⎥ = ⎢0 ⎥ ⎣ pc ⎦ ⎣ f c ⎦ f1
f1
pcf2
c ff2
0
0 T fc1
−
Tmc f 1 T fcf
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
0 −
T fcf
T fcf Tf2c − Tmc f 2 − T fcf
(22)
Here Bc = ψ T Bf ψ , Cc = ψ T Cf I, fc = IT qf , and Dc = ψ T Cf J .I prolongates blocks to cells, that is, Iij =1 if blocks number j contains cell number i and Iij =0 otherwise. Similarly, J prolongates block
186
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Permeability field
(b) Multiscale grid
Fig. 12. Multiscale flux error for different coarse meshes.
(c) EDFM pressure
(d) Multiscale solution
Fig. 11. Pressure distribution of EDFM and multiscale solutions.
faces to fine faces, Jij =1 if coarse face j contains the fine face i and Jij =0 otherwise. Once the Eq. (22) is solved, the fine-scale velocity and pressure can be constructed as (Krogstad et al., 2009, 2012)
vf ≈ ψ vc
(23)
pf ≈ I pc + Dλ vc
(24)
The fracture velocity can be obtained by
v ≈ ψ ff vc ff
(25)
¯ 0 /λ ¯ i for block i, λ ¯0 Here Dλ is a diagonal matrix with entries λ i i stands for the mean mobility in the construction of the basis func¯ i stands for the mean mobility in the current time. tions and λ 4. Approximation of the saturation Fig. 13. Total CPU time used by multiscale method and reference method.
The saturation functions are discretized with a finite volume method.
No-flow
( 0. 2, 0. 25)
0
1. 0
( 0. 8, 0. 25)
No-flow
(a) Fractures model
(c) Reference saturation field
P=0Mpa
P=1Mpa
0. 5
1. 0
(b) Multiscale grid
(d) Multiscale saturation field
Fig. 14. Saturation distribution of EDFM and multiscale solutions.
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Here,
Ffk (Sfw ) =
Afk
187
ffw (Sfw )k (vf · nk + Km λfn · ∇ pfc · nk )dA
(28)
In above Eqs (26)–(28), n stands for time step, Ni denotes number of faces, fiw (S) is fractional flow, which is obtained by upstream scheme. Upstream-weighting method is applied to deal with multiple intersecting fractures. We assume that there are NI fracture elements ei intersect at Iin Fig. 8. Here, fw,ei and vf,ei stand for fractional function and flux respectively. n is defined as 0 < n < NI , such that
(a) matching grid
(b) Reference saturation field
vf,ei ≤ 0, 0 < i ≤ N (efflux ) vf,ei ≥ 0, N < i < NI (influx )
(29)
Considering the law of conservation of volumetric and mass, one gets NI
vf,ei = −
i=N+1 NI
N
vf,ei
(30)
i=1
fw,i vf,ei = −
i=N+1
N
fw,I vf,ei = − fw,I
i=1
N
vf,ei
(31)
i=1
Therefore, the fractional function at the intersections is
NI
fw,I = −
v
f i=N+1 w,i f,ei N i=1 f,ei
v
(32)
5. Numerical experiments
(c) Multiscale grid
(d) Multiscale saturation field Several numerical experiments conducted in this work are given in this part. Firstly, the multiscale mimetic method is validated for fractured reservoirs through a comparison with a fine-scale simulation of fractured media. Then a simulation of fractured media with heterogeneous matrix is presented to investigate efficiency and robustness of the multiscale method. The third example is set to validate the multiscale method for multiphase flow in fractured reservoirs. Finally, we extended our multiscale method to 3D fractured model and demonstrated its powerful industrial potential. 5.1. Verification of the multiscale method for fractured reservoirs
(e) Water-cut curve of irregular fractured media Fig. 15. Saturation distribution of DFM and multiscale solutions.
For the matrix system Nm φmi n+1 1 n n (Smwi − Smw ) + ) θ F (Sn+1 ) + (1 − θ )Fmk (Smw i t |mi | k=1 mk mw n = qmwi (Smw i) +
n qm f w ( Sw ) δmf Vmi
(26)
In the first scenario, a 2D homogeneous reservoir was considered, where two orthogonal fractures are located at the centre of the L × L region. The permeability of both fractures satisfy Kf /Km = 105 . We impose no-flow boundary conditions on the top and bottom while the left side and right side are set to constant pressure boundaries p = 1 and p = 0 respectively. The fine-scale reference pressure is obtained with 40 × 40 grid cells. For the simulation of multiscale mimetic-EDFM, 10 × 10 coarse grid cells are employed (Fig. 9(b)). Numerical results of EDFM and multiscale method are presented in Figs. 9 and 10. From the comparison of the results, we can observe that the pressure distribution obtained by reference solution can be reproduced by the fine-scale results of the multiscale method. 5.2. Efficiency analysis of the multiscale method
For the fracture system Nf φfi n+1 n 1 n (Sfwi − Sfwi ) + ) θ F (Sn+1 ) + (1 − θ )Ffk (Sfw t |fi | k=1 fk fw n n qf f w Sfw qm f w ( Sw ) n = qfw j (Sfw ) − δ − δff mf j
Vf j
Vf j
(27)
The second test was presented to explore the efficiency and robustness of the multiscale method for more general fractured reservoirs. In this test, a 2D fractured reservoir model with a heterogeneous permeability field is considered (Fig. 11). Boundary conditions of this model are same as the previous test. The finescale pressure is obtained with 40 × 40 grid cells and for the multiscale method 10 × 10 coarse grid is employed.
188
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Fine-scale matrix grid
(b) Fine-scale fracture grid
(c) Coarse-scale matrix grid
(d) Coarse-scale fracture grid
Fig. 16. Schematic of multiscale grid and fractured model.
In this section, we use relative L2 norm to measure the velocity discrepancies:
ev =
vf − vms 22 vf 22
(33)
where vf denotes reference flux solutions while vms stands for the flux calculated with multiscale mimetic method. Fig. 12 shows the flux L2 error for different coarse grid. It demonstrates that the coarse grid has a significant impact on the computation accuracy. Our multiscale mimetic method generally gives better results if the coarse grid has approximately the unit aspect ratio. Besides, it seems to perform best on the coarse grid with medium number coarse gridcells. Fig. 13 compares the computing time of traditional numerical (MFD) method and multiscale method. It shows that multiscale method cost different time with different coarse grid. Overall, multiscale method is an efficient numerical method. 5.3. Multiphase flow in fracture network This section is set to validate our multiscale method for multiphase flow in fractured reservoirs. First, a heterogeneous 2L × L media embedded with one fracture is considered (Fig. 14). The permeability of fracture satisfies K f /K¯ m = 105 , where K¯ m is the matrix average permeability. No-flow conditions are employed on the upper and lower boundaries. Constant pressure differential was imposed along the horizontal direction. The fine-scale grid contains 40 × 20 cells and for multiscale method 10 × 5 coarse grid cells were employed. Then one-well injection and one-well production model with irregular shape (Fig. 15) is presented. No-flow conditions are used on all boundaries in this test. The matching grid used by discrete fracture model contains 577gridcells. The mul-
tiscale grid contains 40 coarse-scale gridcells and 525 fine-scale gridcells respectively. An injection well located at the bottom left corner while the production well located at the top right corner. Linear relative permeability is employed for flow in porous media. In this study, water cut is defined as fw =Qw /(Qo +Qw ), where Qw is the water flow rate and Qo is the oil flow rate at the outlet. Figs. 14 and 15 shows saturation distributions obtained with MFD method and multiscale mimetic method. The saturation error of Fig. 14 in L1 norm is 0.0981. Both test cases indicate that the finescale and multiscale results are nearly identical, even though the water-cut curve in reference method presents bigger value. This confirms that the multiscale mimetic-EDFM can track the reference solution very closely in the case of multiphase flow. 5.4. Multiphase flow in 3D fractured reservoirs Finally, we extend the multiscale mimetic method to 3D fractured reservoir model and take a crucial step towards the final real-field simulation. In this case, a 3D homogeneous media containing 4 fractures are considered. The fine-scale grid is 10 × 10 × 4 and the multiscale method was run using coarsening 5 × 5 × 2 grid (Fig. 16). The injection well is located at the bottom left corner (i.e.x = 0 m, y = 0 m) and the production is located in the top right corner (i.e.x = 100 m, y = 100 m). The simulation results depicted in Fig. 17 including the saturation profiles after 0.5 PVI, the pressure distributions and the watercut curve of the fractured reservoir model. These simulation results demonstrate the powerful industrial potential of our multiscale method. Based on the effective multiscale method, we can conduct various analyses such as production prediction, risk assessment and develop effective management strategy.
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Reference saturation maps after 0.5 PVI
(c) Reference pressure distribution
189
(b) Multiscale saturation maps after 0.5 PVI
(d) Multiscale pressure distribution
(e) Water-cut curve of 3D fractured media Fig. 17. Schematic of multiscale grid and fractured model.
6. Concluding remarks In this work we develop a multiscale mimetic-EDFM for multiphase flow in fractured porous media. By using embedded discrete fracture model to construct the multiscale basis functions, we can capture the detailed interactions between the fractures and the background. Moreover, the use of embedded discrete fracture model enables us to avoid the challenges related to unstructured grid which results in a complicated grid generation and computing progress. Our multiscale mimetic-embedded discrete fracture model is under the framework of mimetic finite difference method, the excellent conservation property and robustness of MFD method makes the multiscale method can deal with complex and highly unstructured grids systems (Aarnes et al., 2008, 2014). Some nu-
merical results presented in the paper demonstrate that the multiscale mimetic-embedded discrete fracture is an accurate and efficient method for fractured reservoir. The construction of multiscale basis functions can be completed with parallel computing method to further reduce memory requirements. Several multiscale methods could get an increased accuracy and stability by iterative methods (Krogstad et al., 2009; Hajibeygi et al., 2011, 2008; Kuenze and Lunati, 2011) and this is a good inspiration for us. The design of an iterative method that improves the solution accuracy systematically is our ongoing research. The above lead to the conclusion that the multiscale mimeticembedded discrete fracture method is an important multiscale development for the simulation of fractured reservoirs.
190
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Acknowledgements The authors gratefully acknowledge support from National Natural Science Foundation of China (51404292), the Fundamental Research Funds for the Central Universities (17CX06007), National Science and Technology Major Project (2016ZX05060-010). References Aarnes, J.E., 2006. On the use of a mixed multiscale finite element method for greater flexibility and increased speed or improved accuracy in reservoir simulation. SIAM J. Multiscale Model. Simul. 2 (3), 421–439. Aarnes, J.E., Kippe, V., Lie, K.A., 2005. Mixed multiscale finite elements and streamline methods for reservoir simulation of large geomodels. Adv. Water Res. 28 (3), 257–271. Aarnes, J.E., Krogstad, S., Lie, K.A., 2006. A hierarchical multiscale method for two-phase flow based upon mixed finite elements and nonuniform coarse grids. SIAM J. Multiscale Model. Simul. 5 (2), 337–363. Aarnes, J.E., Krogstad, S., Lie, K.A., 2008. Multiscale mixed/mimetic methods on corner-point grids. Comput. Geosci. 12 (3), 297–315. Aarnes J.E., Lie K.A., Laptev V., Krogstad S. (2014). Multiscale mixed methods on corner-point grids: mimetic versus mixed subgrid solvers. Akkutlu, I., Efendiev, Y., Vasilyeva, M., 2015. Multiscale model reduction for shale gas transport in fractured media. Comput. Geosci. 1–21. Arbogast, T., 2011. Homogenization-based mixed multiscale finite elements for problems with anisotropy. SIAM J. Multiscale Model. Simul. 9 (2), 624–653. Babuska, I., Caloz, G., Osborn, J.E., 1992. Special finite element methods for a class of second order elliptic problems with rough coefficients. SIAM J. Numer. Anal. 31 (4), 945–981. Babuska, I., Osborn, J.E., 1981. Generalized finite element methods: their performance and their relation to mixed methods. SIAM J. Numer. Anal. 20 (3), 510–536. Chen, Z., Hou, T.Y., 2003. A mixed multiscale finite element method for elliptic problems with oscillating coefficients. Math. Comput. 72 (242), 541–576. Efendiev, Y., Galvis, J., Hou, T.Y., 2013. Generalized multiscale finite element methods (GMsFEM). J. Comput. Phys. 251 (23), 116–135. Efendiev, Y., Ginting, V., Hou, T., Ewing, R., 2006. Accurate multiscale finite element methods for two-phase flow simulations. J. Comput. Phys. 220 (1), 155–174. Efendiev, Y., Lee, S., Li, G., Yao, J., Zhang, N., 2015. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method. GEM - Int. J. Geomath. 6 (2), 1–22. Efendiev, Y.R., Wu, X.H., 2002. Multiscale finite element for problems with highly oscillatory coefficients. Numer. Math. 90 (3), 459–486. Ghorayeb, K., Firoozabadi, A., 20 0 0. Numerical study of natural convection and diffusion in fractured porous media. Spe J. 5 (5), 12–20. Gulbransen, A, Hauge, V.L., Lie, K.A., 2009. A multiscale mixed finite-element method for vuggy and naturally fractured reservoirs. SPE J. 15 (15), 395–403. Hajibeygi, H., Bonfigli, G., Hesse, M.A., Jenny, P., 2008. Iterative multiscale finite-volume method. J. Comput. Phys. 227 (19), 8604–8621. Hajibeygi, H., Jenny, P., 2009. Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media. J. Comput. Phys. 228 (14), 5129–5147. 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.
Hauge, V.L., Aarnes, J.E., 2009. Modeling of two-phase flow in fractured porous media on unstructured non-uniformly coarsened grids. Transp. Porous Media 77 (3), 373–398. Hoteit, H., Firoozabadi, A., 2008. An efficient numerical model for incompressible two-phase flow in fractured media. Adv. Water Res. 31 (6), 891–905. Hou, T.Y., Wu, X.H., 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134 (1), 169–189. Huang, Z.Q., Yao, J., Wang, Y.Y., 2011b. Numerical study on waterflooding development of fractured reservoir with discrete-fracture model. Chinese J Comput Phys 28, 148–156. Huang, Z.Q., Yao, J., Wang, Y.Y., Tao, K., 2011a. Numerical study on two-phase flow through fractured porous media. Sci. China Technol. Sci. 54 (9), 2412–2420. Jenny, P., Lee, S.H., Tchelepi, H.A., 2003. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. J. Comput. Phys. 187 (1), 47–67. Jenny, P., Lee, S.H., Tchelepi, H.A., 2004. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. SIAM J. Multiscale Model. Simul. 3 (1), 50–64. Juanes, R., 2005. A variational multiscale finite element method for multiphase flow in porous media. Finite Elements Anal. Des. 41 (7–8), 763–777. Karimi-Fard, M., Durlofsky, L.J., Aziz, K., 2004. Anefficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Karimi-Fard, M., Firoozabadi, A., 2001. Numerical simulation of water injection in 2D fractured media using discrete fracture model. SPE Reservoir Eval. Eng.. Krogstad, S., Lie, K.A., Nilsen, H., Natvig, J., Skaflestad, B., Aarnes, J., 2009. A multiscale mixed finite element solver for three phase black oil flow. SPE Reservoir Simulation Symposium. Krogstad, S., Lie, K.A., Skaflestad, B., 2012. Mixed multiscale methods for compressible flow. Ecmor Xiii, European Conference on the Mathematics of Oil Recovery. Kuenze, R., Lunati, I., 2011. An adaptive iterative multiscale finite volume method to simulate density- driven flow instabilities. In: AGU Fall Meeting, p. 1476. 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. 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. Matei, Tene, 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. Moinfar, A., Varavei, A., Sepehrnoori, K., Johns, R.T., 2012. Development of a novel and computationally-efficient discrete-fracture model to study ior processes in naturally fractured reservoirs. SPE Improved Oil Recovery Symposium. Natvig, J.R., Skaflestad, B., Bratvedt, F., Bratvedt, K., Lie, K.A., Laptev, V., Khataniar, S.K., 2011. Multiscale mimetic solvers for efficient streamline simulation of fractured reservoirs. SPE J. 16 (4), 880–888. Skaflestad, B., Krogstad, S., 2008. Multiscale/mimetic pressure solvers with near-well grid adaptation. Ecmor. Yan, X., Huang, Z., Yao, J., Li, Y., Fan, D., 2016. An efficient embedded discrete fracture model based on mimetic finite difference method. J. Petrol. Sci. Eng. 145, 11–21. Zhang, Q., Yao, J., Huang, Z., 2016. An efficient multiscale mixed finite element method for modelling flow in discrete fractured reservoirs. ECMOR XV, European Conference on the Mathematics of Oil Recovery. Zhou, F., Shi, A., Wang, X., 2014. An efficient finite difference model for multiphase flow in fractured reservoirs. Petrol. Explorat. Develop. 41 (2), 262–266.
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Multiscale mimetic method for two-phase flow in fractured media using embedded discrete fracture model Qingfu Zhang, Zhaoqin Huang, Jun Yao∗, Yueying Wang, Yang Li China University of Petroleum (East China), No. 66, Changjiang West Road, Huangdao District, Qingdao 266580, China
a r t i c l e
i n f o
Article history: Received 21 February 2017 Revised 22 June 2017 Accepted 24 June 2017 Available online 27 June 2017 Keywords: Multiscale mimetic method Embedded discrete fracture model Multiphase flow Fractured reservoirs
a b s t r a c t A multiscale mimetic method is developed for the simulation of multiphase flow in fractured porous media in the context of an embedded discrete fracture model (EDFM). The EDFM constructs independent grids for matrix and fracture system. Therefore, it is an efficient and practical flow model as it avoids the complicated unstructured grid subdivision and computing process. In order to extend the EDFM to field-scale applications, we integrate EDFM into a multiscale mimetic method. In this work, we use the multiscale basis functions to capture the detailed interactions between the fractures and the background. The multiscale basis functions are calculated numerically by solving EDFM on the local fine-grid with mimetic finite difference (MFD) method. The MFD method is conservative and robust, which makes it possible to deal with highly complex grid systems. Through combination of multiscale mimetic method and EDFM, this formulation can generate accurate velocity field and pressure field on the fine-scale grid more efficiently than the traditional methods. Numerical results are presented for verification of this multiscale mimetic approach for embedded discrete fracture media, and demonstrate its computational efficiency. The results show that this method is an accurate and efficient method for flow simulation in real-field fractured heterogeneous reservoirs. © 2017 Published by Elsevier Ltd.
1. Introduction Lots of engineering problems, such as reservoir development, groundwater transport, geothermal exploitation, largely depend on the information provided by simulation of flow through porous media. Apart from their intrinsic heterogeneous properties, these geological formations commonly contain complex fractures, with multiple scales and different conductivity properties. Given their essential influence on flow patterns, the fractures should be represented accurately. Among the methods targeted at describing fractures explicitly, single-porosity model (Ghorayeb and Firoozabadi, 20 0 0) regards fractures as a narrow high-permeability region. Fully discrete single-porosity models may contain up to 100 million grid cells. Hence, the simulation of such models is deemed intractable even with the advent of supercomputers. Discrete fracture model (DFM) (Huang et al., 2011a, 2011b; Karimi-Fard and Firoozabadi, 2001; Hauge and Aarnes, 2009; Hoteit and Firoozabadi, 2008) is based on conforming unstructured grid generation techniques which treat fractures as inner boundaries of the matrix cells. However, this technique would cause difficulties in grid generation and comput-
∗
Corresponding author. E-mail address: [email protected] (J. Yao).
http://dx.doi.org/10.1016/j.advwatres.2017.06.020 0309-1708/© 2017 Published by Elsevier Ltd.
ing progress. Especially when the distances between fractures are small, the grid generation often has a poor quality which will lead to miscalculation. As an alternative, the embedded discrete fracture model (EDFM) (Lee et al., 2001; Li and Lee, 2008; Yan et al., 2016) constructs independent grids for matrix and fracture system. That is, the discrete fracture network is embedded into the matrix structured grid system directly. Matrix and fractures are coupled by transfer functions. Therefore, EDFM is appealing for its ability to bypass the challenges related to unstructured grid (Moinfar et al., 2012; Zhou et al., 2014). Even though EDFM is an efficient flow model due to its ability to avoid the complicated unstructured grid subdivision and computing process, its ability for field-scale simulation is limited by numerical methods used to solve the flow model. Even after homogenization of small-scale fractures, the remaining degrees of freedom still exceed the controllable level of traditional numerical methods. This issue motivated the development of multiscale methods for EDFM. The multiscale methods discussed herein originate from the seminal paper (Hou and Wu, 1997; Efendiev and Wu, 2002), in which a governing equation with zero right hand side (i.e., homogeneous problem) is solved in each element to construct basis functions. It is similar to the ideas proposed earlier by Babuska et al. (1992) and Babuska and Osborn (1981). Then this
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
idea has been expanded by many researchers and spawned a series of relevant methods, including mixed multiscale finite-element (MsMFEM) (Chen and Hou, 20 03; Aarnes, 20 06, 20 05), Generalize multiscale finite element method (GMsFEM) (Efendiev et al., 2013, 2015), multiscale finite-volume (MsFV) method (Jenny et al., 20 03, 20 04; Hajibeygi and Jenny, 2009). All these methods construct basis functions by solving local flow problems to incorporate fine-scale effects into coarse-scale flow equations. Early versions of multiscale methods are proposed as a robust alternative to upscaling methods. Shortly thereafter, these methods have been expanded to model multiphase flow in porous media (Juanes, 2005; Krogstad et al., 2009; Efendiev et al., 2006). During the past few years, multiscale methods have been extended to simulate fluid flow in fractured media. They were initially proposed to model fractured porous media which treat fractures as narrow high-permeability regions (Natvig et al., 2011; Gulbransen et al., 2009). Later, instead of regarding fractures as volumetric element, discrete fracture model was integrated into multiscale finite element method (Zhang et al., 2016) and GMsFEM (Akkutlu et al., 2015). Recently, a combination of iterative multiscale finite volume method and hierarchical fracture model was presented (Hajibeygi et al., 2011; Matei et al., 2016). These multiscale methods have achieved satisfying results and made great contributions to the exploration of multiscale methods toward fractured porous media. In this work, we develop a multiscale mimetic method in the context of an embedded discrete fracture model (EDFM). The final matrix formulation of the multiscale mimetic method is similar to the multiscale mixed finite method. However, the discretization of the multiscale mimetic method is similar to finite volume method. Mimetic finite difference (MFD) method builds computational formulation on separate gridcells by introducing pressures located at the middle of the cell faces. This makes MFD could deal with any complex grid. To summarize, different from existing multiscale formulations, by combing multiscale mimetic method and EDFM, this method could avoid the complicated unstructured grid generation and computing process. Moreover, as discussed in many works (Skaflestad and Krogstad, 2008; Aarnes et al., 2008, 2014), multiscale mimetic method can deal with highly complex grid systems because of its applicability for the complex unstructured grid. Besides, the excellent local conservation property of MFD method makes multiscale mimetic method produce conservative velocity fields which are necessary for flow simulation. This paper proceeds as follows: we start by introducing the embedded discrete fracture model and its discretization. Next, we describe the multiscale mimetic method briefly and then extend it to deal with the EDFM system. In the numerical experiment section, several 2D and 3D numerical test cases are presented to validate the correctness and effectiveness of the multiscale method. Finally the concluding remarks were given in the final section.
181
Fracture system:
−∇ · (Kf λf · ∇ pf ) = qf −
qmf q − ff δff Vf Vf
(2)
Here, Ki (i = m, f) stands for the permeability tensor; pi, Vi represent pressure, volume of grid cells, respectively; qi = qin + qiw is sink/source term (the subscript w denotes the wetting phase and n denotes the non-wetting phase); λi = λin + λiw is the total mobility where λiα = kriα /μα is the mobility of phase α , which depends on relative permeabilitykriα and viscosity μα ; qmf is transfer flow between matrix and fracture; qff is transfer flow between intersecting fractures; δ mf = 1 if matrix cell contains fracture cells, else δ mf = 0; δ ff = 1 if fracture cell intersects with another fracture cells, else δ ff = 0. The pressure in matrix is assumed to be continuous, therefore, the crossflow term can be written in this form
qmf = −Tmf ( pm − pf )
(3)
Here, Tmf = kmf λmf Amf /d, where Amf is interfacial area, d is equivalent distance between matrix cell and fracture cell, λmf is total mobility decided by upstream weight method, and kmf is harmonically averaged permeability. The crossflow term qff of intersected fractures (Fig. 1) is given by (Karimi-Fard et al., 2004)
qff = Tff pfi − pf j
(4)
, T = k d /d ; Here, Tmf = λff Tfi Tfj /(Tfi + Tfj ), where Tfi = kfi dfi /d i j fj fj fj λff is total mobility decided by upstream weight method; df and kf stand for fracture aperture and permeability respectively; d represents the average normal distance between the centre of fracture segments and the intersections of fractures. The saturation equations are written as: Matrix system
φm
∂ Smw q + ∇ · vmw = qmw + mfw δm f ∂t Vm
vmw = fmw [vm + Km λmn · ∇ pmc ]
(5) (6)
Fracture system
φf
∂Sfw q f fw q + ∇ · v f w = q f w − mfw δm f − δ ∂t Vf Vf f f
v f w = f f w [v f + K f λ f n · ∇ p f c ]
(7) (8)
Here, φ i denotes the porosity, fiw = λiw /λi denotes the fractional flow function, and pic is the capillary pressure. In this work, considering the computational efficiency, we use IMPES strategy to solve the coupled system. IMPES means the flow equations are solved implicitly to obtain the velocity field. Then the velocity field is employed to solve the saturation equations. 2.2. Discretization
2. Embedded discrete fracture modelling 2.1. Mathematical model The flow system for incompressible and isothermal two-phase flow without considering the influence of gravity in fractured media is governed by pressure and saturation equations. The pressure equations are written as Matrix system:
−∇ · (Km λm · ∇ pm ) = qm +
qmf δmf Vm
Matrix grid and fracture grid are constructed to solve the flow equations (Fig. 2). The fracture grid is embedded into matrix grid independently. Therefore, this method can alleviate the gridding complexities. To efficiently handle polyhedral cells, MFD is used as our discretization method. As shown in Fig. 3, Ak is the interface between the matrix cells, nk is the normal to the interface, xik is the vector pointing from cell centroid to the interface centroid. f Let vmi = [vmi1 , vmi2 , · · ·,vmin ]T be the fluxes of interfaces of i , f
(1)
pemi and pmk is the average pressures located at the gridcell centres and interfaces respectively. These quantities are related through a transmissibility matrixTmi , that is
182
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
li1 lj1
lj2
li2
Fig. 1. Schematic of the intersecting fractures in embedded discrete fracture model.
(a) matrix grid
(b) fracture grid
Fig. 2. Schematic of the independent grids in embedded discrete fracture model.
The matrix system can be coupled with the fracture system to form the resulting equation
⎡ nk
Ak Ði
x ik
pei
pfk
⎢CmT ⎢ ⎢ T ⎢Dm ⎢ ⎣0
pej
Ðj
vmf i = Tmi · ei pemi − pmf k
(9)
1]T .
Here, ei = [1, 1, ..., By using the divergence theorem, we can derive the following equation from Eq. (1)
k=1
i
qmi d + qm f δm f
(10)
In this work, the dimension of fractures is reduced to (d-1)dimension. By using implicit difference method, the flow equations in the fracture system are discretized as
Tξ i+ 1 pfi+1 − pfi − Tξ i− 1 pfi − pfi−1 = ffi + qmfi + qffi δffi 2
2
Here, Tξ i± 1 = Kf λf 0.5( ξ 2
0
0
Tmf1 + Tmf2
0
−Tmf1
−Tmf2
0
0
0
0
Tmf1
0
Tf1 − Tmf1 − Tff
Tff
Tmf2
0
Tff
Tf2 − Tmf2 − Tff
⎡ ⎤ ⎡ ⎤ vm 0 ⎢ pm ⎥ ⎢ f m ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ×⎢ ⎢πm ⎥ = ⎢0 ⎥ ⎣ pf1 ⎦ ⎣ ff1 ⎦
Fig. 3. Schematic of mimetic finite difference method.
vmf k =
Dm
0
face centroid cell centroid
n
−Cm
Bm
dfi
i±1 + ξ i )
(11)
, and ffi = Vfi qfi , where Vfi and dfi
represent the volume and aperture of i thfracture element, respectively.
pf2
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(12)
ff2
Where Bm and Cm are block diagonal for which the gridblocks are Tm−1 and an Ne × 1 vector with all entries equal one, where Tmi i denotes transmissibility between cell central and cell boundaries of gridcell i . This denotes that the mimetic method works directly on each gridcell. Ne denotes the number of grids. Tmfi is the matrix of transfer coefficient between fracture cell and matrix cell, Tff is the matrix of transfer coefficient between two intersected fractures. The sought unknowns are the face-velocity vm , the matrix gridcell pressures pm , matrix face-pressure π m , the fracture pressure pfi . Developing an efficient approach for Eq. (12) is a challenging problem in field applications. One reason is the system may contain up to several millions of sought unknowns in field-scale applications. Another reason is the great differences between fracture and matrix properties would deteriorate the value of condition numbers. Apparently traditional upscaling method cannot be used to solve this system because of the high resolution fractures, which have essential influence on whole flow patterns. This creates a mo-
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) matrix grid
183
(b) fracture grid
Fig. 4. Schematic of multiscale grids in multiscale mimetic method.
⎧ ⎨ ωi (x ), x ∈ i ∇ · ψi j = −ω j (x ), x ∈ j ⎩ 0, x∈ / i j
In ij with ψ ij • n = 0 ∀x ∈ ∂ ij . No-flow boundaries are imposed. We should notice that it is generally not convergent as the mesh size tends to zero when K is an anisotropic tensor. A homogenization-based mixed multiscale finite element could be applied to avoid this problem (Arbogast, 2011). The ωi (x) must satisfy ωi (x )dx = 1 so that there is a unit i average flow over the ij between two coarse gridcells. The most simplest option is ωi (x) = 1/|i |, in this work, in order to avoid the computational error resulting from homogenization, the ωi (x) is chosen as
Fig. 5. Schematic description of basis functions.
tivation for the multiscale method, which have an attractive advantage of solving coarse problems while honouring fine-scale information. The introduction of this method is presented in the following section.
ω (x )i =
In this part, we briefly describe the multiscale mimetic framework, an efficient multiscale procedure for fractured porous media. This method is formulated using a coarse-scale grid and a finescale grid. The multiscale mimetic method first superposes a finescale grid on both matrix and fractures as shown in Fig. 3. The petrophysical properties of matrix and fracture are given in the fine-scale grid. The coarse-scale grid is introduced on top of the fine-scale grid as seen in Fig. 4. Each block is composed of a connected collection of fine-scale gridcells. In this work, the fracture gridcells are also decomposed to coarse-scale grid and fine-scale grid. Over this coarse-scale grid, two approximation spaces are defined: each coarse gridcell corresponds to a pressure approximation space while each interface corresponds to a flux approximation space. Next, all pairs of adjacent blocks are detected for the construction of multiscale basis functions. This is a crucial step of multiscale method since the fine-scale features and the interaction between the fractures and the matrix are accounted for with multiscale functions. Consider two coarse-scale grid cells i and j (Fig. 5) and let ij consists of i and j , that is ij = i ∪ ij ∪j .
ij is the interface between i and j . The basis functions ϕ ij for pressure and ψ ij for velocity are constructed by solving following local flow equations (Aarnes et al., 2005)
(13)
⎧ ⎪ σ x / σ ξ d ξ , qdx = 0 ( ) ( ) ⎨
⎪ ⎩q(x )/
i
i
q(ξ )dξ
i
i
(15) qdx = 0
Where σ (x) = trace(K(x))/d. By using MFD, the hybrid system of Eqs. (13)–(15) can be written as (Aarnes et al., 2006)
⎡
3. Multiscale mimetic method for fractured porous media
ψi j = −λt K ∇ ϕi j
(14)
−C
B
⎣C T
0
DT
0
⎤ ⎤ ⎡ ψi j 0 0 ⎦⎣ϕi j ⎦ = ⎣ωi θi j ⎦ π 0 0
D
⎤⎡
(16)
Here B and C are block diagonal for which the gridblocks are Tm−1 and an ne × 1 vector with all entries equal one, Tm is transmissibility matrix, ne is the number of gridcells in the block ij , θ ij = 1 if x ∈ i , θ ij = −1 if x ∈ j . The velocity basis function ψ ij is represented as a vector of fluxes on the set of faces inside the coarse gridcells. Similarly, the pressure basis function ϕ ij is represented as a vector of fine-scale gridcells pressure. When the coarse blocks contain fractured regions, the fractures are represented by EDFM. This model is effective for constructions of basis functions, especially when the fractures are close to each other, which may lead to a poor unstructured grid generation and calculation error. The flowing system is solved for basis functions
−∇ · Km λm · ∇ϕimj − ψ f m δmf = ω (x )i
−∇ · Kf λf · ∇ϕifj + ψ m f δm f + ψ f f δ f f = ω (x )i
(17) (18)
The basis functions, including the fracture basis functions ψ fm and ψ ff , are calculated using a linear system like Eq. (12). The fracture velocity basis function ψ fm is represented as a vector of fluxes on the faces between fracture and matrix, the fracture pressure baf sis function ϕi j is represented as a vector of fine-scale fracture cells
184
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Fracture
(a) Coarse grid cell with partially overlapping fracture
(b)
(c)
(d)
(f)
Fig. 6. Illustration of all four pressure basis functions of a coarse cell. (b) is basis function of bottom edge, (c) is the basis function of right edge, (d) is the basis function of upper edge, (f) is the basis function of left edge.
important factors (such as matrix, fractures) in the construction of the global equation system (Fig. 6). If we construct multiscale basis functions in a larger domain, it means the multiscale basis functions are constructed using oversampling technique. Multiscale basis functions normally only have support in the corresponding coarse gridcells. Oversampling means that we extend the support of the multiscale basis functions (see Fig. 7), that is, we include more fine gridcells when computing the basis function. The use of oversampling is particularly attractive in cases of highly heterogeneous reservoirs. In fractured blocks, the use of oversampling technique is more important than its application in the matrix blocks. By using oversampling technique, we can extend the support of the basis functions and lower the impact of boundary. Specifically, the multiscale basis functions ψ Ej (also called oversampling functions) are constructed in the larger domain KE ⊃K (see Fig. 7), that is, the support of the multiscale basis functions is expanded from K to KE . We then form the actual ψ i by linear combination of ψ Ej ,
Fig. 7. Illustration of oversampled region.
pressures, similar to the matrix basis functions. Here, basis function ψ ij is constructed for matrix-matrix impacts, ψ fm is for the matrix-fracture coupling, ψ ff is for fracture-fracture interactions. Multiscale basis functions are applied to capture the impacts of all
ψi =
e
ci j ψ jE
(19)
j=1
where cij are determined by the corresponding boundary conditions of target block.
2-D problem
3-D problem
e5 I
e5 e4
e1
I e3
e2
e4
e1
e3
Fig. 8. Schematic with intersecting fractures at I.
e2
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
185
No-flow
1. 0
( 0. 2, 0. 52)
( 0. 8, 0. 52)
P=0Mpa
P=1Mpa
( 0. 52, 0. 8)
( 0. 52, 0. 2) 0
1. 0
No-flow
(a) Fractures model
(c) EDFM pressure
(b) Multiscale grid
(d) Multiscale coarse-scale results
(e) Multiscale fine-scale results
Fig. 9. Reference pressure solution and multiscale pressure solutions.
Fig. 10. Velocity and pressure distribution along y = 0.52 L of EDFM and multiscale solutions.
Once the multiscale basis functions are constructed, we then assemble the coarse-scale system. To form the coarse system, we split the basis functions into two parts
ψi j = ψiHj − ψ jiH
(20)
Where
organized as ψ ff and ψ fm . Then the global coarse system reads
⎡
c −Cm
Dcm
0
0
Tmc f 1 + Tmc f 2
0
−Tmc f 1
−Tmc f 2
0
0
Bcm
⎢C c T ⎢ m ⎢ cT ⎢D m ⎢ ⎢ 0 ⎣
Tmc f 1
0
−ψi j (E ), E ∈ j ψ i j ( E ) , E ∈ i j \ j H H ψi j ( E ) = · ψ ji (E ) = 0, E∈ / i 0, E∈ / j (21) Next, we assemble all the velocity basis functions ψiHj as
columns in ψ and assembled the pressure basis function as columns in matrix . Likewise, the fractures basis functions are
⎡
vc
⎤
0
Tmc f 2
⎡
0
⎤
⎢ pc ⎥ ⎢ f c ⎥ ⎢ c⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ×⎢ ⎢π ⎥ = ⎢0 ⎥ ⎣ pc ⎦ ⎣ f c ⎦ f1
f1
pcf2
c ff2
0
0 T fc1
−
Tmc f 1 T fcf
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
0 −
T fcf
T fcf Tf2c − Tmc f 2 − T fcf
(22)
Here Bc = ψ T Bf ψ , Cc = ψ T Cf I, fc = IT qf , and Dc = ψ T Cf J .I prolongates blocks to cells, that is, Iij =1 if blocks number j contains cell number i and Iij =0 otherwise. Similarly, J prolongates block
186
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Permeability field
(b) Multiscale grid
Fig. 12. Multiscale flux error for different coarse meshes.
(c) EDFM pressure
(d) Multiscale solution
Fig. 11. Pressure distribution of EDFM and multiscale solutions.
faces to fine faces, Jij =1 if coarse face j contains the fine face i and Jij =0 otherwise. Once the Eq. (22) is solved, the fine-scale velocity and pressure can be constructed as (Krogstad et al., 2009, 2012)
vf ≈ ψ vc
(23)
pf ≈ I pc + Dλ vc
(24)
The fracture velocity can be obtained by
v ≈ ψ ff vc ff
(25)
¯ 0 /λ ¯ i for block i, λ ¯0 Here Dλ is a diagonal matrix with entries λ i i stands for the mean mobility in the construction of the basis func¯ i stands for the mean mobility in the current time. tions and λ 4. Approximation of the saturation Fig. 13. Total CPU time used by multiscale method and reference method.
The saturation functions are discretized with a finite volume method.
No-flow
( 0. 2, 0. 25)
0
1. 0
( 0. 8, 0. 25)
No-flow
(a) Fractures model
(c) Reference saturation field
P=0Mpa
P=1Mpa
0. 5
1. 0
(b) Multiscale grid
(d) Multiscale saturation field
Fig. 14. Saturation distribution of EDFM and multiscale solutions.
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Here,
Ffk (Sfw ) =
Afk
187
ffw (Sfw )k (vf · nk + Km λfn · ∇ pfc · nk )dA
(28)
In above Eqs (26)–(28), n stands for time step, Ni denotes number of faces, fiw (S) is fractional flow, which is obtained by upstream scheme. Upstream-weighting method is applied to deal with multiple intersecting fractures. We assume that there are NI fracture elements ei intersect at Iin Fig. 8. Here, fw,ei and vf,ei stand for fractional function and flux respectively. n is defined as 0 < n < NI , such that
(a) matching grid
(b) Reference saturation field
vf,ei ≤ 0, 0 < i ≤ N (efflux ) vf,ei ≥ 0, N < i < NI (influx )
(29)
Considering the law of conservation of volumetric and mass, one gets NI
vf,ei = −
i=N+1 NI
N
vf,ei
(30)
i=1
fw,i vf,ei = −
i=N+1
N
fw,I vf,ei = − fw,I
i=1
N
vf,ei
(31)
i=1
Therefore, the fractional function at the intersections is
NI
fw,I = −
v
f i=N+1 w,i f,ei N i=1 f,ei
v
(32)
5. Numerical experiments
(c) Multiscale grid
(d) Multiscale saturation field Several numerical experiments conducted in this work are given in this part. Firstly, the multiscale mimetic method is validated for fractured reservoirs through a comparison with a fine-scale simulation of fractured media. Then a simulation of fractured media with heterogeneous matrix is presented to investigate efficiency and robustness of the multiscale method. The third example is set to validate the multiscale method for multiphase flow in fractured reservoirs. Finally, we extended our multiscale method to 3D fractured model and demonstrated its powerful industrial potential. 5.1. Verification of the multiscale method for fractured reservoirs
(e) Water-cut curve of irregular fractured media Fig. 15. Saturation distribution of DFM and multiscale solutions.
For the matrix system Nm φmi n+1 1 n n (Smwi − Smw ) + ) θ F (Sn+1 ) + (1 − θ )Fmk (Smw i t |mi | k=1 mk mw n = qmwi (Smw i) +
n qm f w ( Sw ) δmf Vmi
(26)
In the first scenario, a 2D homogeneous reservoir was considered, where two orthogonal fractures are located at the centre of the L × L region. The permeability of both fractures satisfy Kf /Km = 105 . We impose no-flow boundary conditions on the top and bottom while the left side and right side are set to constant pressure boundaries p = 1 and p = 0 respectively. The fine-scale reference pressure is obtained with 40 × 40 grid cells. For the simulation of multiscale mimetic-EDFM, 10 × 10 coarse grid cells are employed (Fig. 9(b)). Numerical results of EDFM and multiscale method are presented in Figs. 9 and 10. From the comparison of the results, we can observe that the pressure distribution obtained by reference solution can be reproduced by the fine-scale results of the multiscale method. 5.2. Efficiency analysis of the multiscale method
For the fracture system Nf φfi n+1 n 1 n (Sfwi − Sfwi ) + ) θ F (Sn+1 ) + (1 − θ )Ffk (Sfw t |fi | k=1 fk fw n n qf f w Sfw qm f w ( Sw ) n = qfw j (Sfw ) − δ − δff mf j
Vf j
Vf j
(27)
The second test was presented to explore the efficiency and robustness of the multiscale method for more general fractured reservoirs. In this test, a 2D fractured reservoir model with a heterogeneous permeability field is considered (Fig. 11). Boundary conditions of this model are same as the previous test. The finescale pressure is obtained with 40 × 40 grid cells and for the multiscale method 10 × 10 coarse grid is employed.
188
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Fine-scale matrix grid
(b) Fine-scale fracture grid
(c) Coarse-scale matrix grid
(d) Coarse-scale fracture grid
Fig. 16. Schematic of multiscale grid and fractured model.
In this section, we use relative L2 norm to measure the velocity discrepancies:
ev =
vf − vms 22 vf 22
(33)
where vf denotes reference flux solutions while vms stands for the flux calculated with multiscale mimetic method. Fig. 12 shows the flux L2 error for different coarse grid. It demonstrates that the coarse grid has a significant impact on the computation accuracy. Our multiscale mimetic method generally gives better results if the coarse grid has approximately the unit aspect ratio. Besides, it seems to perform best on the coarse grid with medium number coarse gridcells. Fig. 13 compares the computing time of traditional numerical (MFD) method and multiscale method. It shows that multiscale method cost different time with different coarse grid. Overall, multiscale method is an efficient numerical method. 5.3. Multiphase flow in fracture network This section is set to validate our multiscale method for multiphase flow in fractured reservoirs. First, a heterogeneous 2L × L media embedded with one fracture is considered (Fig. 14). The permeability of fracture satisfies K f /K¯ m = 105 , where K¯ m is the matrix average permeability. No-flow conditions are employed on the upper and lower boundaries. Constant pressure differential was imposed along the horizontal direction. The fine-scale grid contains 40 × 20 cells and for multiscale method 10 × 5 coarse grid cells were employed. Then one-well injection and one-well production model with irregular shape (Fig. 15) is presented. No-flow conditions are used on all boundaries in this test. The matching grid used by discrete fracture model contains 577gridcells. The mul-
tiscale grid contains 40 coarse-scale gridcells and 525 fine-scale gridcells respectively. An injection well located at the bottom left corner while the production well located at the top right corner. Linear relative permeability is employed for flow in porous media. In this study, water cut is defined as fw =Qw /(Qo +Qw ), where Qw is the water flow rate and Qo is the oil flow rate at the outlet. Figs. 14 and 15 shows saturation distributions obtained with MFD method and multiscale mimetic method. The saturation error of Fig. 14 in L1 norm is 0.0981. Both test cases indicate that the finescale and multiscale results are nearly identical, even though the water-cut curve in reference method presents bigger value. This confirms that the multiscale mimetic-EDFM can track the reference solution very closely in the case of multiphase flow. 5.4. Multiphase flow in 3D fractured reservoirs Finally, we extend the multiscale mimetic method to 3D fractured reservoir model and take a crucial step towards the final real-field simulation. In this case, a 3D homogeneous media containing 4 fractures are considered. The fine-scale grid is 10 × 10 × 4 and the multiscale method was run using coarsening 5 × 5 × 2 grid (Fig. 16). The injection well is located at the bottom left corner (i.e.x = 0 m, y = 0 m) and the production is located in the top right corner (i.e.x = 100 m, y = 100 m). The simulation results depicted in Fig. 17 including the saturation profiles after 0.5 PVI, the pressure distributions and the watercut curve of the fractured reservoir model. These simulation results demonstrate the powerful industrial potential of our multiscale method. Based on the effective multiscale method, we can conduct various analyses such as production prediction, risk assessment and develop effective management strategy.
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
(a) Reference saturation maps after 0.5 PVI
(c) Reference pressure distribution
189
(b) Multiscale saturation maps after 0.5 PVI
(d) Multiscale pressure distribution
(e) Water-cut curve of 3D fractured media Fig. 17. Schematic of multiscale grid and fractured model.
6. Concluding remarks In this work we develop a multiscale mimetic-EDFM for multiphase flow in fractured porous media. By using embedded discrete fracture model to construct the multiscale basis functions, we can capture the detailed interactions between the fractures and the background. Moreover, the use of embedded discrete fracture model enables us to avoid the challenges related to unstructured grid which results in a complicated grid generation and computing progress. Our multiscale mimetic-embedded discrete fracture model is under the framework of mimetic finite difference method, the excellent conservation property and robustness of MFD method makes the multiscale method can deal with complex and highly unstructured grids systems (Aarnes et al., 2008, 2014). Some nu-
merical results presented in the paper demonstrate that the multiscale mimetic-embedded discrete fracture is an accurate and efficient method for fractured reservoir. The construction of multiscale basis functions can be completed with parallel computing method to further reduce memory requirements. Several multiscale methods could get an increased accuracy and stability by iterative methods (Krogstad et al., 2009; Hajibeygi et al., 2011, 2008; Kuenze and Lunati, 2011) and this is a good inspiration for us. The design of an iterative method that improves the solution accuracy systematically is our ongoing research. The above lead to the conclusion that the multiscale mimeticembedded discrete fracture method is an important multiscale development for the simulation of fractured reservoirs.
190
Q. Zhang et al. / Advances in Water Resources 107 (2017) 180–190
Acknowledgements The authors gratefully acknowledge support from National Natural Science Foundation of China (51404292), the Fundamental Research Funds for the Central Universities (17CX06007), National Science and Technology Major Project (2016ZX05060-010). References Aarnes, J.E., 2006. On the use of a mixed multiscale finite element method for greater flexibility and increased speed or improved accuracy in reservoir simulation. SIAM J. Multiscale Model. Simul. 2 (3), 421–439. Aarnes, J.E., Kippe, V., Lie, K.A., 2005. Mixed multiscale finite elements and streamline methods for reservoir simulation of large geomodels. Adv. Water Res. 28 (3), 257–271. Aarnes, J.E., Krogstad, S., Lie, K.A., 2006. A hierarchical multiscale method for two-phase flow based upon mixed finite elements and nonuniform coarse grids. SIAM J. Multiscale Model. Simul. 5 (2), 337–363. Aarnes, J.E., Krogstad, S., Lie, K.A., 2008. Multiscale mixed/mimetic methods on corner-point grids. Comput. Geosci. 12 (3), 297–315. Aarnes J.E., Lie K.A., Laptev V., Krogstad S. (2014). Multiscale mixed methods on corner-point grids: mimetic versus mixed subgrid solvers. Akkutlu, I., Efendiev, Y., Vasilyeva, M., 2015. Multiscale model reduction for shale gas transport in fractured media. Comput. Geosci. 1–21. Arbogast, T., 2011. Homogenization-based mixed multiscale finite elements for problems with anisotropy. SIAM J. Multiscale Model. Simul. 9 (2), 624–653. Babuska, I., Caloz, G., Osborn, J.E., 1992. Special finite element methods for a class of second order elliptic problems with rough coefficients. SIAM J. Numer. Anal. 31 (4), 945–981. Babuska, I., Osborn, J.E., 1981. Generalized finite element methods: their performance and their relation to mixed methods. SIAM J. Numer. Anal. 20 (3), 510–536. Chen, Z., Hou, T.Y., 2003. A mixed multiscale finite element method for elliptic problems with oscillating coefficients. Math. Comput. 72 (242), 541–576. Efendiev, Y., Galvis, J., Hou, T.Y., 2013. Generalized multiscale finite element methods (GMsFEM). J. Comput. Phys. 251 (23), 116–135. Efendiev, Y., Ginting, V., Hou, T., Ewing, R., 2006. Accurate multiscale finite element methods for two-phase flow simulations. J. Comput. Phys. 220 (1), 155–174. Efendiev, Y., Lee, S., Li, G., Yao, J., Zhang, N., 2015. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method. GEM - Int. J. Geomath. 6 (2), 1–22. Efendiev, Y.R., Wu, X.H., 2002. Multiscale finite element for problems with highly oscillatory coefficients. Numer. Math. 90 (3), 459–486. Ghorayeb, K., Firoozabadi, A., 20 0 0. Numerical study of natural convection and diffusion in fractured porous media. Spe J. 5 (5), 12–20. Gulbransen, A, Hauge, V.L., Lie, K.A., 2009. A multiscale mixed finite-element method for vuggy and naturally fractured reservoirs. SPE J. 15 (15), 395–403. Hajibeygi, H., Bonfigli, G., Hesse, M.A., Jenny, P., 2008. Iterative multiscale finite-volume method. J. Comput. Phys. 227 (19), 8604–8621. Hajibeygi, H., Jenny, P., 2009. Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media. J. Comput. Phys. 228 (14), 5129–5147. 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.
Hauge, V.L., Aarnes, J.E., 2009. Modeling of two-phase flow in fractured porous media on unstructured non-uniformly coarsened grids. Transp. Porous Media 77 (3), 373–398. Hoteit, H., Firoozabadi, A., 2008. An efficient numerical model for incompressible two-phase flow in fractured media. Adv. Water Res. 31 (6), 891–905. Hou, T.Y., Wu, X.H., 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134 (1), 169–189. Huang, Z.Q., Yao, J., Wang, Y.Y., 2011b. Numerical study on waterflooding development of fractured reservoir with discrete-fracture model. Chinese J Comput Phys 28, 148–156. Huang, Z.Q., Yao, J., Wang, Y.Y., Tao, K., 2011a. Numerical study on two-phase flow through fractured porous media. Sci. China Technol. Sci. 54 (9), 2412–2420. Jenny, P., Lee, S.H., Tchelepi, H.A., 2003. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. J. Comput. Phys. 187 (1), 47–67. Jenny, P., Lee, S.H., Tchelepi, H.A., 2004. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. SIAM J. Multiscale Model. Simul. 3 (1), 50–64. Juanes, R., 2005. A variational multiscale finite element method for multiphase flow in porous media. Finite Elements Anal. Des. 41 (7–8), 763–777. Karimi-Fard, M., Durlofsky, L.J., Aziz, K., 2004. Anefficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Karimi-Fard, M., Firoozabadi, A., 2001. Numerical simulation of water injection in 2D fractured media using discrete fracture model. SPE Reservoir Eval. Eng.. Krogstad, S., Lie, K.A., Nilsen, H., Natvig, J., Skaflestad, B., Aarnes, J., 2009. A multiscale mixed finite element solver for three phase black oil flow. SPE Reservoir Simulation Symposium. Krogstad, S., Lie, K.A., Skaflestad, B., 2012. Mixed multiscale methods for compressible flow. Ecmor Xiii, European Conference on the Mathematics of Oil Recovery. Kuenze, R., Lunati, I., 2011. An adaptive iterative multiscale finite volume method to simulate density- driven flow instabilities. In: AGU Fall Meeting, p. 1476. 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. 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. Matei, Tene, 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. Moinfar, A., Varavei, A., Sepehrnoori, K., Johns, R.T., 2012. Development of a novel and computationally-efficient discrete-fracture model to study ior processes in naturally fractured reservoirs. SPE Improved Oil Recovery Symposium. Natvig, J.R., Skaflestad, B., Bratvedt, F., Bratvedt, K., Lie, K.A., Laptev, V., Khataniar, S.K., 2011. Multiscale mimetic solvers for efficient streamline simulation of fractured reservoirs. SPE J. 16 (4), 880–888. Skaflestad, B., Krogstad, S., 2008. Multiscale/mimetic pressure solvers with near-well grid adaptation. Ecmor. Yan, X., Huang, Z., Yao, J., Li, Y., Fan, D., 2016. An efficient embedded discrete fracture model based on mimetic finite difference method. J. Petrol. Sci. Eng. 145, 11–21. Zhang, Q., Yao, J., Huang, Z., 2016. An efficient multiscale mixed finite element method for modelling flow in discrete fractured reservoirs. ECMOR XV, European Conference on the Mathematics of Oil Recovery. Zhou, F., Shi, A., Wang, X., 2014. An efficient finite difference model for multiphase flow in fractured reservoirs. Petrol. Explorat. Develop. 41 (2), 262–266.