Numerical convergence of the family of flux-continuous schemes with variable quadrature (qu1,qu2) for single phase flow in porous media

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:

and the closed integral can be written as:

where e.g. ${∆}_{\xi }F$ is the difference in net flux with respect to $\xi$ and $\tilde{F}=T_{11}{p}_{\xi }+T_{12}{p}_{\eta }$, $\tilde{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 non-K-Orthogonal 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.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 control-volume 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 control-volume cell sub-face 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 sub-face, please see Fig. 1(b) $(0<q_{u_1},q_{u_2}\le 1]$. For any given control-volume 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 flux-continuous. 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 sub-triangles are defined with local triangular support imposed within each quarter (sub-cell) of the dual-cell as shown in Fig. 1(b). Pressure $p$, in local cell coordinates, is piecewise linear over each triangle. The physical space flux-continuity conditions for cells 1 to 4, sharing a common grid vertex inside the dual-cell are expressed as:

The parametric variation in $q=(q_1,q_2)$ is illustrated further using the sub-cell example of Fig. 1(c), with sub-cell containing sub-triangle $(1,S,W)$. Let $r_1=(x_1,y_1)$ denote the coordinates of the cell-centre 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 sub-triangles belonging to the sub-cells of the dual-cell as shown in Fig. 1(b). The local linear pressure $p$, is expanded in sub-triangle coordinates. The Darcy flux approximation for sub-triangle $(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 sub-cell (physical-space) 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 sub-cell 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 dual-cell 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 dual-cell flux and coefficient matrix:

Thus, the cell-face 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, 21-23] for details. The physical space formulation does not posses a symmetric discretization matrix for arbitrary quadrilaterals, however, transform space (cell and sub-cell) formulations that are symmetric positive definite are presented in [1, 4]. Flux continuity in the case of a general-tensor 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 cell-face midpoint quadrature resulting in a 2-point 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:

where $i+1/2$, $j+1/2$ denote the “integer” coordinates of the top right hand side dual-cell, Fig. 1(a)). The unstructured formulation is presented in [22, 23].

4.1. M-matrix conditions

The following set of conditions are often easier to verify.

$\mathbf{A}$ is an M-matrix if:

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 M-matrix analysis for schemes of this type is presented in [1], where conditions for nine-node flux continuous schemes to be an M-Matrix 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 M-Matrix criteria of Eq. (23) and expose the M-Matrix limit.

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:

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: non-uniform 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.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 M-matrix 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 M-matrix conditions are satisfied. Using a reduced support scheme (7-point – 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 M-matrix 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 M-matrix are also violated.

Fig. 10Numerical pressure contours free of instabilities like oscillations for double quadrature qu1= 0.005025125, qu2= 1

Fig. 11a) Numerical pressure 7-pt scheme tensor not aligned with scheme stencil, b) numerical pressure contours 7-pt scheme tensor aligned with scheme stencil

a)

b)

Fig. 12a) Numerical pressure contours 9-pt 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 M-matrix 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 quasi-monotonic.

Fig. 13a) Numerical pressure contours 7-pt scheme tough tensor not aligned with scheme stencil, b) numerical pressure 7-pt scheme tough tensor aligned with scheme stencil

a)

b)