Upwinding and artificial viscosity for robust discontinuous Galerkin schemes of two-phase flow in mass conservation form

Abstract

High-order discretizations have become increasingly popular across a wide range of applications, including reservoir simulation. However, the lack of stability and robustness of these discretizations for advection-dominant problems prevent them from being widely adopted. This paper presents work towards improving the stability and robustness of the discontinuous Galerkin (DG) finite element scheme, for advection-dominant two-phase flow problems in particular. A linearized analysis of the two-phase flow equations is used to show that a standard DG discretization of the two-phase flow equations in mass conservation form results in a neutrally stable semi-discrete system in the advection-dominant limit. Furthermore, the analysis is also used to propose additional terms to the DG method which linearly stabilize the discretization. These additional terms are derived by comparing the linearized equations in mass conservation form against an upwinded pressure-saturation form of the equations. Next, a partial differential equation-based artificial viscosity method is proposed for the Buckley-Leverett and two-phase flow equations, as a means of mitigating Gibbs oscillations in high-order discretizations and ensuring convergence to physical solutions. The modified DG method with artificial viscosity is demonstrated on a two-phase flow problem with heterogeneous rock permeabilities, where the high-order discretizations significantly outperform a conventional first-order approach in terms of computational cost required to achieve a given level of error in an output of interest.

Appendices

Appendix 1: Discrete linearized analysis

Consider the following discontinuous Galerkin weak form for the incompressible two-phase flow equations in Eqs. 1 and 2, with additional stabilization terms g_w and g_n that are yet to be determined. For simplicity, boundary conditions are ignored as the goal is the development of an upwinded interior discretization. The DG method seeks a discrete solution $\mathbf {u}_{h,p} = [p_{n}, S_{w}] \in \mathcal {V}_{h,p}$ that satisfies:

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} v \rho_{w} \phi \frac{\partial S_{w}}{\partial t} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \left( \rho_{w} \lambda_{w} \mathbf{K} {\nabla} p_{w} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \lambda_{w} \mathbf{K} {\nabla} p_{w} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \lambda_{w} \mathbf{K} \left( \eta_{f} \vec{r}_{p}\left( \llbracket p_{n} \rrbracket \right) \right) \right\} \ d{\Gamma} \\ \end{array} $$

$$ \begin{array}{@{}rcl@{}} - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \lambda_{w} \mathbf{K} \left( - \eta_{f} p_{c_{S}} \vec{r}_{S}\left( \llbracket S_{w} \rrbracket \right) \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ \rho_{w} \lambda_{w} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket p_{n} \rrbracket \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ - \rho_{w} \lambda_{w} p_{c_{S}} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket S_{w} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} g_{w}(v^{\pm}, \mathbf{u}_{h,p}^{\pm}, {\nabla} \mathbf{u}_{h,p}^{\pm}; \vec{n}^{+} ) \ d{\Gamma} &= 0, \\ \forall v \in \mathcal{V}_{h,p}, \end{array} $$

(69)

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} - v \rho_{n} \phi \frac{\partial S_{w}}{\partial t} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \left( \rho_{n} \lambda_{n} \mathbf{K} {\nabla} p_{n} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{n} \lambda_{n} \mathbf{K} \left( {\nabla} p_{n} + \eta_{f} \vec{r}_{p}\left( \llbracket p_{n} \rrbracket\right) \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ \rho_{n} \lambda_{n} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket p_{n} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} g_{n}(v^{\pm}, \mathbf{u}_{h,p}^{\pm}, {\nabla} \mathbf{u}_{h,p}^{\pm}; \vec{n}^{+} ) \ d{\Gamma} &= 0, \\ \forall v \in \mathcal{V}_{h,p} . \end{array} $$

(70)

All of the spatial flux terms in 69 and 70, except for the g_w and g_n terms, are obtained by expanding the BR2 operator given in Eq. 14 for the wetting and non-wetting equations separately. The lifting operators for the primary variables p_n and S_w are represented by $\vec {r}_{p}$ and $\vec {r}_{S}$ respectively.

The DG weak form given by 69 and 70 is then linearized about some base pressure distribution $\bar {p}_{n}(\vec {x},t)$ and a constant saturation value $\bar {S}_{w}$, yielding the linearized weak form of the incompressible two-phase flow equations,

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} v \rho_{w} \phi \frac{\partial S_{w}^{\prime}}{\partial t} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \left( \rho_{w} \bar{\lambda}_{w_{S}} \mathbf{K} {\nabla} \bar{p}_{w} S_{w}^{\prime} + \rho_{w} \bar{\lambda}_{w} \mathbf{K} {\nabla} p_{w}^{\prime} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \bar{\lambda}_{w_{S}} \mathbf{K} {\nabla} \bar{p}_{w} S_{w}^{\prime} + \rho_{w} \bar{\lambda}_{w} \mathbf{K} {\nabla} p_{w}^{\prime} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \bar{\lambda}_{w} \mathbf{K} \left( \eta_{f} \vec{r}_{p}^{\prime}\left( \llbracket p_{n}^{\prime} \rrbracket \right) \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \bar{\lambda}_{w} \mathbf{K} \left( - \eta_{f} \bar{p}_{c_{S}} \vec{r}_{S}^{\prime}\left( \llbracket S_{w}^{\prime} \rrbracket \right) \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ \rho_{w} \bar{\lambda}_{w} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket p_{n}^{\prime} \rrbracket \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ - \rho_{w} \bar{\lambda}_{w} \bar{p}_{c_{S}} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket S_{w}^{\prime} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} g_{w}^{\prime} \ d{\Gamma} = 0 &, \\ \forall v \in \mathcal{V}_{h,p}, \end{array} $$

(71)

and,

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} - v \rho_{n} \phi \frac{\partial S_{w}^{\prime}}{\partial t} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \left( \rho_{n} \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} S_{w}^{\prime} + \rho_{n} \bar{\lambda}_{n} \mathbf{K} {\nabla} p_{n}^{\prime} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{n} \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} S_{w}^{\prime} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{n} \bar{\lambda}_{n} \mathbf{K} \left( {\nabla} p_{n}^{\prime} + \eta_{f} \vec{r}_{p}^{\prime}\left( \llbracket p_{n}^{\prime} \rrbracket\right) \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ \rho_{n} \bar{\lambda}_{n} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket p_{n}^{\prime} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} g_{n}^{\prime} \ d{\Gamma} &= 0, \\ \forall v \in \mathcal{V}_{h,p}, \end{array} $$

(72)

where it is assumed that $\llbracket \bar {p}_{n} \rrbracket = \llbracket \bar {S}_{w} \rrbracket = 0$, and therefore $\vec {\bar {r}}_{p}(\llbracket \bar {p}_{n} \rrbracket ) = \vec {\bar {r}}_{S}(\llbracket \bar {S}_{w} \rrbracket ) = 0$.

Taking the weighted sum of the linearized weak form equations, ρ_n ×(Eq. 71) + ρ_w × (Eq. 72), produces the discrete weak form of the “pressure” equation:

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \rho_{w} \rho_{n} \left( (\bar{\lambda}_{w_{S}} \mathbf{K} {\nabla} \bar{p}_{w} + \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} ) S_{w}^{\prime} \right) \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \rho_{w} \rho_{n} \left( \bar{\lambda}_{w} \mathbf{K} {\nabla} p_{w}^{\prime} + \bar{\lambda}_{n} \mathbf{K} {\nabla} p_{n}^{\prime} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} (\bar{\lambda}_{w_{S}} \mathbf{K} {\nabla} \bar{p}_{w} + \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} ) S_{w}^{\prime} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} (\bar{\lambda}_{w} \mathbf{K} {\nabla} p_{w}^{\prime} + \bar{\lambda}_{n} \mathbf{K} {\nabla} p_{n}^{\prime} ) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} (\bar{\lambda}_{w} + \bar{\lambda}_{n}) \mathbf{K} \left( \eta_{f} \vec{r}_{p}^{\prime} \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} \bar{\lambda}_{w} \mathbf{K} \left( - \eta_{f} \bar{p}_{c_{S}} \vec{r}_{S}^{\prime} \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ \rho_{w} \rho_{n} (\bar{\lambda}_{w} + \bar{\lambda}_{n}) \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket p_{n}^{\prime} \rrbracket \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ - \rho_{w} \rho_{n} \bar{\lambda}_{w} \bar{p}_{c_{S}} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket S_{w}^{\prime} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \rho_{n} g_{w}^{\prime} + \rho_{w} g_{n}^{\prime} \ d{\Gamma} = 0 &, \\ \forall v \in \mathcal{V}_{h,p}. \end{array} $$

(73)

Similarly, the weighted difference of the linearized weak form equations, $\rho _{n} \bar {\lambda }_{n} \times \text {(Eq.~71)} - \rho _{w} \bar {\lambda }_{w} \times \text {(Eq.~72)}$, produces the discrete weak form of the “saturation” equation,

$$ \begin{array}{@{}rcl@{}} \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} v \rho_{w} \rho_{n} \phi (\bar{\lambda}_{w} + \bar{\lambda}_{n}) \frac{\partial S_{w}^{\prime}}{\partial t} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \rho_{w} \rho_{n} \left( \bar{\lambda}_{w_{S}} \bar{\lambda}_{n} \mathbf{K} {\nabla} \bar{p}_{w} \right) S_{w}^{\prime} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \rho_{w} \rho_{n} \left( - \bar{\lambda}_{w} \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} \right) S_{w}^{\prime} \ d\varOmega \\ + \sum\limits_{\kappa \in \mathcal{T}_{h}} {\int}_{\kappa} {\nabla} v \cdot \rho_{w} \rho_{n} \left( - \bar{\lambda}_{w} \bar{\lambda}_{n} \bar{p}_{c_{S}} \mathbf{K} {\nabla} S_{w}^{\prime} \right) \ d\varOmega \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} \left( \bar{\lambda}_{w_{S}} \bar{\lambda}_{n} \mathbf{K} {\nabla} \bar{p}_{w} \right) S_{w}^{\prime} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ \rho_{w} \rho_{n} \left( - \bar{\lambda}_{w} \bar{\lambda}_{n_{S}} \mathbf{K} {\nabla} \bar{p}_{n} \right) S_{w}^{\prime} \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \llbracket{v} \rrbracket \cdot \left\{ - \rho_{w} \rho_{n} \bar{\lambda}_{w} \bar{\lambda}_{n} \bar{p}_{c_{S}} \mathbf{K} \left( {\nabla} S_{w}^{\prime} + \eta_{f} \vec{r}_{S}^{\prime} \right) \right\} \ d{\Gamma} \\ - \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \left\{ - \rho_{w} \rho_{n} \bar{\lambda}_{w} \bar{\lambda}_{n} \bar{p}_{c_{S}} \mathbf{K}^{T} {\nabla} v \right\} \cdot \llbracket S_{w}^{\prime} \rrbracket \ d{\Gamma} \\ + \sum\limits_{f \in {\Gamma}_{I}} {\int}_{f} \rho_{n} \bar{\lambda}_{n} g_{w}^{\prime} - \rho_{w} \bar{\lambda}_{w} g_{n}^{\prime} \ d{\Gamma} = 0 &, \\ \forall v \in \mathcal{V}_{h,p}. \end{array} $$

(74)

Appendix 2: Stability of linear ordinary differential equations

Consider the following system of linear ordinary differential equations (ODEs):

$$ \begin{array}{@{}rcl@{}} \mathbf{M} \dot{\mathbf{u}} + \mathbf{A} \mathbf{u} = \mathbf{b}, \end{array} $$

(75)

where $\mathbf {u}(t) \in \mathbb {R}^{m}$ is the solution vector, $\mathbf {M} \in \mathbb {R}^{m \times m}$ and $\mathbf {A} \in \mathbb {R}^{m \times m}$ are constant matrices, and $\mathbf {b}(t) \in \mathbb {R}^{m}$ is a forcing vector. Substituting a perturbed solution $\tilde {\mathbf {u}}(t) = \mathbf {u}(t) + \mathbf {u}^{\prime }(t)$ into the equation above and simplifying shows that the perturbations need to satisfy the following homogeneous equation:

$$ \begin{array}{@{}rcl@{}} \mathbf{M} \dot{\mathbf{u}}^{\prime} + \mathbf{A} \mathbf{u}^{\prime} = \mathbf{0}. \end{array} $$

(76)

Solving this system of linear ODEs shows that the evolution of the perturbations, $\mathbf {u}^{\prime }(t)$, is given by:

$$ \begin{array}{@{}rcl@{}} \mathbf{u}^{\prime}(t) = \sum\limits_{i} \hat{\mathbf{u}}^{\prime}_{i} e^{\omega_{i} t}, \end{array} $$

(77)

where ω_i and $\hat {\mathbf {u}}^{\prime }_{i}$ are the generalized eigenvalues and eigenvectors, respectively, of the following generalized eigenvalue problem:

$$ \begin{array}{@{}rcl@{}} \mathbf{A} \hat{\mathbf{u}}^{\prime} = - \omega \mathbf{M} \hat{\mathbf{u}}^{\prime}. \end{array} $$

(78)

It is easy to see from Eq. 77 that all generalized eigenvalues ω_i need to have negative real components for the solution perturbations to decay over time. In other words, the system of linear ODEs is stable only if Re(ω_i) < 0, for i = 1 to m.