An efficient multi-point flux approximation method for Discrete Fracture–Matrix simulations

https://doi.org/10.1016/j.jcp.2012.01.023Get rights and content

Abstract

We consider a control volume discretization with a multi-point flux approximation to model Discrete Fracture–Matrix systems for anisotropic and fractured porous media in two and three spatial dimensions. Inspired by a recently introduced approach based on a two-point flux approximation, we explicitly account for the fractures by representing them as hybrid cells between the matrix cells. As well as simplifying the grid generation, our hybrid approach excludes small cells in the intersection of the fractures and hence avoids severe time-step restrictions associated with small cells. Excluding the small cells also reduces the condition number of the discretization matrix. For examples involving realistic anisotropy ratios in the permeability, numerical results show significant improvement compared to existing methods based on two-point flux approximations. We also investigate the hybrid method by studying the convergence rates for different apertures and fracture/matrix permeability ratios. Finally, the effect of removing the cells in the intersections of the fractures are studied. Together, these examples demonstrate the efficiency, flexibility and robustness of our new approach.

Introduction

In response to external stresses, systems of discontinuities, i.e., joints, fissures and faults, occur in geological porous media. We term all these systems of discontinuities as fractures. Fractures occur on different scales, with different geometries, and may behave as either channels or barriers for the fluid flow. The contrast in the permeability between the fractures and the porous medium frequently spans orders of magnitudes; therefore, fractures yield a significant influence on the fluid flow in many situations. Applications where fractures play an important role in the fluid flow span all engineered geological flows, including energy production and waste disposal. The geometric complexity of the fractures, and the multiple length and time scales involved, has led to the development of a variety of methods targeting fluid flow in fractured porous media.

The first, and most commonly used approach, is to model the fluid flow through effective parameters defined on equivalent continua, either as single, dual or multi continuum models [1], [2]. This approach incorporates the effect of micro-scale variations without resorting to high spatial and temporal resolution. With this concept the interaction between the continua is assumed to be a steady-state process and represented through exchange terms, often approximated to be proportional to the difference between average pressure in the fracture and the matrix [3]. For iso-thermal and single phase flow this approximation may give satisfactory results. However, for problems involving transport of for example, heat, contaminants, or multiple phases, the difficulties associated with correctly characterizing the exchange terms reduce the applicability of the concept. Remedies include improved transfer function with higher order terms [4], [5], and sub-gridding of the matrix as done in the Multiple INteracting Continua (MINC) approach [6], [7].

Discrete Fracture–Matrix (DFM) models attempt to overcome the deficiencies of multi continuum models by explicitly resolving the fractures [8]. This is appropriate in situations where a small number of fractures dominate the system. Most DFM models are based on gridding the fractures explicitly, either with a hybrid formulation where fractures are modeled as lines in two dimensions and planes in three dimensions [9], [10] or with an equi-dimensional formulation where fractures have the same dimension as the matrix [11], [12], [13]. In these approaches, the fracture and the matrix grid coincide, and an unstructured grid is needed in order to honor the explicit fracture geometry. Due to the complex geometry of fracture networks, grid generation is considered a challenge for DFM models [14]. An alternative approach is to embed fractures into a coarse structured grid, as done in [15], where fractures are modelled through transport indices. A recent comparison of some of the models used for fractured porous media can be found in [16]. In [16] the authors compare an embedded DFM model, a conforming DFM model and a dual continuum model for different relevant scenarios.

For many realistic applications, neither a multi continuum model nor a DFM model are feasible, as a DFM model accounting for all the fractures is computationally too demanding and no effective parameters can be defined on relevant scales for the continuum model. Moreover, knowledge about the fracture network is limited, so multiple DFM realizations are needed in order to obtain statistically significant results. One approach to overcome this deficiency is to use DFM models to compute effective parameters for a MINC type model [17]. In this way a link is established between the underlying fracture network and the continua. Another approach is based on combining the DFM model with the continuum model. While long fractures are included as discrete elements, short fractures are included as effective parameters on a continuum scale. Although such scale separation does not strictly exist, a combined approach may often provide reasonable approximation to the fluid flow [18]. An example of a combined approach is given in [15], where a hierarchical classification is used for the fractures. As the fracture orientation in general is anisotropic, the effective permeability will consequently be anisotropic. A space discretization handling anisotropy in the permeability and discrete fracture networks is therefore necessary in order to use the combined approach.

Herein, we introduce a DFM model valid for anisotropic permeabilities on general unstructured grids. We will assume that an effective permeability tensor exists for the small-scale fractures, and focus on the discretization of the DFM model by introducing a control volume finite difference (CVFD) method with a multi-point flux approximation (MPFA). Inspired by ideas presented in [19] for a two-point flux approximation (TPFA), we consider a hybrid formulation where the fractures are lower dimensional in the grid, but equi-dimensional in the computational domain. Our new approach combines the simplicity and efficiency of the method introduced in [19], while providing a consistent discretization for anisotropic permeabilities by using a MPFA scheme. Alternative mass conservative methods used for DFM are control volume finite elements [9], [20] and discontinuous Galerkin mixed finite elements [21]. Our emphasis is on CVFD methods as these methods are used almost exclusively in commercial simulators.

We start in Section 2 by presenting the equations used in this work, before the hybrid grid is presented in Section 3. The main part of the article is in Section 4, where the CVFD discretization is described. A short discussion on how the upwind scheme can be modified to improve the accuracy in the intersection of fractures is included in Section 5, before numerical examples are shown in Section 6 and conclusions based on our preliminary experience with the method are drawn in Section 7.

Section snippets

Model equations

For the derivation of the MPFA method for the DFM model we consider single phase flow. The method derived herein is, however, also applicable to multiphase flow problems.

Single phase flow in porous media is described by the following equation·v=f,xΩ,where f is the source term and v the velocity given by Darcy’s lawv=-Kufor a permeability tensor K and pressure u. In cases where rapid fluid flow occurs in the fractures, Darcy’s law may no longer be valid. Forchheimer’s law corrects for high

The hybrid grid

In order to obtain a CVFD method for DFM systems, Karimi-Fard et.al. [19] presents a hybrid approach where the fractures are lower-dimensional in the geometric grid, but expanded in the computational domain. As well as simplifying the grid generation, the hybrid approach excludes the intermediate cell in the intersection of the fractures, yielding a more robust and efficient discretization of the DFM system.

The first step in creating the hybrid grid is to partition the domain with the

Hybrid CVFD methods for anisotropic media

In this section we formulate CVFD methods valid for the hybrid grid. For Matrix–Matrix connections the standard methods will apply, but adjustment is needed for both Fracture–Matrix and Fracture–Fracture connections. We present the new MPFA method as well as a short description of the TPFA based method as presented in [19].

In order to obtain a CVFD formulation of Eq. (1), we integrate the equation over each cell i and use Gauss’ theorem, obtainingΩifdσ=Ωi-·Kudσ=-Ωin·KudS,where n is the

Upwind schemes for fracture intersections

In this section we describe a modified first order upwind scheme which implicitly accounts for the intermediate cell. The derivation is done only for the linear transport equation as the derivation for the time-of-flight equation is similar.

Integrating the transport equation in (3) over a domain Ω and applying Gauss’ theorem for the second part givesΩϕctdσ+Ω(n·v)cdS=Ωfdσ.To solve (17) using an implicit first order upwind scheme, the discrete version of (17) is written in terms of its

Numerical testing of the hybrid method

To show various aspects of the hybrid method, we present investigations of three test cases. The first case investigates when a hybrid method can be used by studying the convergence of pressure and flux for different apertures and permeabilities in a two dimensional domain with a single fracture. The effect of removing the intermediate cell in the intersection of the fractures is then studied in the second case for an idealized domain with two crossing fractures. In the third case, the TPFA and

Conclusions

In this work we have introduced a method for discretizing Fracture–Matrix systems based on a multi-point flux approximation. Our investigation leads to three main conclusions:

  • 1.

    Multi-point flux approximation methods can successfully be developed for Discrete Fracture–Matrix systems where the fractures are treated as lower-dimensional objects. Key components in the development include expansion of fractures and avoiding explicit representation of fracture intersections.

  • 2.

    Hybrid Discrete

Acknowledgement

The authors thank Brad Mallison at Chevron Energy Technology Company for providing the grid used in the last example and Sintef Applied Mathematics for providing the MATLAB Reservoir Simulation Toolbox (MRST). This work was funded by the Research Council of Norway (190761/s60).

References (26)

  • K. Pruess et al.

    On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal reservoirs

    J. Geophys. Res.

    (1982)
  • K. Pruess et al.

    A practical method for modeling fluid and heat flow in fractured porous media

    (1985)
  • P. Dietrich et al.

    Flow and Transport in Fractured Porous Media

    (2005)
  • Cited by (262)

    View all citing articles on Scopus
    View full text