High-Order Hybridizable Discontinuous Galerkin Formulation for One-Phase Flow Through Porous Media

Abstract

We present a stable high-order hybridizable discontinuous Galerkin (HDG) formulation coupled with high-order diagonal implicit Runge–Kuta (DIRK) schemes to simulate slightly compressible one-phase flow through porous media. The HDG stability depends on the selection of a single parameter and its definition is crucial to ensure the stability and to achieve the high-order properties of the method. Thus, we extend the work of Nguyen et al. in J Comput Phys 228, 8841–8855, 2009 to deduce an analytical expression for the stabilization parameter using the material parameters of the problem and the Engquist-Osher monotone flux scheme. The formulation is high-order accurate for the pressure, the flux and the velocity with the same convergence rate of P+1, being P the polynomial degree of the approximation. This is important because high-order methods have the potential to reduce the computational cost while obtaining more accurate solutions with less dissipation and dispersion errors than low order methods. The formulation can use unstructured meshes to capture the heterogeneous properties of the reservoir. In addition, it is conservative at the element level, which is important when solving PDE’s in conservative form. Moreover, a hybridization procedure can be applied to reduce the size of the global linear system. To keep these advantages, we use DIRK schemes to perform the time integration. DIRK schemes are high-order accurate and have a low memory footprint. We show numerical evidence of the optimal convergence rates obtained with the proposed formulation. Finally, we present several examples to illustrate the capabilities of the formulation.

Appendices

A Jacobian Terms

In this “Appendix”, we write the Jacobian terms required to solve Eq. (25) of the Newton–Raphson method. First, we deduce the partial derivatives of the numerical convective flux, ${\hat{\mathbf {F}}}_{h}\cdot \mathbf {n}$, Eq. (12), respect to the pressure, ${p_{}}_{h}$, and its trace, $\hat{p}_{h}$

$$\begin{aligned} { \dfrac{\partial {\hat{\mathbf {F}}_{h}}\cdot \mathbf {n}}{\partial {{p_{}}_j}} = \left\{ \begin{array}{ll} \left( \dfrac{\mathbf {F}({p_{}}_h)\cdot \mathbf {n}-{\hat{\mathbf {F}}_{h}}\cdot \mathbf {n}}{({p_{}}_h-\hat{p}_{h})} \right) N_{j} &{} \mathrm {if\ } \left( \mathbf {K}\mathbf {g}\right) \cdot \mathbf {n}\ge 0, \\ 0 &{} \mathrm {if\ } \left( \mathbf {K}\mathbf {g}\right) \cdot \mathbf {n}< 0, \end{array} \right. } \end{aligned}$$

(36)

$$\begin{aligned} { \dfrac{\partial {\hat{\mathbf {F}}_{h}}\cdot \mathbf {n}}{\partial {{\hat{p}_{j}}}} = \left\{ \begin{array}{ll} \left( \dfrac{-\mathbf {F}(\hat{p}_{h})\cdot \mathbf {n}+{\hat{\mathbf {F}}_{h}}\cdot \mathbf {n}}{({p_{}}_h-\hat{p}_{h})} \right) N_{j}^{f} &{} \mathrm {if\ } \left( \mathbf {K}\mathbf {g}\right) \cdot \mathbf {n}\ge 0, \\ \mathbf {F}'(\hat{p}_{h}) N_{j}^{f} &{} \mathrm {if\ } \left( \mathbf {K}\mathbf {g}\right) \cdot \mathbf {n}< 0, \end{array} \right. } \end{aligned}$$

(37)

where ${p_{}}_{h}$ and $\hat{p}_{h}$ are defined in Eqs. (16) and (17), respectively. Note that, for $p_{h} = \hat{p}_{h}$ the first equations of Eqs. (36) and (37) becomes

$$\begin{aligned} \lim _{p_{h} \rightarrow \hat{p}_{h}} \dfrac{\partial {{\hat{\mathbf {F}}}_{h}}\cdot \mathbf {n}}{\partial {{p_{}}_j}} = \lim _{p_{h} \rightarrow \hat{p}_{h}} \dfrac{\partial {{\hat{\mathbf {F}}}_{h}}\cdot \mathbf {n}}{\partial {{\hat{p}_{j}}}} = \dfrac{1}{2} \mathbf {F}'(\hat{p}_{h}) . \end{aligned}$$

The partial derivatives of Eqs. (16) to (18) are

$$\begin{aligned} \dfrac{\partial {\mathbf {q}_{}}_{h} }{\partial q_{ij}} = N_{i}{\mathbf {e}}_{j} ,\qquad \dfrac{\partial {p_{}}_{h}}{\partial {p_{}}_i} =N_{i} ,\qquad \dfrac{\partial p_{}}{\partial {{\dot{p_{}}}}_j} =N_{j}\Delta t a_{ii} ,\qquad \dfrac{\partial {{{\hat{p_{}}}}_{h}}}{\partial {p_{}}_j} =N_{j}^{f}. \end{aligned}$$

Finally, we obtain the partial derivative of ${{\hat{\mathbf {F}}}_{h}}\cdot \mathbf {n}$ respect to ${\dot{p_{}}}_j$ as

$$\begin{aligned} \dfrac{\partial {\hat{\mathbf {F}}}_{h}\cdot \mathbf {n}}{\partial {\dot{p_{}}}_j} = \dfrac{\partial {\hat{\mathbf {F}}}_{h}\cdot \mathbf {n}}{\partial {p_{}}_j} \cdot \dfrac{\partial {p_{}}_{h}}{\partial {{\dot{p_{}}}}_i} = \Delta t a_{ii}\dfrac{\partial {\hat{\mathbf {F}}}_{h}\cdot \mathbf {n}}{\partial {p_{}}_j}. \end{aligned}$$

Thus, the elemental Jacobian coefficients are

$$\begin{aligned} ({{{\mathbf {J}}}}_{{q}_{}{q}_{}}^{e})_{ij,kl}= & {} \dfrac{\partial }{\partial {{q}_{}}_{kl}}\left( {\mathbf {A}}^{-1}{\mathbf {q}_{}}_{h}, {N_{i}{\mathbf {e}}_{j}} \right) _e =({\mathbf {A}}^{-1}N_{k}{\mathbf {e}}_{l}, N_{i}{\mathbf {e}}_{j})_e \\ ({{{\mathbf {J}}}}_{{q}_{}\dot{p_{}}}^{e})_{ij,k}= & {} \dfrac{\partial }{\partial {\dot{p_{}}_k}}\left( ({{p_{}}_{i}}, \nabla \cdot N_{i}{\mathbf {e}}_{j})_e + {({{\mathbf {A}}}^{-1}{{\mathbf {q}_{}}}_{h}, N_{i}{\mathbf {e}}_{j})_e} \right) \\= & {} \Delta t a_{ii}\left( \bigg (-{\mathbf {q}_{}}_{h}\mu _{}{\mathbf {K}}^{-1}\dfrac{{\rho _{}}_{ref}c_{f}N_{k}}{{\rho _{}({p_{}}_{i})}^{2}}, N_{i}{\mathbf {e}}_{j}\bigg )_e-\bigg (N_{k}, \dfrac{\partial N_{i}}{\partial x_j}\bigg )_e \right) \\ ({{{\mathbf {J}}}}_{{q}_{}\hat{p_{}}}^{e,f})_{ij,k}= & {} \dfrac{\partial }{\partial {\hat{p_{}}}_{k}} {\langle {\hat{p_{}}}_{h}, N_{i}{\mathbf {e}}_{j}\cdot {\mathbf {n}}\rangle _{\partial e} } = \langle N_{k}^{f}, N_{i}{\mathbf {n}}_j\rangle _{\partial e} \\ ({{{\mathbf {J}}}}_{\dot{p_{}}{q}_{}}^{e})_{i,jk}= & {} \dfrac{\partial }{\partial {{q}_{}}_{jk}}\left( -({{\mathbf {q}_{}}}_{h},\nabla N_{i})_e +\langle {\mathbf {q}_{}}_{h}\cdot \mathbf {n},\ N_{i}\rangle _{\partial e} \right) = -(N_{j}{\mathbf {e}}_{k},\nabla N_{i})_e + \langle N_{j}{\mathbf {n}}_{k}, N_{i}\rangle _{\partial e} \\ ({{{\mathbf {J}}}}_{\dot{p_{}}\dot{p_{}}}^{e})_{i,j}= & {} \dfrac{\partial }{\partial {\dot{p_{}}_j}}\left( {(s{\dot{{p_{}}}}_{h}, N_{i})_e} - ({\mathbf {F}}_{h},\nabla N_{i})_e + \langle \hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}\rangle _{\partial e} + \langle \tau _{\text {diff}}\ {p_{}}_{i},\ N_{i}\rangle _{\partial e} - {\langle \tau _{\text {diff}}\ {{\hat{p_{}}}_{h},\ N_{i}\rangle _{\partial e}}} \right) \\= & {} \bigg ( \big ( c_{t}\Delta t a_{ii}(c_{r}{\phi }_{ref}\rho _{}+ \phi c_{f}{\rho _{}}_{ref})N_{j} , N_{i}\big )_e +\bigg (\phi \rho _{}({p_{}}_{i})c_{t}\dfrac{\partial {\dot{{p_{}}}}_{}}{\partial {p_{}}_{j}},N_{i}\bigg )_e \bigg ) \\&\quad -\, \bigg ( \Delta t a_{ii}\dfrac{\mathbf {K}\nabla z\mathrm {g}}{\mu _{}}2\rho _{}c_{f}{\rho _{}}_{ref}N_{j}, \nabla N_{i} \bigg )_e + {\langle \dfrac{\partial }{\partial {\dot{p_{}}}_{j}}\hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}\rangle _{\partial e}} \\&\quad +\, \bigg \langle { {\Delta t a_{ii}}\left( \dfrac{\mathbf {K}}{\mu _{}l}c_{f}{\rho _{}}_{ref}{p_{}}_{i}+\tau _{\text {diff}}\right) N_{j},N_{i}} \bigg \rangle _{\partial e} + \bigg \langle {\Delta t a_{ii}}\dfrac{\mathbf {K}}{\mu _{}l}c_{f}{\rho _{}}_{ref}\hat{p}_{h}, N_{i}\bigg \rangle _{\partial e} \\ ({{{\mathbf {J}}}}_{\dot{p_{}}\hat{p_{}}}^{e,f})_{i,j}= & {} \dfrac{\partial }{\partial {\hat{p_{}}}_{j}}\left( \langle \hat{\mathbf {F}}_{h}\cdot \mathbf {n}, N_{i} \rangle _{\partial e} - {\langle {\tau _{\text {diff}}\ {\hat{p_{}}}_{h}, N_{i} \rangle _{\partial e}}} \right) = \langle \dfrac{\partial }{\partial {\hat{p_{}}}_{j}}\hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}\rangle _{\partial e} - {\langle {\tau _{\text {diff}}N_{j}^{f}, N_{i}\rangle _{\partial e}}} \\ ({{{\mathbf {J}}}}_{\hat{p_{}}{q}_{}}^{e,f})_{i,jk}= & {} \dfrac{\partial }{\partial {{q}_{}}_{jk}} {\langle {\mathbf {q}_{}}_{h}\cdot \mathbf {n},\ N_{i}^{f} \rangle _{\partial e} } = \langle N_{j}{\mathbf {n}}_{k}, N_{i}^{f} \rangle _{\partial e} \\ ({{{\mathbf {J}}}}_{\hat{p_{}}\dot{p_{}}}^{e,f})_{i,j}= & {} \dfrac{\partial }{\partial {\dot{p_{}}_j}}\left( \langle \hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}^{f} \rangle _{\partial e} + \langle {\tau _{\text {diff}}\ {{p_{}}_{i}},\ N_{i}^{f} \rangle _{\partial e}} - \langle {\tau _{\text {diff}}\ {\hat{p}_{h}},\ N_{i}^{f} \rangle _{\partial e}} \right) \\= & {} {\langle \dfrac{\partial }{\partial {\dot{p_{}}_j}} \hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}^{f}\rangle _{\partial e}} + \bigg \langle { {\Delta t a_{ii}}\left( \dfrac{\mathbf {K}}{\mu _{}l}c_{f}{\rho _{}}_{ref}{p_{}}_{i}+\tau _{\text {diff}}\right) N_{j},N_{i}} \bigg \rangle _{\partial e} \\&\quad -\, \bigg \langle \Delta t a_{ii}{{\dfrac{\mathbf {K}}{\mu _{}l}c_{f}{\rho _{}}_{ref}N_{j}\hat{p}_{h}},\ N_{i}^{f}} \bigg \rangle _{\partial e} \\ ({{{\mathbf {J}}}}_{\hat{p_{}}\hat{p_{}}}^{e,f})_{i,j}= & {} \dfrac{\partial }{\partial {\hat{p_{}}}_{j}}\left( \langle \hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}^{f}\rangle _{\partial e} - \langle {\tau _{\text {diff}}\ {\hat{p_{}}}_{h},\ N_{i}^{f}\rangle _{\partial e}} \right) = \langle \dfrac{\partial }{\partial {\hat{p_{}}}_{j}}\hat{\mathbf {F}}_{h}\cdot \mathbf {n},\ N_{i}^{f}\rangle _{\partial e} - \langle {\tau _{\text {diff}}\ {N_{j}^{f}},\ N_{i}^{f} \rangle _{\partial e}} \end{aligned}$$

B Temporal Discretization Schemes

This second part of the “Appendix” contains the Butcher’s tables of the temporal schemes used in the examples. The stability conditions of each one are

Scheme	Stability
Backward	L-stable
DIRK2s2	A-stable
DIRK3s3	L-stable
DIRK4s6	A-stable

see [7, 23, 24, 33] for more details.

Table 6 Butcher’s table for the backward scheme

Table 7 Butcher’s table for the DIRK2-s2 scheme

Table 8 Butcher’s table for the DIRK3-s3 scheme

Table 9 Butcher’s table for the DIRK4-s6 scheme