Abstract
Finitevolume schemes, which honor pressure and fluxcontinuity conditions, is developed using double quadrature $({q}_{{u}_{1}},{q}_{{u}_{2}})$, referred as double family scheme. The scheme is applicable to solve the elliptic pressure equation used in reservoir simulation. Schemes are applicable on both regular cartesian and unstructured triangular meshes. The scheme is defined over a controlvolume distributed formulation. The scheme can be applied to both diagonal and full permeability tensor elliptic pressure equation with discontinuous coefficients. The scheme removes the first order errors, which are introduced by standard reservoir simulation schemes when applied to full tensor flow. The scheme is quantified with help of a quadrature rule. When the scheme is applied to highly heterogeneous and anisotropic porous media it does not honor maximum principle resulting in unstable solution with oscillatory behavior. The numerical solution is termed nonmonotonicity for high anisotropy ratios with results showing oscillations in the numerical pressure solution. In this paper a double $({q}_{{u}_{1}},{q}_{{u}_{2}})$ quadrature flux continuous schemes is presented, which with specific choice of quadrature $({q}_{{u}_{1}},{q}_{{u}_{2}})$ helps in improved stability of the numerical solutions. Numerical convergence of the scheme is also demonstrated with help of a number of numerical test cases and schemes impact on monotonicity behavior is also demonstrated with numerical examples.
Highlights
 Double quadrature
 Gridding
 Numerical Solution
 Numerical Convergence
1. Introduction
Underground porous media reservoirs of hydrogeological flows have complex description in terms of both their physical geometry and geology. Rapid variation in permeability with strong anisotropy is common occurrence in subsurface reservoir simulation, be it for hydrocarbon flows or for underground hydrological studies. Subsurface reservoir simulation requires discretization of continuous flux and pressure equation in order to honor correct local physical interface conditions between the simulation grid blocks, with strong discontinuities and anisotropy in subsurface properties, primarily permeabilities. The derivation of the algebraic flux continuity conditions for full permeability tensor discretization operators has lead to an efficient and robust locally conservative fluxcontinuous control volume distributed (CVD) finitevolume schemes for determining the discrete pressure and velocity fields in subsurface reservoirs [19]. These schemes are also known as multipoint flux approximation or MPFA [10]. Further schemes of this type are presented in [1113] and more recently in [14, 15]. These numerical schemes could be applied to the diagonal and full tensor elliptic pressure equation with rapid variations in permeability and remove the first order error in numerical approximation in typical reservoir simulation results, which are often based on twopoint flux approximation, when applied to full tensor flow approximation. Some other numerical schemes of similar types that also preserve continuity of fluxes have also been developed using mixed finite element methods [1618] and discontinuous galerkin finite element methods [19, 20]. Numerical convergence of these numerical scheme on unstructured and structured meshes has also been presented earlier in a number of papers [9, 2123].
When these numerical schemes are applied to strongly anisotropic heterogeneous porous media they fail to satisfy a maximum principle (similar to finite eminent, mixed finite element and other finitevolume methods) and result in loss of solution stability, particularly, for very high anisotropy ratios, which often lead to spurious oscillations in the numerical pressure solution. Some work has been carried out previously with aim of preserving stability of the numerical solutions for highly heterogeneous cases with skewed meshes, please see references [3, 6, 24].
Monotonicity conditions were first derived in [1] for an Mmatrix and later in [25] for a monotone matrix. Numerical treatment based on mesh/grid optimization methods have been used to help improve stability of the discrete system [26]. Still, all of the published numerical schemes show only a limited impact and their stability is often subject to the order of magnitude of the permeability anisotropy ratios. A more promising strategy based on Fluxsplitting techniques has been presented in [7, 8, 27], which helps in locally imposing maximum principle and resulting numerical solutions are oscillation free. Therefore, guaranteeing a more stable numerical pressure solution.
Conditions for the (one parameter) family of fluxcontinuous schemes to possess symmetric positive definite matrices and for such schemes to have Mmatrices with diagonal dominance and negative offdiagonals are presented in [1, 2]. The discretized system obtained for the schemes in the case of a full tensor are shown to be conditionally diagonally dominant and thus the Mmatrix property is conditional [1, 28]. For very large anisotropic ratios with skewed meshes, the resulting numerical system for these schemes can be unstable (as in the case of other numerical techniques like FEM & FVM) and the numerical pressure solutions can consequently exhibit oscillatory behavior. This paper is focused on seeking the most robust numerical approximation with help of double parameter $({q}_{{u}_{1}},{q}_{{u}_{2}})$ scheme. Numerical convergence of the double parameter schemes is demonstrated with help of series of numerical test cases.
This paper is organized as follows: Description of the single phase flow problem encountered in reservoir simulation with respect to the general tensor pressure equation is given in Section 2. Details of the numerical construction of the double parameter scheme with discretization in physical space are presented in Section 3. Section 4 details the precise conditions required to obtain a stable numerical solution and shows the performance of the scheme with respect to stability conditions. Section 5 presents the numerical convergence results with help of a series of test cases. Section 6 presents numerical examples that demonstrate the impact of scheme quadrature choice on monotonicity behavior with help of some numerical examples. Finally, conclusions are presented in Section 7.
2. Problem description and governing equations
The typical formulation of Darcy’s law states that the problem is to find the pressure $p$ satisfying:
over a given domain $\mathrm{\Omega}$, subjected to a given (Neumann/Dirichlet) boundary conditions on a given boundary $\partial \mathrm{\Omega}$. The right hand side term $\mathbf{M}$ represents a flow rate and $\nabla \mathrm{}=({\partial}_{x},{\partial}_{y})$. The permeability Matrix $\mathbf{K}$ can either be a diagonal or full cartesian tensor, which has the general form:
The resulting pressure equation with full permeability matrix is assumed to be elliptic such that:
The permeability tensor can be discontinuous across internal boundaries of the domain $\mathrm{\Omega}$. The boundary conditions, which are imposed here are of Dirichlet and Neumann type. For incompressible fluid flow in porous media – pressure need to be specified at atleast at one point in the domain to obtain a numerical solutions. For subsurface reservoir simulation, Neumann boundary conditions on $\partial \mathrm{\Omega}$ requires noflux on the domain boundary such that $(K\nabla p)\cdot \widehat{n}=0$, where $\widehat{n}$ is the outward normal vector to $\partial \mathrm{\Omega}$.
2.1. Elliptic pressure equation
The pressure equation presented in the section above is given in the physical space w.r.t the classical Cartesian coordinate system. In this section we present the elliptic pressure equation in a general curvilinear coordinate system that is defined with respect to a uniform dimensionless transform space with a $(\xi ,\eta )$ coordinate system. Choosing domain ${\mathrm{\Omega}}_{p}$ to represent a control volume comprised of surfaces that are tangential to constant $(\xi ,\eta )$ respectively, Eq. (1) is integrated over ${\mathrm{\Omega}}_{p}$ via the Gauss divergence theorem to yield:
where $\partial {\mathrm{\Omega}}_{p}$ is the boundary of ${\mathrm{\Omega}}_{p}$ and $\widehat{n}$ is the unit outward normal. Spatial derivatives are computed using:
where $J(x,y)={x}_{\xi}{y}_{\eta}{x}_{\eta}{y}_{\xi}$ is the Jacobian. Resolving the $x$, $y$ components of velocity along the unit normals to the curvilinear coordinates $(\xi ,\eta )$, e.g., for $\xi =$ constant, $\widehat{\mathbf{n}}ds=({y}_{\eta},{x}_{\eta})d\eta $ gives rise to the general tensor flux components:
where general tensor $\mathbf{T}$ has elements defined by:
${T}_{22}=\frac{{k}_{11}{y}_{\xi}^{2}+{k}_{22}{x}_{\xi}2{k}_{12}{x}_{\xi}{y}_{\xi}}{J},$
${T}_{12}=\frac{{k}_{12}\left({x}_{\xi}{y}_{\eta}+{x}_{\eta}{y}_{\xi}\right)\left({k}_{11}{y}_{\eta}{y}_{\xi}+{k}_{22}{x}_{\eta}{x}_{\xi}\right)}{J},$
and the closed integral can be written as:
where e.g. ${\u2206}_{\xi}F$ is the difference in net flux with respect to $\xi $ and $\stackrel{~}{F}={T}_{11}{p}_{\xi}+{T}_{12}{p}_{\eta}$, $\stackrel{~}{G}={T}_{12}{p}_{\xi}+{T}_{22}{p}_{\eta}$. This condition demonstrates that any such scheme applicable to a full matrix permeability tensor could also be applied to nonKOrthogonal grids. Note that ${T}_{11}$, ${T}_{22}\ge 0$ and ellipticity of $\mathbf{T}$ follows from Eqs. (3 and 7). Full permeability tensors can arise from upscaling, and local orientation of the grid and the permeability field. For example, by Eq. (7), a diagonal anisotropic Cartesian tensor leads to a full tensor on a curvilinear orthogonal grid [23] (a grid is called orthogonal if all grid lines intersect at a right angle).
3. Fluxcontinuous finite volume schemes
Finite volume schemes, which honor local pressure and flux continuity conditions and are locally conservative have been presented in [16]. The schemes have been developed for different grid types structured and unstructured grids in physical space and transform space. Numerical convergence for a range of quadrature in physical and transform space have been presented in [5]. The schemes have continuity of pressure and flux across primal grid cell controlvolume interfaces (Control Volume is defined as the volume over which mass can transfer in and out across its boundary (Pal 2006)). Here we present a short summary of the structured cell centred quadrilateral formulation. The nine point support of the numerical scheme is shown in Fig. 1(a)). The numerical scheme has flow and rock variables (porosity, permeability etc.) defined at the cell centre, with this definition, the discrete approximation points are at then located at the centres of the primal grid cells, which are also the controlvolumes. For each group of four nodes, four triangles are drawn as in Fig. 1(a). Each group of nodes are defined within a dualcell (each group of four cellcentered nodes surrounding a primal grid vertex defines the fundamental corners of a dual cell) which is obtained by joining cell centres with cell edge midpoints as indicated by the dashed contour in Fig. 1(b). The dual cells partition the primal cells (or controlvolumes) into subcells (the dual cells partition the primal quadrilateral grid cells (or control volumes) into subquadrilateral cells, which are called subcells). Two faces of each subcell also define sub faces of two faces of the parent controlvolume.
Fig. 1a) Fluxcontinuity at specific north, south, east and west location n,s,e,w on the controlvolume subcell faces, b) location of quadrature points (qu1,qu2) on the subcell face
a)
b)
3.1. $({\mathit{q}}_{{\mathit{u}}_{1}},{\mathit{q}}_{{\mathit{u}}_{2}})$ quadrature schemes
Flux continuous schemes are formulated using the principle of pressure and flux continuity conditions on the controlvolume subcell faces. The continuity conditions are shown conceptually with help of the Fig. 1(a). It can be seen in the figure that the continuity conditions is met at the four positions $(n,s,e,w)$, corresponding to north, south, east and west faces. The area is shown with help of shaded triangles and their point of contact in the Fig. 1(b). On each controlvolume cell subface the point of continuity is parameterized with help of the variable quadrature $q=({q}_{{u}_{1}},{q}_{{u}_{2}})$ depending on the location of continuity at the subface, please see Fig. 1(b) $(0<{q}_{{u}_{1}},{q}_{{u}_{2}}\le 1]$. For any given controlvolume subcell, the points of continuity of flux and pressure could be located anywhere in the ranges $(0<{q}_{{u}_{1}},{q}_{{u}_{2}}\le 1]$ on the two faces of each subcell. The precise values of $({q}_{{u}_{1}},{q}_{{u}_{2}})$ define the quadrature points, thereby, resulting in the formulation corresponding to families schemes, which are fluxcontinuous. Fig. 1(b) shows a conceptual plot for quadrature points located at ${q}_{{u}_{1}}=0.01$ and ${q}_{{u}_{2}}=1$. Cell face pressures ${P}_{n}$, ${P}_{e}$, ${P}_{s}$, ${P}_{w}$ are introduced at $n$, $s$, $e$, $w$ locations. Pressure subtriangles are defined with local triangular support imposed within each quarter (subcell) of the dualcell as shown in Fig. 1(b). Pressure $p$, in local cell coordinates, is piecewise linear over each triangle. The physical space fluxcontinuity conditions for cells 1 to 4, sharing a common grid vertex inside the dualcell are expressed as:
${F}_{S}=\frac{1}{2}({T}_{11}{p}_{\xi}+{T}_{12}{p}_{\eta}){}_{S}^{1}=\frac{1}{2}({T}_{11}{p}_{\xi}+{T}_{12}{p}_{\eta}){}_{S}^{2},$
${F}_{E}=\frac{1}{2}({T}_{12}{p}_{\xi}+{T}_{22}{p}_{\eta}){}_{E}^{2}=\frac{1}{2}({T}_{12}{p}_{\xi}+{T}_{22}{p}_{\eta}){}_{E}^{3},$
${F}_{W}=\frac{1}{2}\left({T}_{12}{p}_{\xi}+{T}_{22}{p}_{\eta}\right){}_{W}^{1}=\frac{1}{2}\left({T}_{12}{p}_{\xi}+{T}_{22}{p}_{\eta}\right){}_{W}^{4}.$
The parametric variation in $q=({q}_{1},{q}_{2})$ is illustrated further using the subcell example of Fig. 1(c), with subcell containing subtriangle $(1,S,W)$. Let ${r}_{1}=({x}_{1},{y}_{1})$ denote the coordinates of the cellcentre and ${r}_{S}=({x}_{S},{y}_{S})$,${r}_{W}=({x}_{W},{y}_{W})$ denote the local continuity coordinates. Then it is understood that the continuity position is a function of $q$ with ${r}_{S}({q}_{{u}_{1}},{q}_{{u}_{2}})$ and ${r}_{W}({q}_{{u}_{1}},{q}_{{u}_{2}})$.
Piecewise constant fluxes are then computed on every pressure subtriangles belonging to the subcells of the dualcell as shown in Fig. 1(b). The local linear pressure $p$, is expanded in subtriangle coordinates. The Darcy flux approximation for subtriangle $(1,S,W)$ is given below:
and:
Using Eqs. (10, 11) the discrete Darcy velocity is defined as:
where $\mathbf{K}$ is the local permeability tensor of cell 1 and dependency of $\nabla {p}_{h}$ on quadrature point arises through:
where approximate ${r}_{\xi}\left(q\right)$ and ${r}_{\eta}\left(q\right)$ are defined by Eq. (11). The normal flux at the left hand side of $S$ (Fig. 1(b)) is resolved along the outward normal vector $d{L}_{S}=(\mathrm{\Delta}{y}_{{r}_{3},S},\mathrm{\Delta}{x}_{{r}_{3},S})=\frac{1}{2}(\mathrm{\Delta}{y}_{32},\mathrm{\Delta}{x}_{32})$ (Fig. 1(b) and is defined in terms of the general tensor $T=T\left(q\right)$ as:
and the resulting coefficients of $\frac{1}{2}({p}_{\xi},{p}_{\eta}){}_{S}^{1}$ denoted by ${T}_{11}{}_{S}^{1}$ and ${T}_{12}{}_{S}^{2}$ are subcell (physicalspace) approximations of the general tensor components (Eq. (9)) at the left hand face at S, and are functions of $q$. A similar expression for flux is obtained at the right hand side of S from cell 2 (Fig. 1(b). Similarly subcell fluxes are resolved on the two sides of the other faces at W, N and E. Flux continuity is then imposed across the four cell interfaces at the four positions N, S, E and W (Fig. 1(a) which are specified according to quadrature point $\mathbf{q}$. The local physical space flux continuity conditions are now defined in the dual cell and expressed as follows.
The above system of Eq. (9) is now expressed as:
where $F=({F}_{N},{F}_{S},{F}_{E},{F}_{W}{)}^{T}$ are the fluxes defined in the dualcell and ${P}_{f}=({p}_{N},{p}_{S},{p}_{E},{p}_{W}{)}^{T}$ are the interface pressures. Similarly ${P}_{v}=({p}_{1},{p}_{2},{p}_{3},{p}_{4}{)}^{T}$ are the cell centered pressures. Thus the four interface pressures are expressed in terms of the four cell centered pressures. Using Eq. (15), ${P}_{f}$ is now expressed in terms of ${P}_{v}$ to obtain the dualcell flux and coefficient matrix:
Thus, the cellface pressures are eliminated from the flux by being determined locally in terms of the cell centered pressures in a preprocessing step, avoiding introduction of the interface pressure equations into the assembled discretization matrix. The Eq. (16) can also be written as:
where the entries of matrix $A$ are accumulated inverse tensor elements and $\mathrm{\Delta}{P}_{v}=({p}_{21},{p}_{32},{p}_{34},{p}_{41}{)}^{T}$ are the differences of vertex pressures. Consistency of the formulation follows from Eq. (17) which shows that flux is zero for constant potential. The relationship with the mixed method can also be deduced from Eq. (17), [2] from which it also follows that the above systems Eqs. (16), (17) represent the generalisation of the standard flux with harmonic coefficients to general elements with families of schemes defined by quadrature point $q$, see [2, 4, 9, 2123] for details. The physical space formulation does not posses a symmetric discretization matrix for arbitrary quadrilaterals, however, transform space (cell and subcell) formulations that are symmetric positive definite are presented in [1, 4]. Flux continuity in the case of a generaltensor is obtained while maintaining the standard single degree of freedom per cell. Since the continuity equations depend on both ${p}_{\xi}$ and ${p}_{\eta}$ (unless a diagonal tensor is assumed with cellface midpoint quadrature resulting in a 2point flux), the interface pressures ${p}_{f}=({p}_{N},{p}_{S},{p}_{E},{p}_{W}{)}^{T}$ are locally coupled and each group of four interface pressures is determined simultaneously in terms of the group of four cell centered pressures whose union contains the continuity positions. Finally for a structured grid the scheme is defined by:
where $i$, $j$ are the integer coordinates of the central quadrilateral cell, Fig. 1(a)) and:
${F}_{i,j+1/2}={F}_{{E}_{i1/2,j+1/2}}+{F}_{{W}_{i+1/2,j+1/2}},$
where $i+1/2$, $j+1/2$ denote the “integer” coordinates of the top right hand side dualcell, Fig. 1(a)). The unstructured formulation is presented in [22, 23].
4. Monotonicity
Fluxcontinuous schemes results in a discrete matrix which has five to nine diagonals in two dimensions and three to twentyseven diagonals in three dimensions on a structured mesh. The discrete matrix can be expressed as:
where $\mathbf{A}$ is the matrix of discrete operator, $p$ is the unknown pressure and $l$ is the source term. Ideally the discrete system obtained in this way should results in an Eq. (20), which should be a nonsingular matrix and should satisfy maximum principle that is comparable to the continuous counterpart of the discrete problem. The condition ensures that the numerical solution obtained in this way is free from any nonphysical oscillations. The maximum principle is satisfied for the discrete system given in 20 if the discrete matrix operator $\mathbf{A}$ is a monotone. The matrix operator, of the discrete system, $\mathbf{A}$ is a nonsingular monotone matrix if it obeys the following condition [29]:
where $\mathbf{O}$ is a null matrix. A sufficient condition for a maximum principle that can ensure that nonphysical oscillations do not occur in the resulting discrete solution is that $\mathbf{A}$ is an Mmatrix, i.e. monotone or positive definite with ${a}_{i,j}\le 0$, $i\ne j$.
4.1. Mmatrix conditions
The following set of conditions are often easier to verify.
$\mathbf{A}$ is an Mmatrix if:
${a}_{i,j}\le 0,\forall i,j,i\ne j,$
${\mathrm{\Sigma}}_{j}{a}_{i,j}\ge 0,\forall i.$
In addition, $\mathbf{A}$ must be either strictly diagonally dominant (strict inequality in the latter of Eq. (22) or weakly diagonally dominant with strict inequality for at least one row, $\mathbf{A}$ must also be irreducible. The first Mmatrix analysis for schemes of this type is presented in [1], where conditions for ninenode flux continuous schemes to be an MMatrix are:
and $\eta \left(q\right)$ is a function of quadrature point. One of the essential conditions here is that:
which is only sufficient for ellipticity [1] and therefore quite limiting on the range of tensors that are applicable. Tensors that are elliptic with:
and are such that $\left{T}_{\mathrm{1,2}}\right>min({T}_{\mathrm{1,1}},{T}_{\mathrm{2,2}})$ violate the MMatrix criteria of Eq. (23) and expose the MMatrix limit.
5. Numerical results: numerical convergence study
While numerical convergence tests for transform and physical space formulations on cartesian and triangular unstructured meshes in twoand three dimensions have previously been performed by [5, 2123] for the default member of the family of fluxcontinuous MPFA schemes i.e. for $q=$1.
The aim of the numerical convergence study presented here is to test the effect of variable quadrature points ${q}_{1}$, ${q}_{2}$ on the numerical convergence of the scheme. Although a number of quadrature combinations exist for ${q}_{1}$, ${q}_{2}$, only two specific combination of quadrature points are evaluated ${q}_{{u}_{1}}=$ 0.5, ${q}_{{u}_{2}}=$ 0.1 and ${q}_{{u}_{1}}=$ 0.1, ${q}_{{u}_{2}}=$ 0.5.
In this section numerical convergence results are presented cartesian grids in two dimensions. In all test cases the permeability field remains unchanged under grid refinement, ensuring that each problem is invariant with respect to each grid level for the convergence study. The classical ${L}_{2}$ norm is used to measure errors in fluxes and pressure, which is defined for fluxes and pressure as:
where, $f=\mathbf{v}\cdot \mathbf{n}$ (where $\mathbf{v}=\mathbf{K}\nabla p$) is the edge normal flow velocity. Subscript $h$ refers to numerical solution. Further ${V}_{i}$ is the volume of the grid cell $i$, and ${Q}_{j}$ is the volume associated with edge $j$ (where two cells are separated by edge $j$). The grid refinement levels used for the ${L}_{2}$ norm calculation where mesh sizes of 16×16, 32×32, 64×64, and 128×128 are used for all test cases in 2D. In each case Dirichlet boundary conditions are prescribed via the exact solutions.
5.1. Numerical test 1: uniformly discontinuous permeability
In this case the pressure field is piecewise quadratically varying over the domain shown in Fig. 2(a). Fig. 2(b) shows analytical pressure results on this test case. Fig. 2(c) shows typical mesh used for the numerical convergence in this case. The domain is discontinuous with discontinuity aligned along the line $x=1/2$, and the analytical solution is given as:
$K=\left\{\begin{array}{ll}\left(\begin{array}{ll}50& 0\\ 0& 1\end{array}\right),& x<1/2,\\ \left(\begin{array}{ll}1& 0\\ 0& 10\end{array}\right),& x\ge 1/2,\end{array}\right.$
$\alpha =\frac{{K}_{11}{}_{r}}{{K}_{11}{}_{l}},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\beta =\frac{{K}_{22}{}_{l}}{{K}_{22}{}_{l}},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{a}_{r}=1,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}f=\frac{4{a}_{r}}{\left(\alpha 2\right)\beta +1},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{b}_{r}=(\beta 1)f,$
${c}_{r}=f,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{d}_{r}=\frac{{c}_{r}{K}_{11}{}_{r}}{{K}_{22}{}_{r}},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{c}_{l}=\alpha \beta {c}_{r},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{d}_{l}={d}_{r}.$
The imposed top boundary flux is also discontinuous at $x=1/2$, resulting in a discontinuous tangential flux across the domain. The computed numerical solution and the plots showing ${L}_{2}$ norm of pressure for the varying quadrature ${q}_{1}$ and ${q}_{2}$ are shown in Fig. 3(a) and (b), respectively. ${O}^{{h}^{1.0179}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.5}$ and ${q}_{{u}_{2}}=\text{0.1}$ and ${O}^{{h}^{1.0496}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.1}$ and ${q}_{{u}_{2}}=\text{0.5}$.
Fig. 2Case 1: a) Discontinuity of the domain under consideration, b) exact analytical solution for pressure, c) sample mesh used for numerical convergence
a)
b)
c)
Fig. 3Case 1: a) convergence plots for qu1, qu2. b) exact numerical pressure – physical space
a)
b)
5.2. Numerical test 2: nonuniform discontinuous permeability field
In this case the medium is divided into two parts as, see Fig. 4(a). The permeability is discontinuous and permeability ratio is 1/10 across the medium discontinuity. The discontinuity is aligned along the line $rx+sy=0$, where $r=\mathrm{t}\mathrm{a}\mathrm{n}(\pi /3)/(1+\mathrm{t}\mathrm{a}\mathrm{n}(\pi /3\left)\right)$ and $s=1/(1+\mathrm{t}\mathrm{a}\mathrm{n}(\pi /3\left)\right)$.The pressure is piecewise linear and varies as:
The diagonal permeability tensor $\mathbf{K}=c\mathbf{I}$, where $c=10$ for $rx+sy<0$ and $c=1$ for $rx+sy>0$. The numerical solution shown in Fig. 4(b) was obtained using a grid aligned along the discontinuity. Typical mesh used for the numerical convergence is shown in Fig. 4(c). Numerical convergence for varying ${q}_{{u}_{1}}$ and ${q}_{{u}_{2}}$ quadrature is shown in Fig. 5(a) and the numerical pressure Physical space for the varying quadrature is shown in Fig. 5(b). ${O}^{{h}^{2.0218}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=0.5$ and ${q}_{{u}_{2}}=0.1$ and ${O}^{{h}^{1.9949}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=0.1$ and ${q}_{{u}_{2}}=0.5$.
Fig. 4Case 1: a) discontinuity of the domain under consideration, b) exact analytical solution for pressure, c) sample mesh used for numerical convergence
a)
b)
c)
Fig. 5Case 1: a) convergence plots for qu1, qu2, b) Exact numerical pressure – physical space
a)
b)
5.3. Numerical test 3: permeability field with internal discontinuity
This case involves an internal discontinuity. The problem involves a rectangular domain with internal discontinuous permeability variation as indicated in Fig. 6(a). The exact solution in each case takes the form:
The exact analytical pressure solution is shown in Fig. 6(b) and the cartesian mesh used in this case is shown in Fig. 6(c). Numerical convergence for varying ${q}_{{u}_{1}}$ and ${q}_{{u}_{2}}$ quadrature is shown in Fig. 7(a) and the numerical pressure Physical space for the varying quadrature is shown in Fig.7(b). ${O}^{{h}^{0.9955}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.5}$ and ${q}_{{u}_{2}}=\text{0.1}$ and ${O}^{{h}^{0.9859}}$ accuracy is observed in numerical pressure solution for ${q}_{{u}_{1}}=\text{0.1}$ and ${q}_{{u}_{2}}=\text{0.5}$.
Fig. 6a) Internal discontinuity along θ=π/2, b) exact analytical pressure, c) sample mesh used for numerical convergence
a)
b)
c)
Fig. 7Case 1: a) convergence plots for q1, q2, b) exact numerical pressure – physical space
a)
b)
6. Numerical results: test of monotonicity behavior
In this section a series of numerical test are presented to test the quasimonotonic property of the variable quadrature points ${q}_{{u}_{1}}\text{,}$${q}_{{u}_{2}}$ scheme. The test examples involve subsurface permeability discontinuity with highly heterogeneous and anisotropic permeability tensor.
6.1. Test case 1: v shape permeability distribution
In this test case, a two dimensional domain with a permeability tensor that changes direction in anisotropy halfway across the domain going from $\mathbf{K}$= [1, 0.99, 0.99, 1] to $\mathbf{K}$ = [1, –0.99, –0.99, 1], as shown in Fig. 8, is investigated. A point source is placed in the centre of the 2D domain (of unit length in $x$ and $y$ direction) and Dirichlet conditions are applied. In this test case due to discontinuity in permeability tensor Mmatrix conditions are tested under variable coefficient conditions. Fig. 9(a) shows the results with ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=1$ with visible spurious oscillations. In the second case using the different choice of quadrature parameter with ${q}_{{u}_{1}}={q}_{{u}_{2}}=q=$ 0.01 oscillation free results are obtained, see Fig. 9(b).
Fig. 8Change in direction of anisotropy in Permeability Tensor at y=0.5
Fig. 9a) Numerical pressure contours with instabilities like oscillations qu1=qu2=q=1, b) numerical pressure contours free of instabilities e.g., no oscillation. qu1=qu2=q=0.01
a)
b)
Results using the double family with ${q}_{{u}_{1}}=$ 0.005025125, ${q}_{{u}_{2}}=\text{1}$ are shown in Fig. 10. In these cases the Mmatrix conditions are satisfied. Using a reduced support scheme (7point – constant tensor) motivated by the above analysis yields the result in Fig. 11(a). If the alternative orientation is selected the scheme still has reduced support and the scheme does not posses an Mmatrix and the result obtained is shown in Fig. 11(b). There are clear oscillations and spread in the solution in this case, which serves as a test of the angled condition.
6.2. Test case 2: strong anisotropic permeability distribution
In this case the same boundary conditions apply as in case 2. The permeability tensor changes direction in anisotropy halfway across the domain as before, now with stronger tensor field, going from $\mathbf{K}=$ [750.25 432.58 432.58 250.75] to $\mathbf{K}=$[750.25 –432.58 –432.58 250.75] as indicated in Fig. 8. In this case the tensor has a principal anisotropy ratio of 1/1000, while the tensor is elliptic, the sufficient condition for ellipticity, is clearly violated. Consequently, the conditions for an Mmatrix are also violated.
Fig. 10Numerical pressure contours free of instabilities like oscillations for double quadrature qu1= 0.005025125, qu2= 1
Fig. 11a) Numerical pressure 7pt scheme tensor not aligned with scheme stencil, b) numerical pressure contours 7pt scheme tensor aligned with scheme stencil
a)
b)
Fig. 12a) Numerical pressure contours 9pt scheme tough tensor, q=1
Results are presented for ${q}_{{u}_{1}}=\text{1}$ and ${q}_{{u}_{2}}=\text{1}$ Fig. 12 for a 64×64 grid. There are strong oscillations in the solution showing clear violation of the maximum principle. Next, we turn to the reduced support scheme with angled approximation with orientation of opposite sign to the cross terms, results are in Fig. 13(a) which show severe oscillations. Finally, we employ the reduced support scheme with angled approximation and orientation that matches the sign of the cross terms, consistent with the Mmatrix analysis, for the same 64×64 grid Fig. 13(b). While oscillations are present, they are seen to be reduced. The results clearly show a reduction in oscillations in each case, however at each grid level and as convergence is approached, the consistent reduced support scheme yields almost oscillation free results or quasimonotonic.
Fig. 13a) Numerical pressure contours 7pt scheme tough tensor not aligned with scheme stencil, b) numerical pressure 7pt scheme tough tensor aligned with scheme stencil
a)
b)
7. Conclusions
A two parameter $({q}_{{u}_{1}},{q}_{{u}_{2}})$ scheme is presented, which is locally conservative and grantees continuity of flux and pressure. Numerical convergence of the scheme is also presented, which shows numerical accuracy of order ${O}^{h}{O}^{{h}^{2}}$ on the tested examples.
When a classical scheme is used that violates the governing conditions oscillations appear in the solutions. In this paper results are presented for highly anisotropic fulltensor fields that verify that double family scheme offers quasimonotonic behavior. When the governing conditions are satisfied the for given choice of quadrature the numerical pressure solution is free oscillations.
The variable support $({q}_{{u}_{1}},{q}_{{u}_{2}})$ schemes are able to compute numerical pressure solutions with minimal or no oscillations and has much improved resolution compared to schemes relying on fixed quadrature rules. The scheme also proves effective for high anisotropy ratios.
References

M. G. Edwards and C. F. Rogers, “Finite volume discretization with imposed flux continuity for the general tensor pressure equation,” Computational Geosciences, Vol. 2, No. 4, pp. 259–290, 1998, https://doi.org/10.1023/a:1011510505406

M. G. Edwards, “Unstructured, controlvolume distributed, fulltensor finitevolume schemes with flow based grids,” Computational Geosciences, Vol. 6, No. 3/4, pp. 433–452, 2002, https://doi.org/10.1023/a:1021243231313

M. G. Edwards, “MMatrix flux splitting for general full tensor discretization operators on structured and unstructured grids,” Journal of Computational Physics, Vol. 160, No. 1, pp. 1–28, May 2000, https://doi.org/10.1006/jcph.2000.6418

M. G. Edwards and M. Pal, “Symmetric positive definite general tensor discretization: a family of subcell flux continuous schemes on cell centred quadrilateral grids,” (in press).

M. Pal, M. G. Edwards, and A. R. Lamb, “Convergence study of a family of flux‐continuous, finite‐volume schemes for the general tensor pressure equation,” (in Press), International Journal for Numerical Methods in Fluids, Vol. 51, No. 910, pp. 1177–1203, Jul. 2006, https://doi.org/10.1002/fld.1211

Mayur Pal and Michael G. Edwards, “Effective upscaling using a family of fluxcontinuous finitevolume schemes for the pressure equation,” in ACME UK, Apr. 2006.

M. Pal and M. G. Edwards, “Family of fluxcontinuous finitevolume schemes with improved monotonicity,” in Proceedings of the 10th European Conference on the Mathematics of Oil Recovery, 2006.

M. Pal and M. G. Edwards, “FluxSplitting Schemes for improved montonicity of discrete solution of elliptic equation with highly anisotropic coefficients,” in ECCOMAS CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, 2006.

M. Pal and M. G. Edwards, “Non‐linear flux‐splitting schemes with imposed discrete maximum principle for elliptic equations with highly anisotropic coefficients,” International Journal for Numerical Methods in Fluids, Vol. 66, No. 3, pp. 299–323, May 2011, https://doi.org/10.1002/fld.2258

I. Aavatsmark, “Introduction to multipoint flux approximation for quadrilateral grids,” Computational Geosciences, Vol. 6, No. 3/4, pp. 405–432, 2002, https://doi.org/10.1023/a:1021291114475

S. H. Lee, P. Jenny, and H. A. Tchelepi, “A finitevolume method with hexahedral multiblock grids for modeling flow in porous media,” Computational Geosciences, Vol. 6, No. 3/4, pp. 353–379, 2002, https://doi.org/10.1023/a:1021287013566

S. Verma, “Flexible grids for reservoir simulation,” Ph.D. Thesis, Stanford University, 1996.

S. K. Khattri, “Analyzing finitevolume for singlephase flow in porous media,” Journal of Porous Media, Vol. 10, No. 2, pp. 109–123, 2007.

O. Iliev and I. Rybak, On Approximation Property of Multipoint Flux Approximation Method. 2007.

A. A. Lukyanov, “A meshless multipoint flux approximation method,” in Numerical Analysis and Applied Mathematics ICNAAM 2012: International Conference of Numerical Analysis and Applied Mathematics, 2012, https://doi.org/10.1063/1.4756377

L. J. Durlofsky, “A Triangle Based Mixed Finite ElementFinite Volume Technique for Modeling Two Phase Flow through Porous Media,” Journal of Computational Physics, Vol. 105, No. 2, pp. 252–266, Apr. 1993, https://doi.org/10.1006/jcph.1993.1072

A. S. Abushaikha, D. V. Voskov, and H. A. Tchelepi, “Fully implicit mixedhybrid finiteelement discretization for general purpose subsurface reservoir simulation,” Journal of Computational Physics, Vol. 346, pp. 514–538, Oct. 2017, https://doi.org/10.1016/j.jcp.2017.06.034

A. S. Abushaikha and K. M. Terekhov, “A fully implicit mimetic finite difference scheme for general purpose subsurface reservoir simulation with full tensor permeability,” Journal of Computational Physics, Vol. 406, p. 109194, Apr. 2020, https://doi.org/10.1016/j.jcp.2019.109194

B. Rivière, M. F. Wheeler, and K. Banaś, “Discontinuous Glaerkin method applied to a single phase flow in porous media,” Computational Geosciences, Vol. 4, No. 4, pp. 337–349, 2000, https://doi.org/10.1023/a:1011546411957

T. F. Russell and M. Wheeler, “Mathematics of reservoir simulation: finite element and finite difference methods for continuous flows in porous media,” in Frontiers in Applied Mathematics SIAM, 1983, p. 35.

M. Pal and M. G. Edwards, “The competing effects of discretization and upscaling – a study using the qfamily of CVDMPFA,” 11th European Conference on the Mathematics of Oil Recovery, Vol. 68, No. 1, Sep. 2008, https://doi.org/10.3997/22144609.20146372

M. Pal and M. G. Edwards, “A family of multipoint flux approximation schemes for general element types in two and three dimensions with convergence performance,” International Journal for Numerical Methods in Fluids, Vol. 69, No. 11, pp. 1797–1817, Aug. 2012, https://doi.org/10.1002/fld.2665

Mayur Pal, “Families of controlvolume distributed CVD(MPFA) finite volume schemes for the porous medium pressure equation on structured and unstructured grids,” Ph.D. thesis, Swansea University, 2007.

J. M. Nordbotten and G. T. Eigestad, “Discretization on quadrilateral grids with improved monotonicity properties,” Journal of Computational Physics, Vol. 203, No. 2, pp. 744–760, Mar. 2005, https://doi.org/10.1016/j.jcp.2004.10.002

J. Nordbotten and I. Aavatsmark, “Monotonicity conditions for control volume methods on uniform parallelogram grids in homogeneous media,” Computational Geosciences, Vol. 9, pp. 61–72, 2005.

M. J. Mlacnik and L. J. Durlofsky, “Unstructured grid optimization for improved monotonicity of discrete solutions of elliptic equations with highly anisotropic coefficients,” Journal of Computational Physics, Vol. 216, No. 1, pp. 337–361, Jul. 2006, https://doi.org/10.1016/j.jcp.2005.12.007

T. Aadland, I. Aavatsmark, and H. K. Dahle, “Performance of a flux splitting scheme when solving the singlephase pressure equation discretized by MPFA,” Computational Geosciences, Vol. 8, No. 4, pp. 325–340, Dec. 2004, https://doi.org/10.1007/s105960043774y

G. T. Eigestad, I. Aavatsmark, and M. Espedal, “Symmetry and MMatrix issues for the OMethod on an unstructured grid,” Computational Geosciences, Vol. 6, No. 3/4, pp. 381–404, 2002, https://doi.org/10.1023/a:1021239130404

O. Axelsson, Iterative Solution Methods. Cambridge University Press, 1994, https://doi.org/10.1017/cbo9780511624100

I. Aavatsmark, G. T. Eigestad, B. T. Mallison, and J. M. Nordbotten, “A compact multipoint flux approximation method with improved robustness,” Numerical Methods for Partial Differential Equations, Vol. 24, No. 5, pp. 1329–1360, Sep. 2008, https://doi.org/10.1002/num.20320