Enriched multi-point flux approximation for general grids

Abstract

It is well known that the two-point flux approximation, a numerical scheme used in most commercial reservoir simulators, has O(1) error when grids are not K-orthogonal. In the last decade, the multi-point flux approximations have been developed as a remedy. However, non-physical oscillations can appear when the anisotropy is really strong. We found out the oscillations are closely related to the poor approximation of pressure gradient in the flux computation.

In this paper, we propose the control volume enriched multi-point flux approximation (EMPFA) for general diffusion problems on polygonal and polyhedral meshes. Non-physical oscillations are not observed for realistic and strongly anisotropic heterogeneous material properties described by a full tensor. Exact linear solutions are recovered for grids with non-planar interfaces, and a first and second order convergence are achieved for the flux and scalar unknowns, respectively.

Introduction

Many physical processes, for example, the Darcy flow in porous media and the heat transfer in Basin modeling, can be modeled by a diffusion equation in the form

$$-\nabla ·\mathbf{K}\nabla p=f\text{.}$$

The diffusion coefficient K is often a space dependent full tensor because of the strongly anisotropic and inhomogeneous media, which imposes a great challenge for numerical schemes, especially when the principal directions of K are not aligned with the grids. This is usually the case for unstructured meshes. In the following, we will state everything in the context of reservoir simulation, simply because the initial motive is to develop a more robust discretization for simulating multi-phase flows in porous media, although the proposed scheme can apply to general diffusion problems described by Eq. (1). In this paper, we assume K is a piecewise constant tensor. The scalar unknown p will be called pressure for convenience.

There has been extensive research on developing numerical schemes for Eq. (1) on general grids. Some highly desirable properties of the discretization beyond the classical stability and accuracy, include local mass conservation, harmonic average transmissibilities, discrete maximum principle, and cost-efficiency. For historic reasons, the two-point flux approximation (TPFA), a scheme that is inconsistent when the grids are not K-orthogonal, is still used in most commercial reservoir simulators. Roughly speaking, K-orthogonality means that the flux −K∇p · n can be approximated by a factor of the pressure difference from two neighbouring control volumes (see (8) for one example). Unfortunately, K-orthogonality is often lost in the grids that honor the geologic features such as the sloping faults and channels, and other important features such as nearly horizontal wells [18]. Recently, Wu and Parashkevov studied the effect of non-orthogonality error of deviated grids on the flow solutions from the two-point flux, control volume method [35]. They concluded that for most practical cases the errors in horizontal flow are relatively small and the errors in vertical flow can be rather significant.

To obtain a consistent scheme, different groups of researchers have independently developed the so-called control volume multi-point flux approximation (MPFA), such as Aavatsmark et al [3], Edwards and Rogers [15], Verma and Aziz [34]. Since their pioneer work, a whole class multi-point flux methods have been developed [2], [19], [26]. In Ref. [22], the authors implement the multi-point flux methods on a hexahedral grid. More recently, Klausen and Russell [19] compare the multi-point flux approximation with several other locally conservative schemes which handle discontinuous coefficients, such as the classical mixed finite element method and the support operators method (also called mimetic finite difference method). Klausen and Winther also proved the convergence of the multi-point flux approximations on smooth quadrilateral grids [20].

In most cases, good results are obtained by the multi-point flux methods, but non-physical oscillations can appear for strongly anisotropic fields because of non-monotonicity of the matrix (see, e.g. [24], [27]). This means that the numerical solution does not satisfy the maximum principle, although the analytical solution does. In order to improve the monotonicity property, Nordbotten and Eigestad recently developed a MPFA Z-method for quadrilateral grids by choosing a different stencil when computing the flux [29]. But it is not clear how to extend their method to unstructured polyhedral meshes. Some other very recent research on this issue includes [23], [28], [17], [4].

The mimetic finite difference schemes proposed by Shashkov et al. [33] are very promising in terms of dealing with highly distorted grids and heterogeneous media. The convergence and super convergence [6], [7] are also established for smooth problems on smooth meshes by rewriting the scheme into the form of mixed finite element methods. (A more general proof was provided in Ref. [8].) Furthermore, in Ref. [9], Brezzi et al. employ an innovative technique to generalize it to a family of schemes with parameters which might be tuned to achieve the discrete maximum principle. But by now such a scheme still does not exist, and the computational cost of the mimetic finite difference schemes is an issue too.

There is also considerable advance in mixed finite element type methods for distorted general meshes [10], [12], but they are quite expensive and the 3-D case is still under development [25], [21]. A worthy alternative, the control volume finite element (CVFE) methods (often called finite volume element methods) [5], [11], [32], are locally mass conservative and applicable on flexible grids, but they do not handle discontinuous diffusion coefficients well [14].

In this paper, we develop a new multi-point flux approximation method, called enriched multi-point flux approximation (EMPFA), based on more consistent pressure approximations in the interaction region. The proposed method does not produce oscillatory solutions and works for general matching or non-matching polygonal and polyhedral meshes, including meshes with non-planar interfaces. Expected convergence rate is also achieved for the new method, Moreover, the enriched multi-point flux approximation is equivalent to the two-point flux approximation for the K-orthogonal grids.

The remaining of the paper is organized as follows. In Section 2, we restate the original multi-point flux approximation on quadrilateral grids. Section 3 describes our new multi-point flux approximation for two-dimensional grids, shows its numerical convergence, and compares the numerical results with the multi-point flux approximation method. In Section 4, we extend the technique in Section 3 to polyhedral meshes. Finally, we summarize the paper and make some concluding remarks in Section 5.