Polytopic discontinuous Galerkin methods for the numerical modelling of flow in porous media with networks of intersecting fractures

https://doi.org/10.1016/j.camwa.2021.08.015Get rights and content

Abstract

We present a numerical approximation of Darcy's flow through a porous medium that incorporates networks of fractures with non empty intersection. Our scheme employs PolyDG methods, i.e. discontinuous Galerkin methods on general polygonal and polyhedral (polytopic, for short) grids, featuring elements with edges/faces that may be in arbitrary number (potentially unlimited) and whose measure may be arbitrarily small. Our approach is then very well suited to tame the geometrical complexity featured by most of applications in the computational geoscience field. From the modelling point of view, we adopt a reduction strategy that treats fractures as manifolds of codimension one and we employ the primal version of Darcy's law to describe the flow in both the bulk and in the fracture network. In addition, some physically consistent conditions couple the two problems, allowing for jump of pressure at their interface, and they as well prescribe the behaviour of the fluid along the intersections, imposing pressure continuity and flux conservation. Both the bulk and fracture discretizations are obtained employing the Symmetric Interior Penalty DG method extended to the polytopic setting. The key instrument to obtain a polyDG approximation of the problem in the fracture network is the generalization of the concepts of jump and average at the intersection, so that the contribution from all the fractures is taken into account. We prove the well-posedness of the discrete formulation and perform an error analysis obtaining a priori hp-error estimates. All our theoretical results are validated performing preliminary numerical tests with known analytical solution.

Introduction

This work is concerned with the simulation of Darcean flows through porous media that incorporate networks of fractures with non empty intersection. The focus is on the development and analysis of a numerical approximation that employs PolyDG methods, i.e. discontinuous Galerkin methods on general polygonal and polyhedral (polytopic, for short) grids. In the past decades, increased attention has been given to the efficient implementation of numerical methods for fractured reservoir simulations. The analysis and prediction of the flow is indeed fundamental in many environmental and energy engineering applications, which include petroleum extraction, CO2 storage in depleted oil fields, isolation of radioactive waste and geothermal energy production, for example. In all the aforementioned applications, fractures severely affect the flow, since they can act as barriers for the fluid (when they are filled with low permeable material), or as conduits (when they have higher permeability than the surrounding medium). Moreover, in many cases the geometry of the fault system can be highly intricate, featuring thousands of fractures, which may also intersect with small angles or be nearly coincident [30].

A first step towards a reduction in the complexity of the simulation is usually taken in the conceptual modelling of the flow: since fractures usually present a small width-to-length ratio, as well as a small relative size with respect to the domain, a popular choice consists in treating them as manifolds of codimension one. The development of this kind of reduced models has been addressed for single-phase flows in several works, for example in [1], [2], [31], [34]. Our main reference will be the model described in [34], which considers for simplicity the case of a single fracture that is not immersed in the bulk domain. In the model, the flow in the bulk is governed by Darcy's law, whereas a suitable dimensionally reduced version of the law is formulated on the surface modelling the fracture. Moreover, the exchange of fluid between the fracture and the porous medium is described via some physically consistent coupling conditions. We also remark that both low and large permeable fractures can be handled.

Although the use of dimensionally reduced models avoids the requirement of extremely refined grids inside the fracture domains, realistic simulations still call for high mesh resolution in these areas, so that all the small geometrical features can be captured without resorting to low-quality elements in the classical sense. The mesh generation process within a classical finite element approach can then represent a bottleneck in the whole simulation, especially in 3D, as only computational grids composed by tetrahedral/hexahedral/prismatic elements are supported. The same issue can be encountered in many application areas, ranging from fluid-structure interaction, to wave propagation problems, to name a few. This has motivated a huge effort in the past years in the design of numerical methods supporting meshes made of general polytopic elements. A huge reduction on the computational cost may be achieved by resorting to hybrid mesh generation techniques, for example: first a (possibly structured) grid is generated independently of geometric features (e.g., fractures); secondly the elements are cut according to the required pattern. It follows that the final mesh contains arbitrarily shaped elements in the surrounding of such features and is regular far from them. In addition to the simplicity of this procedure, polytopic meshes present, on average, a much lower number of elements, even on relatively simple geometries, without committing a variational crime [10], [11].

Within this framework, many numerical methods have been developed on top of polytopic meshes in the context of flows in fractured porous media. In particular, we mention [9], [30], where a mixed approximation based on Mimetic Finite Differences was applied; the works [14], [15], where virtual elements were employed to deal with flows in Discrete Fracture Networks, and [23], which uses the Hybrid High-Order method. We also remark that an important alternative is given by the use of non-conforming discretizations, in which the bulk grid can be chosen fairly regular independently of the fractures, since these are considered immersed in the geometry. We refer to [25], [27], [33] for the use of the eXtended Finite Element Method and to [18] for the Cut Finite Element Method.

In [7] we presented an approximation of the coupled bulk-fracture problem that employs PolyDG methods. The inherited flexibility of DG methods in handling arbitrarily shaped, non-necessarily matching, grids and elementwise variable polynomial orders represents, in fact, the ideal setting to handle such kind of problems that typically feature a high-level of geometrical complexity. In particular, since they employ local polynomial spaces defined elementwise without any continuity constraint, DG methods feature a high-level of intrinsic parallelism. Furthermore, the lack of continuity between neighbouring elements allows for the employment of extremely broad families of meshes containing elements with edges/faces that may be in arbitrary number (potentially unlimited) and whose measure may be arbitrarily small; cf. [6] for a comprehensive review on PolyDG methods on polyhedral grids for geophysical applications, including seismic wave propagation and fractured reservoir simulations. The geometric flexibility highlighted so far is not the only motivation to employ such techniques in the context of fractured porous media. A more physically motivated argument is given by the discontinuous nature of the solution at the matrix-fractures interface, which can be intrinsically captured in the choice of the discrete spaces. In addition, the bulk and fractures coupling can be easily reformulated by means of the jump and average operators, which are a fundamental tool in the development of DG formulations, and thus naturally incorporated into the variational formulation. Finally, the abstract setting proposed in [13], based on the flux formulation, allows for the introduction of a unified framework where primal or mixed formulations can be chosen independently for both bulk and fractures, depending on the application at hand and on the quantities of interest. We refer to [8] for further details on the unified analysis and to [7] for a focus on the primal-primal framework. In both cases, our analysis was carried out in the simplified setting of a single, non-immersed fracture. The purpose of the present work is to extend our formulation to networks of intersecting fractures. For simplicity, we consider, as in [7], the primal-primal setting, so that we can mainly focus on handling the intersections. To this aim, we supplement the mathematical model [34] with some suitable physical conditions at the intersections, prescribing the behaviour of the fluid. Following [15], [17], [30], we impose that:

  • pressure between fractures is continuous along the intersections;

  • flux is conserved, so that no exchange of fluid between bulk and fracture network takes place along the intersections.

We mention that more general conditions, where the angle between fractures is taken into account and jumps of pressure across the intersection are allowed, may be imposed. Some examples can be found in [16], [24], [28], [32], [35]. We also mention that the analysis of the mixed-mixed setting in the case of a totally immersed network of fractures has been addressed in [30].

From the DG-discretization point of view, the key instrument for dealing with intersections is the generalization of the concepts of jump and average. If we assume that the fracture network may be approximated by the union of NΓ fractures γk, each of which is a one co-dimensional manifold with zero curvature, i.e. Γ=k=1NΓγk, the intersections correspond to lines when d=3 and to points when d=2. Let us focus for simplicity on the case d=3, see Fig. 1 for an example. Here, the intersection line is denoted by I and each fracture γk, k=1,2,3,4, is characterised by the outward normal vector τk at the intersection, which belongs to the plane containing the fracture. In order to describe the pressure field in all the network, we employ the global variable pΓ=(pΓ1,,pΓNΓ), defined in a suitable product space of all the local fracture spaces. Our aim is to introduce some operators that are able to capture the behaviour of the function pΓ across the intersection line, taking into account the contribution from all the fractures, similarly to how classic jump and average operators [13] describe the discontinuity of a piecewise-continuous function across elemental interfaces. The main difference with respect to the standard case is that the normal vectors, contained into the definition of the operators, are not aligned. This is related to the linear DG approximation of elliptic PDEs on surfaces presented in [26], then extended to high order in [5]. Here, the surface is approximated by a piecewise linear surface composed of planar triangles, so that a new definition of jump and average operators is needed, to take into account the fact that the outward normal vectors of two neighbouring triangles are not, in general, opposite. Our definition is a further generalization, since it considers the intersection of an arbitrary number of planar surfaces. Using the newly defined jump and average operators we are able to define a DG approximation for the problem in the bulk combined with a DG approximation for the problem in the fracture network, where the conditions at the intersection are imposed “in the spirit of DG methods ”. In particular, this means that continuity is enforced penalizing the jump of the pressure (after a suitable definition of the penalization coefficient at the intersection), while balance of fluxes is imposed naturally, similarly to how homogeneous Neumann boundary conditions are usually enforced. Both the bulk and fracture discretizations are obtained employing the SIPDG method extended to the polytopic setting.

The rest of the paper is structured as follows. In Section 2, we present the mathematical model. In Section 3, we introduce the weak formulation of the problem and prove its well-posedness. Section 4 contains the polyDG discretization of the coupled system based on the new definition of jump and average operators at intersections, which is also introduced in Section 4. Finally, Sections 5 and 6 enclose the stability and error analysis of the discrete method. We conclude with Section 7, where we present some preliminary numerical experiments with known analytical solution, so that we are able to verify the obtained convergence rates.

Section snippets

Mathematical model

We consider the domain ΩRd, with d=2,3, representing the porous medium. We assume that the fracture network may be approximated by a collection of one co-dimensional planar manifolds ΓRd1, adopting the reduced model introduced in [34] and extended to fracture networks in [15], [17], [30].

In particular, we consider Γ to be the union of NΓ fractures γk,Γ=k=1NΓγk, with every γk being an open, bounded, connected, planar (d1)-dimensional orientable manifold. Each γk is, in fact, the

Weak formulation

In this section we introduce the weak formulation of the model problem (3)-(5)-(6)-(8a), (8b) and prove its well-posedness.

For the sake of simplicity we will assume that homogeneous Dirichlet boundary conditions are imposed for both the bulk and fracture problems, i.e., gD=0 and gΓ=0. The extension to the general non-homogeneous case is straightforward.

First, we introduce the functional setting. We will employ the following notation. For an open, bounded domain DRd, d=2,3, we will denote by Hs(

PolyDG discretization

In this section we present a numerical discretization for the coupled bulk-network problem that is based on DG methods on polytopic grids. In particular, we discretize both the bulk and fracture network problems in primal form, employing the Symmetric Interior Penalty DG method [12], [38]. The key idea to obtain a DG discretization will be the generalization of the concepts of jump and average at the intersection point/line, so that we will be able to impose the conditions at the intersection

Well-posedness of the discrete formulation

In this section, we prove that formulation (20) is well-posed.

We recall that, for simplicity in the analysis, we are assuming the permeability tensors ν and ΓνΓτ to be piecewise constant. We will employ the following notation ν¯E=|ν|E|22 and ν¯Fτ=|ΓνΓτ|F|22, where ||2 denotes the l2-norm.

In order to work in a polytopic framework, we need to introduce some technical tools as in [4], [19], [20], [21], [22]. The first tool consists in trace inverse estimates, so that the norm of a polynomial on

Error analysis

In this section, we derive a-priori error estimates for the discrete problem (20). To this aim, in the following, we summarize the results contained in [4], [19], [20], [21], [22], where standard hp-approximation bounds on simplices are extended to arbitrary polytopic elements. These results are indeed the basic tool for the error analysis of DG methods.

Numerical experiments

In this section we present several numerical examples, with increasing complexity, in order to validate the performance of our method. In the experiments the analytical solution is known, so that we are able to verify the convergence rates obtained in Theorem 6.6. We point out that choice of the model coefficients is here made only with the aim of testing the effectiveness of the numerical method and it does not intend to have any physical meaning. We remark that, in all the presented test

References (38)

  • P.F. Antonietti et al.

    High–Order Discontinuous Galerkin Methods on Polyhedral Grids for Geophysical Applications: Seismic Wave Propagation and Fractured Reservoir Simulations

    (2021)
  • P.F. Antonietti et al.

    Discontinuous Galerkin approximation of flows in fractured porous media on polytopic grids

    SIAM J. Sci. Comput.

    (2019)
  • P.F. Antonietti et al.

    Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids

    Math. Eng.

    (2020)
  • P.F. Antonietti et al.

    Mimetic finite difference approximation of flows in fractured porous media

    ESAIM: Math. Model. Numer. Anal.

    (2016)
  • P.F. Antonietti et al.

    hp-version composite discontinuous Galerkin methods for elliptic problems on complicated domains

    SIAM J. Sci. Comput.

    (2013)
  • P.F. Antonietti et al.

    Domain decomposition preconditioners for discontinuous Galerkin methods for elliptic problems on complicated domains

    J. Sci. Comput.

    (2014)
  • D.N. Arnold

    An interior penalty finite element method with discontinuous elements

    SIAM J. Numer. Anal.

    (1982)
  • D.N. Arnold et al.

    Unified analysis of discontinuous Galerkin methods for elliptic problems

    SIAM J. Numer. Anal.

    (2001/2002)
  • W.M. Boon et al.

    Robust discretization of flow in fractured porous media

    SIAM J. Numer. Anal.

    (2018)
  • Cited by (3)

    P.F.A. and M.V. have been partially supported by PRIN Italian research grant n. 201744KLJL funded by MIUR and by INdAM-GNCS.

    View full text