Polytopic discontinuous Galerkin methods for the numerical modelling of flow in porous media with networks of intersecting fractures☆
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.
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 fractures , each of which is a one co-dimensional manifold with zero curvature, i.e. , the intersections correspond to lines when and to points when . Let us focus for simplicity on the case , see Fig. 1 for an example. Here, the intersection line is denoted by and each fracture , , is characterised by the outward normal vector 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 , 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 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 , with , representing the porous medium. We assume that the fracture network may be approximated by a collection of one co-dimensional planar manifolds , 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 fractures , with every being an open, bounded, connected, planar -dimensional orientable manifold. Each 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., and . 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 , , we will denote by
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 and , where denotes the -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)
- et al.
The virtual element method for discrete fracture network simulations
Comput. Methods Appl. Mech. Eng.
(2014) - et al.
A globally conforming method for solving flow in discrete fracture networks using the virtual element method
Finite Elem. Anal. Des.
(2016) - et al.
Cut finite elements for convection in fractured domains
Comput. Fluids
(2019) - et al.
Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations
J. Comput. Phys.
(2019) - et al.
A numerical method for two-phase flow in fractured porous media with non-matching grids
Adv. Water Resour.
(2013) - et al.
Modeling fractures as interfaces for flow and transport in porous media
- et al.
Domain decomposition for some transmission problems in flow in porous media
- et al.
Asymptotic and numerical modelling of flows in fractured porous media
ESAIM: Math. Model. Numer. Anal.
(2009) - et al.
Review of Discontinuous Galerkin Finite Element Methods for Partial Differential Equations on Complicated Domains
(2016) - et al.
High order discontinuous Galerkin methods for elliptic problems on surfaces
SIAM J. Numer. Anal.
(2015)
High–Order Discontinuous Galerkin Methods on Polyhedral Grids for Geophysical Applications: Seismic Wave Propagation and Fractured Reservoir Simulations
Discontinuous Galerkin approximation of flows in fractured porous media on polytopic grids
SIAM J. Sci. Comput.
Unified analysis of discontinuous Galerkin approximations of flows in fractured porous media on polygonal and polyhedral grids
Math. Eng.
Mimetic finite difference approximation of flows in fractured porous media
ESAIM: Math. Model. Numer. Anal.
hp-version composite discontinuous Galerkin methods for elliptic problems on complicated domains
SIAM J. Sci. Comput.
Domain decomposition preconditioners for discontinuous Galerkin methods for elliptic problems on complicated domains
J. Sci. Comput.
An interior penalty finite element method with discontinuous elements
SIAM J. Numer. Anal.
Unified analysis of discontinuous Galerkin methods for elliptic problems
SIAM J. Numer. Anal.
Robust discretization of flow in fractured porous media
SIAM J. Numer. Anal.
Cited by (3)
A Stokes–Darcy–Darcy model and its discontinuous Galerkin method on polytopic grids
2024, Journal of Computational PhysicsDiscontinuous Galerkin method for hybrid-dimensional fracture models of two-phase flow
2023, Journal of Computational PhysicsAdaptive discontinuous Galerkin finite element methods for the Allen-Cahn equation on polygonal meshes
2024, Numerical Algorithms
- ☆
P.F.A. and M.V. have been partially supported by PRIN Italian research grant n. 201744KLJL funded by MIUR and by INdAM-GNCS.