A discontinuous control volume finite element method for multi-phase flow in heterogeneous porous media

Article history: Received 24 November 2016 Received in revised form 26 June 2017 Accepted 29 September 2017


Introduction
Many geologic porous media are highly heterogeneous: they contain domains of contrasting material properties such as porosity, permeability and capillary pressure. Moreover, these domains often have complex three dimensional (3D) geometries and spatial distribution, and the changes in material properties across the domain boundaries are effectively discontinuous, because they occur over a few tens or hundreds of grains (e.g. [ (2) and (3) lie within a high permeability domain; FEs (4), (5) and (6) lie within a low permeability domain. If the phase saturation in the high permeability domain changes, the saturation change spreads artificially into the low permeability domain via the shared CV. (B) Face-centred CV formulation implemented by [2]. Each CV is shared across at most two FEs so the spreading of saturation from the high into the low permeability domain is reduced, but not eliminated.
architectures (e.g. Fig. 1), with the capabilities of control-volume methods to produce stable, locally mass-conservative solutions ([50]). Typically, pressure and velocity, and material properties, such as permeability and porosity, are represented FE-wise, but mass conservative properties such as saturation and concentration are represented CV-wise.
In standard CVFE methods, the CVs span the boundaries between contrasting material properties ( Fig. 2(A)). This creates a problem when there are contrasts in permeability or capillary entry pressure across the boundary. A number of studies have considered the case of contrasting capillary entry pressure, where appropriate conditions need to be implemented at the boundary to ensure the correct physical representation (see [ Fig. 1(B)). Fluid saturation and component concentration are therefore smeared across neighbouring elements ([2]). This issue has led some authors to argue that, despite the advantages outlined above, CVFE methods are not suitable for modelling multiphase flow in highly heterogeneous porous media (e.g. [21], Fig. 2(A)).
Several recent papers have addressed the problem of CVs that span domain boundaries in heterogeneous porous media.
[44] separated the CVs at FE interfaces that correspond to domain boundaries, and employed the average flux between elements at the interface. However, this approach does not guarantee local mass conservation.
[7] used the boundary FEs as CVs and enforced mass conservation in a post-processing step to calculate the fluxes between CVs. However, this approach is numerically expensive and inconsistent, and was also tested on 2D elements with a radial construction of the local mass conservation equations, which is not the case in general 3D elements.
[2] constructed CVs around FE boundaries rather than FE nodes, which reduces the 'smearing' because each CV spans at most two FEs (Fig. 2(B)). However, the approach is numerically expensive, as the number of degrees of freedom at FE interfaces can be more than nine times the number of degrees of freedom at FE nodes in a tetrahedral mesh. Furthermore, the smearing effect is not eliminated, as CVs still span more than one FE. Thus, higher mesh resolution is still required at domain boundaries to reduce the 'smearing' effect.
[52] introduced a new element pair which eliminates the 'smearing' effect by employing CVs that are discontinuous between FEs. The method is numerically consistent because pressure is also discontinuous. Their results were promising, but the method was implemented only in 2D. Most models of interest in porous media flow are 3D.
The first key aim of this paper is to develop and implement a new CVFE method in 3D that does not require the use of control volumes (CVs) that span domain boundaries in heterogeneous porous media. We implement discontinuous CVs here by allowing pressure to be discontinuous between elements, using the element pair first presented by [52] which has a discontinuous 1st-order representation for pressure, and a discontinuous 2nd-order representation for velocity. The (discontinuous) pressure and velocity fields are solved using a FE approach and saturation is updated within CVs, but the CVs do