A locally mass-conservative enriched Petrov-Galerkin method without penalty for the Darcy flow in porous media

In this work, we consider the coupled Darcy flow and species transport within a single flowing phase in porous media. The velocity computed from the Darcy flow will be used for simulations of the species transport in the porous media to solve for the concentration.

For the Darcy flow, it has the following form:

where $\Omega$ is a bounded polygonal/polyhedral domain in $\mathbb{R}^{d}$ ($d$ = 2 or 3) with boundary $\partial\Omega$, and n denotes the unit outward normal vector on $\partial\Omega$. The unknowns are the pressure $p$ and the velocity u. Assume that the conductivity $\boldsymbol{K}$ (the ratio of the permeability and the viscosity) is uniformly symmetric positive definite and bounded from above and below. For simplicity, we assume $\boldsymbol{K}=K\boldsymbol{I}$ in the following, where $\boldsymbol{I}$ is the identity matrix, and $K$ is a piecewisely positive and bounded constant. We impose a Dirichlet boundary condition on $\Gamma_{D}$ and a Neumann boundary condition on $\Gamma_{N}$ ($\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$).

The species transport equation in porous media is described as follows:

where $\phi$ is porosity of porous media, $\widehat{f}$ is the source or sink term, and $T_{f}$ denotes the final simulation time. Let $\partial\Omega=\Gamma_{in}\cup\Gamma_{out}$ with $\Gamma_{in}=\{x\in\partial\Omega:\textit{{u}}\cdot\textit{{n}}<0\}$ as the inflow boundary and $\Gamma_{out}=\{x\in\partial\Omega:\textit{{u}}\cdot\textit{{n}}\geq 0\}$ as the outflow/no-flow boundary. $c_{I}$ denotes the inflow concentration and an initial concentration $c_{0}$ is specified to close the system.

Let $\mathcal{T}_{h}$ be a family of partitions of $\Omega$ composed of triangles if $d$ = 2 or tetrahedrons if $d$ = 3. Our method can be extended to quadrilaterals if $d$ = 2 or hexahedra if $d$ = 3, but the space of $\mathbb{P}_{k}(T)$ below will need to be replaced by $\mathbb{Q}_{k}(T)$. For convenience of presentation, the following necessary notations are given:

where $M^{k}(\mathcal{T}_{h})$ denotes the space of piecewise discontinuous polynomials of degree $k$, while $M_{0}^{k}(\mathcal{T}_{h})$ consists of continuous polynomials. $M^{b}(\mathcal{T}_{h})$ is the space of piecewise bubble functions and $\alpha_{T}$ means the degree of freedom.

We now define the average and jump for $v\in H^{s}(\mathcal{T}_{h})$. Let $T_{i},T_{j}\in\mathcal{T}_{h}$ and $e\in\partial T_{i}\cap\partial T_{j}$ with $n_{e}$ exterior to $T_{i}$. Denote

We comment that the average $\{v\}$ or the jump $[v]$ actually means the function value $v$ itself if it is on the Dirichlet boundaries. Next we denote the upwind value of concentration $c^{*}|_{e}$ as follows:

We now propose the EPG method for the Darcy flow, which enriches the approximation trial space of the CG method with bubble functions and enriches the approximation test space of the CG method with piecewise constant functions. The trial and test finite element spaces $V_{h}^{1}$ and $V_{h}^{2}$ for the EPG method are denoted as

For the Darcy flow, we introduce the bilinear form $a(p_{h},q_{h})$ and the linear functional $l(q_{h})$ as follows:

where the set of all faces ($d=3$) or edges ($d=2$) on interior boundaries, Dirichlet boundaries, Neumann boundaries for $\mathcal{T}_{h}$ is denoted by $\Gamma_{h},\Gamma_{h,D}$ and $\Gamma_{h,N}$, respectively. On each face, a unit normal vector $\textit{{n}}_{e}$ is uniquely chosen. It is obvious that the EPG method uses a weak formulation similar to the EG method, but it does not require any penalty term due to the continuity of the trial space.

The EPG method for the Darcy flow is to seek $p_{h}\in V_{h}^{1}$ satisfying

Now we present a decoupled algorithm for the EPG method, which contains two steps. In the first step, we consider the continuous situation. In the second step, we solve the bubble function, which is regarded as a post-processing correction of the first step.

We decompose $p_{h}$ and $q_{h}$ in (2.10) as follows:

The EPG method can also be rewritten as follows:

If we choose $q_{h}^{c}|_{\Gamma_{h,D}}=0$, then we have

The bubble function $p_{h}^{b}$ could be carefully chosen such that $a(p_{h}^{b},q_{h}^{c})=0$. More details on the construction of bubble functions will be provided in the subsequent sections. Now the decoupled algorithm of the EPG method (2.10) is briefly summarized as follows:

We note that the decoupled algorithm of the EPG method divides the solution $p_{h}$ into two parts. The first part of the solution is the standard solution of the CG method and the well-posedness of $p_{h}^{c}$ is well-known. Therefore, in order to verify the well-posedness of $p_{h}$, we need to focus on the well-posedness of $p_{h}^{b}$, which is in the second step of the algorithm. The Babuška theory (cf. [7, 11, 47]) reveals that the problem (2.10) or (2.12) is well-posed if and only if the continuity and coercivity (inf-sup) condition of the bilinear forms are satisfied.

Since the well-posedness of $p^{c}_{h}$ is well-known, we will provide the well-posedness of $p_{h}^{b}$ in the next section. Now we define the energy norm as follows:

The energy norm defined by (2.13) can be applied to both $V_{h}^{1}$ and $V_{h}^{2}$ spaces. Since $p_{h}\in V_{h}^{1}$ is globally continuous, the jump term vanishes on $\Gamma_{h}$. Therefore, for any $p_{h}\in V_{h}^{1}$ and $q_{h}\in V_{h}^{2}$, we can describe the energy norms as follows:

In fact, based on the arbitrariness of $p_{h}$ and $q_{h}$ in (3.1), we can take $p_{h}$ as $p_{h}^{b}$ and $q_{h}$ as $q_{h}^{0}$, which directly shows the continuity condition of the bilinear form in the second step of the decoupled algorithm as follows:

In order to show the well-posedness of (2.12), next we only need to prove the inf-sup condition of $a(p_{h}^{b},q_{h}^{0})$.

Denote $p_{h}^{b}=\sum\limits_{T\in\mathcal{T}_{h}}\alpha_{T}b_{T}$ where $p_{h}^{b}$ is a global bubble function while $b_{T}$ is its local bubble basis function and $\alpha_{T}$ is a constant on $T$. Let $\alpha$ be a piecewise constant function with $\alpha|_{T}=\alpha_{T}$ for any $T\in\mathcal{T}_{h}$. The definition of $[\alpha]$ on the face ($d=3$) or the edge ($d=2$) $e$ is based on the previous definition in (2.7). For any $p^{b}_{h}\in M^{b}(\mathcal{T}_{h})$ and $q^{0}_{h}\in M^{0}(\mathcal{T}_{h})$, we denote $\interleave p_{h}^{b}\interleave_{M^{b}(\mathcal{T}_{h})}:=\interleave p_{h}^{b}\interleave_{V^{1}_{h}}$ and $\interleave q_{h}^{0}\interleave_{M^{0}(\mathcal{T}_{h})}:=\interleave q_{h}^{0}\interleave_{V^{2}_{h}}$. Obviously, we have $\interleave p_{h}^{b}\interleave_{M^{b}(\mathcal{T}_{h})}=\|K^{1/2}\nabla p_{h}^{b}\|_{0,\mathcal{T}_{h}}$. Besides, we define another energy norm $\interleave\cdot\interleave_{*}$ for $p_{h}^{b}$ by

Next we start to construct the bubble function such that the orthogonality $a(p_{h}^{b},q_{h}^{c})=0$ is satisfied for any $p_{h}^{b}\in M^{b}(\mathcal{T}_{h})$ and $q_{h}^{0}\in M^{0}(\mathcal{T}_{h})$. For any $p_{h}^{b}=\sum\limits_{T\in\mathcal{T}_{h}}\alpha_{T}b_{T}$, we define $b_{T}$ in any $T\in\mathcal{T}_{h}$ as follows:

where $J=C_{k+d-2}^{d}$ and $b_{T,i}$ is a $d$-dimensional one-sided bubble function with the undetermined coefficient $\beta_{i}$, i.e., $b_{T,i}=\beta_{i}\frac{\lambda_{1}^{2}\cdots\lambda_{d+1}^{2}}{\lambda_{i}}$. Here, $\gamma_{j}$ is also an undetermined coefficient and $\psi_{j}$ is a basis function in the barycentric coordinate form of $M^{k-2}(T)$.

We should illustrate that $\beta_{i}$ and $\gamma_{j}$ are well-defined. Firstly, we set the undetermined coefficient $\beta_{i}$ satisfying

Thus, $\beta_{i}$ is well-defined by $\beta_{i}=\frac{1}{\int_{e_{i}}K\nabla\lambda_{i}\cdot\textit{{n}}_{i}\frac{\lambda_{1}^{2}\cdots\lambda_{d+1}^{2}}{\lambda_{i}^{2}}d\sigma},$ where the denominator can be easily observed as nonzero. Secondly, it is easy to get that $a(p_{h}^{b},q_{h}^{c})=0$ when $k=1$. For $k\geq 2$, we will show that the construction of $b_{T}$ in (3.6) satisfies the orthogonality property $a(p_{h}^{b},q_{h}^{c})=0$ if the parameter $\{\gamma_{j}\}^{J}_{j=1}$ are well chosen. If $a(p_{h}^{b},q_{h}^{c})=0$ is satisfied, for any $T\in\mathcal{T}_{h}$, the above $b_{T}$ satisfies that

where $(\cdot,\cdot)_{T}$ means the inner product on $T$. Then, we get the linear system for $\{\gamma_{j}\}^{J}_{j=1}$ as follows:

where $m_{ij}=(K\lambda_{1}^{2}\cdots\lambda_{d+1}^{2}\psi_{j},\psi_{i})_{T}$. Next, we will show that $\{\gamma_{j}\}^{J}_{j=1}$ is well defined by the above linear system. Define a new space $W_{h}(T)$ as follows:

where $\overline{\psi}_{i}:=K^{1/2}\lambda_{1}\cdots\lambda_{d+1}\psi_{i}$ is the basis function for $W_{h}(T)$. The coefficient matrix can be simplified as follows:

which yields a nonsingular coefficient matrix because $\{\overline{\psi}_{i}\}_{i=1}^{J}$ are linearly independent. Besides, the right-hand term of the above linear system is nonzero, which shows that $\{\gamma_{j}\}^{J}_{j=1}$ can be well determined. Thus, we can suitably choose the bubble function such that the orthogonality property $a(p_{h}^{b},q_{h}^{c})=0$ holds true.

Next we show the coercivity estimate for the bilinear form $a(p_{h}^{b},q_{h}^{0})$.

The continuity and the two inf-sup conditions for the bilinear form $a(p_{h}^{b},q_{h}^{0})$ have been proved through Theorem 3.1 and Theorem 3.3. Thus, we can directly see that the solution $p_{h}^{b}$ of (2.12) is well-posed in the second step of the decoupled algorithm by the Babuška theory (cf. [7, 11, 47]). Since the well-posedness of the solution $p_{h}^{c}$ in the first step to solve is well-known, we can easily conclude that the solution $p_{h}$ of (2.10) for the EPG method is well-posed. By the Babuška theory (cf. [7, 11, 47]) again, the following two inf-sup conditions for the bilinear form $a(p_{h},q_{h})$ are satisfied:

We now introduce several necessary interpolations, which will be used in the analysis of convergence. Firstly, let $I_{h}^{c}:L^{1}(\Omega)\rightarrow M_{0}^{k}(\mathcal{T}_{h})$ be the Clément interpolation operator [17]. The interpolation error estimates of $I_{h}^{c}$ have been well-known in the literature (cf. [17]) as follows: for any $v\in H^{s}(\Omega)$, $s>1$, there hold