Weak Galerkin finite element methods for Darcy flow: Anisotropy and heterogeneity

Abstract

This paper presents a family of weak Galerkin finite element methods (WGFEMs) for Darcy flow computation. The WGFEMs are new numerical methods that rely on the novel concept of discrete weak gradients. The WGFEMs solve for pressure unknowns both in element interiors and on the mesh skeleton. The numerical velocity is then obtained from the discrete weak gradient of the numerical pressure. The new methods are quite different than many existing numerical methods in that they are locally conservative by design, the resulting discrete linear systems are symmetric and positive-definite, and there is no need for tuning problem-dependent penalty factors. We test the WGFEMs on benchmark problems to demonstrate the strong potential of these new methods in handling strong anisotropy and heterogeneity in Darcy flow.

Introduction

For convenience of presentation, we concentrate on the two-dimensional Darcy equation, which is usually formulated as

$$\{\begin{array}{cc}\mathrm{\nabla }\cdot (-\mathbf{K}\mathrm{\nabla }p)\equiv \mathrm{\nabla }\cdot \mathbf{u}=f,\hfill & \mathbf{x}\in \Omega ,\hfill \\ p=p_D,\hfill & \mathbf{x}\in {\Gamma }^D,\hfill \\ \mathbf{u}\cdot \mathbf{n}=u_N,\hfill & \mathbf{x}\in {\Gamma }^N,\hfill \end{array}$$

where $\Omega \subset {\mathbb{R}}^2$ is a bounded polygonal domain, p the unknown pressure, u the Darcy velocity, K a permeability tensor that is uniformly symmetric positive-definite, f a source term, $p_D,u_N$ are respectively Dirichlet and Neumann boundary data, n the unit outward normal vector on ∂Ω, which has a nonoverlapping decomposition ${\Gamma }^D\cup {\Gamma }^N$. When ${\Gamma }^D\ne \varnothing$, the problem has a unique solution.

Here is a list of commonly used finite element methods for solving the Darcy equation, or more generally, elliptic boundary value problems.

• Continuous Galerkin (CG) FEMs [7]. There are different opinions on whether CGFEMs are not locally conservative [21]. CGFEMs generally have smaller numbers of unknowns. CGFEMs assume continuous (polynomial) approximants and are not flexible for certain problems that have solutions with low regularity.

• Discontinuous Galerkin (DG) FEMs [2] (and the references therein). DGFEMs utilize discontinuous shape functions that exhibit great flexibility in handling complicated domain geometry and problems with low regularity. DGFEMs are parallelizable by design. But for the symmetric formulations, one has to choose sufficiently large penalty factors to ensure stability of the numerical methods [18]. Choosing appropriate penalty factors is problem-dependent and hence a drawback in applications of DGFEMs.

• DGFEMs with weak overpenalization (WOP) [6]. The inconvenience of choosing penalty factors motivates the weakly over-penalized DGFEMs, which avoid choosing penalty factors by means of weak over-penalization. However, these methods result in large condition numbers, even though simple block-diagonal preconditioners are available [6].

• Mixed FEMs [9]. The mixed methods solve for the pressure and Darcy velocity simultaneously. But it is well known that the mixed FEMs result in saddle-point problems, whose solvers are more complicated than those for definite systems [5].

• Hybridized FEMs [13] (and the references therein). Hybridization introduces new unknowns (Lagrangian multipliers) and relax the continuity constraints across element interfaces. Hybridization proliferates (conceptually) the numbers of unknowns, although degrees of freedom are actually smaller after certain variables are eliminated.

The weak Galerkin finite element methods introduced in [34] adopt a completely different approach. They rely on novel concepts such as weak functions, the weak gradient operator, discrete weak functions, and discrete weak gradients. As is well known, the variational formulation of a second order elliptic problem relies on the duality of the classical gradient operator. Locality of a finite element space implies relative independence of elementwise shape functions and local conservation. However, shape functions in element interiors could be related to their values on the element interfaces through integration by parts. The weak gradient operator characterizes exactly this connection, see Eq. (3) in Section 2. Discrete weak gradient operators inherit the connection for finite dimensional (Galerkin type) polynomial approximation subspaces, see Eq. (5) in Section 2.

The Darcy's law has its fundamental importance in flow and transport in porous media [12], [15], [16], [17], [19], [23], [24], [27]. Although Darcy's law can be essentially treated as an elliptic problem, the quantities of interest are pressure and velocity (flux) and the focus is placed on flow features. In this paper, we investigate the performance of the novel WGFEMs on handling the numerical difficulties brought up by anisotropy and heterogeneity in permeability.

The rest of this paper is organized as follows. In Section 2, we explain a line of new concepts such as weak functions, weak gradients, discrete weak functions, and discrete weak gradients. Based on these, we develop weak Galerkin finite element methods for Darcy flow. Section 3 addresses implementation issues of the WGFEMs. Section 4 briefly reviews other existing numerical methods for Darcy flow. Section 5 presents numerical simulations of the WGFEMs on several benchmark problems for Darcy flow. Section 6 concludes the paper with some remarks.