A conforming discontinuous Galerkin finite element method for Brinkman equations

Haoning Dang Qilong Zhai [email protected] Zhongshu Zhao School of Mathematics, Jilin University, Changchun, China. School of Mathematical Sciences and Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai, China.

Abstract

In this paper, we present a conforming discontinuous Galerkin (CDG) finite element method for Brinkman equations. The velocity stabilizer is removed by employing the higher degree polynomials to compute the weak gradient. The theoretical analysis shows that the CDG method is actually stable and accurate for the Brinkman equations. Optimal order error estimates are established in $H^{1}$ and $L^{2}$ norm. Finally, numerical experiments verify the stability and accuracy of the CDG numerical scheme.

keywords:

Brinkman equations, discontinuous Galerkin, discrete weak gradient operators, polyhedral meshes.

^†^†journal: journal of computational and applied mathematics

1 Introduction

The Brinkman equations frequently appears in modeling the incompressible flow of a viscous fluid in complex porous medium. These equations extend Darcy’s law to describe the dissipation of kinetic energy caused by viscous forces, similarly to the Navier-Stokes equations [5]. Brinkman equations are also applied in many other fields, such as environmental science, geophysics, petroleum engineering, biotechnology, and so on [3, 9, 10, 14, 18].

Mathematically speaking, the Brinkman equations combine the Darcy equations and the Stokes equations by a highly varied parameter. For simplicity, we consider the following Brinkman equations in a bounded polygonal domain $\Omega\in\mathbb{R}^{2}$ : find the unknown fluid velocity ${\mathbf{u}}$ and pressure $p$ satisfying

$\displaystyle-\mu\Delta{\mathbf{u}}+\mu\kappa^{-1}{\mathbf{u}}+\nabla p$	$\displaystyle={\mathbf{f}}\quad\text{in}~{}\Omega,$	(1.1)
$\displaystyle\nabla\cdot{\mathbf{u}}$	$\displaystyle=0\quad\text{in}~{}\Omega,$	(1.2)
$\displaystyle{\mathbf{u}}$	$\displaystyle=\boldsymbol{0}\quad\text{on}~{}\partial\Omega,$	(1.3)

where $\kappa$ is permeability tensor, and $\mu$ is the fluid viscosity coefficient. For the convenience of analysis, we assume the permeability tensor $\kappa$ is piecewise constant. Since $\mu$ can be used to scale the solution ${\mathbf{u}}$ , we can take $\mu=1$ for simplicity. ${\mathbf{f}}$ represents the momentum source term. In addition, assume that there exist two positive constants $\lambda_{1}$ and $\lambda_{2}$ such that

\displaystyle\lambda_{1}\xi^{t}\xi\leq\xi^{t}\kappa^{-1}\xi\leq\lambda_{2}\xi^{t}\xi,\quad\forall\xi\in\mathbb{R}^{d}.

(1.4)

The main challenge of designing the numerical algorithm comes from the different regularity requirements of the velocity in such two extreme cases: the variational form of Brinkman equations in the Stokes limit requires $H^{1}$ -regularity, while the Darcy limit requires $H(\text{div})$ -regularity. To overcome this difficulty, many scholars have made extensive research. A natural attempt is to modify the existing Stokes or Darcy elements. We can easily find corresponding work related to Stokes based elements [4, 8, 25] and Darcy based elements [7, 12]. Another strategy is to design new formulations for Brinkman equations, such as the dual-mixed formulation [11, 13], the pseudostress-velocity formulation [17] and the vorticity-velocity-pressure formulation [2]. In addition, some new numerical methods are introduced to Brinkman equations, such as weak Galerkin methods [16, 29], virtual element methods [6, 19, 30] and much more.

The purpose of this paper is to introduce a new conforming discontinuous Galerkin (CDG) method for Brinkman equations. The CDG method based on the weak Galerkin (WG) method proposed in [20], is first proposed by Ye and Zhang in 2020 [26]. It retains the key idea of WG method, which uses the weak differential operators to approximate the classical differential operators in the variational form. In [26], the authors prove that no stabilizer is required for Poisson problem when the local Raviart-Thomas (RT) element is used to approximate the classic gradient operator. However, we know that the RT elements are only applicable to triangular and rectangular meshes. In subsequent work [27], they found that the stabilizer term can be removed from the numerical scheme constructed on polygon meshes by raising the degree of polynomial approximating the discrete weak gradient operator. The CDG method has been applied to Stokes equations [28], elliptic interface problem [23] and linear elasticity interface problem [24].

In this paper, we consider the same variational form based on gradient-gradient operators as [29]: find the unknown functions ${\mathbf{u}}\in[H^{1}_{0}(\Omega)]^{d}$ and $p\in L^{2}_{0}(\Omega)$ satisfying

	$\displaystyle(\nabla{\mathbf{u}},\nabla{\mathbf{v}})+(\kappa^{-1}{\mathbf{u}},{\mathbf{v}})+(\nabla p,{\mathbf{v}})$	$\displaystyle=({\mathbf{f}},{\mathbf{v}}),$	$\displaystyle\forall{\mathbf{v}}\in[H^{1}_{0}(\Omega)]^{d},$		(1.5)
	$\displaystyle(\nabla q,{\mathbf{u}})$	$\displaystyle=0,$	$\displaystyle\forall q\in L^{2}_{0}(\Omega).$		(1.6)

We construct the conforming discontinuous Galerkin scheme to discrete this variational form on polygon meshes. The stable term for the velocity ${\mathbf{u}}$ is removed by a new definition of weak gradient operator, thus our numerical formulation is simpler compared with the standard WG method [29]. Furthermore, we prove the well-posedness of the CDG scheme and derive the optimal error estimates for velocity and pressure, which implies the optimal convergence order for both the Stokes and Darcy dominated problems. Some numerical experiments are provided to verify our theoretical analysis.

The rest of the paper is organized as follows. In Section 2, we define two discrete weak gradient operators and construct the conforming discontinuous Galerkin scheme for Brinkman equations. Then, the well-posedness is proved in Section 3. The error equations for the CDG scheme are established in Section 4. And we prove optimal error estimates for both velocity and pressure in $H^{1}$ and $L^{2}$ norms in Section 5. In Section 6, we present some numerical experiments to verify the stability and accuracy of the CDG scheme.

2 Discrete Weak Gradient Operators

In this section, we define two discrete weak gradient operators that we’re going to use later.

Let ${\mathcal{T}}_{h}$ be a polygonal or polyhedral partition of the domain $\Omega$ and ${\mathcal{E}}_{h}$ be the set of all edges or faces in ${\mathcal{T}}_{h}$ . Assume that all cells in ${\mathcal{T}}_{h}$ are closed and simply connected, and satisfy some specific shape regular conditions in [21]. Denote the set of all interior edges or faces by ${\mathcal{E}}_{h}^{0}={\mathcal{E}}_{h}\backslash\partial\Omega$ . For each $T\in{\mathcal{T}}_{h}$ , $e\subset\partial T$ , let $h_{T}$ and $h_{e}$ be the diameter of $T$ and $e$ , respectively. And we define the size of ${\mathcal{T}}_{h}$ as $h=\max\limits_{T\in{\mathcal{T}}_{h}}h_{T}$ .

For a given integer $k\geq 1$ , the space of polynomial with degree no more than $k$ on a cell $T$ denotes by $P_{k}(T)$ . We define the space for the vector-valued functions as

\displaystyle V_{h}=\left\{{\mathbf{v}}\in[L^{2}(\Omega)]^{d}:{\mathbf{v}}|_{T}\in[P_{k}(T)]^{d},\forall T\in{\mathcal{T}}_{h}\right\}.

Denote by $V_{h}^{0}$ the subspace of $V_{h}$ that

\displaystyle V_{h}^{0}=\left\{{\mathbf{v}}:{\mathbf{v}}\in V_{h},{\mathbf{v}}=\boldsymbol{0}~{}\text{on}~{}\partial\Omega\right\}.

For the scalar-valued functions, we define

\displaystyle W_{h}=\left\{q:q\in L^{2}_{0}(\Omega),q|_{T}\in P_{k-1}(T)\right\}.

Let $T_{1}$ and $T_{2}$ be two cells in ${\mathcal{T}}_{h}$ sharing $e\in{\mathcal{E}}_{h}^{0}$ , ${\mathbf{n}}_{1}$ and ${\mathbf{n}}_{2}$ be the unit outward normal vectors of $T_{1}$ and $T_{2}$ on $e$ . In particular, when $e\subset\partial\Omega$ , we denote the unit outward normal vector of $T$ on $e$ by ${\mathbf{n}}_{e}$ . For a vector-valued function ${\mathbf{v}}\in V_{h}+[H^{1}(\Omega)]^{d}$ , we define the average $\left\{\cdot\right\}$ and jump $[\cdot]$ as follows

\left\{{\mathbf{v}}\right\}=\left\{\begin{array}[]{ll}\frac{1}{2}({\mathbf{v}}|_{T_{1}}+{\mathbf{v}}|_{T_{2}})&~{}e\in{\mathcal{E}}_{h}^{0},\\ {\mathbf{v}}|_{e}&~{}e\subset\partial\Omega,\end{array}\right.\quad[{\mathbf{v}}]=\left\{\begin{array}[]{ll}{\mathbf{v}}|_{T_{1}}\cdot{\mathbf{n}}_{1}+{\mathbf{v}}|_{T_{2}}\cdot{\mathbf{n}}_{2}&~{}e\in{\mathcal{E}}_{h}^{0},\\ {\mathbf{v}}|_{e}\cdot{\mathbf{n}}_{e}&~{}e\subset\partial\Omega.\end{array}\right.

For a scalar-valued function $q\in W_{h}$ , the average $\left\{\cdot\right\}$ and the jump $[\![\cdot]\!]$ are

\left\{q\right\}=\left\{\begin{array}[]{ll}\frac{1}{2}(q|_{T_{1}}+q|_{T_{2}})&~{}e\in{\mathcal{E}}_{h}^{0},\\ q|_{e}&~{}e\subset\partial\Omega,\end{array}\right.\quad[\![q]\!]=\left\{\begin{array}[]{ll}q|_{T_{1}}{\mathbf{n}}_{1}+q|_{T_{2}}{\mathbf{n}}_{2}&~{}e\in{\mathcal{E}}_{h}^{0},\\ q|_{e}{\mathbf{n}}_{e}&~{}e\subset\partial\Omega.\end{array}\right.

In addition, for ${\mathbf{v}}\in V_{h}^{0}+[H^{1}_{0}(\Omega)]^{d}$ , the following equations hold true:

\displaystyle\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\|_{e}=\|[{\mathbf{v}}]\|_{e},\quad\text{if}~{}e\subset\partial\Omega,\quad\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\|_{e}=\frac{1}{2}\|[{\mathbf{v}}]\|_{e},\quad\text{if}~{}e\in\mathcal{E}_{h}^{0}.

(2.1)

Similarly, for $q\in W_{h}$ , we have

\displaystyle\|q-\left\{q\right\}\|_{e}=0,\quad\text{if}~{}e\subset\partial\Omega,\quad\|q-\left\{q\right\}\|_{e}=\frac{1}{2}\|[\![q]\!]\|_{e},\quad\text{if}~{}e\in\mathcal{E}_{h}^{0}.

(2.2)

Then we give the definition of the discrete weak gradient operators. For a vector-valued function ${\mathbf{v}}\in V_{h}+[H^{1}(\Omega)]^{d}$ , the discrete weak gradient $\nabla_{w}{\mathbf{v}}$ on each cell $T$ is a unique polynomial function in $[P_{j}(T)]^{d\times d}(j>k)$ satisfying

\displaystyle(\nabla_{w}{\mathbf{v}},\tau)_{T}=-({\mathbf{v}},\nabla\cdot\tau)_{T}+\left<\left\{{\mathbf{v}}\right\},\tau{\mathbf{n}}\right>_{\partial T},\quad\forall\tau\in[P_{j}(T)]^{d\times d}.

(2.3)

Note that $j$ depends on $k$ and $n$ is the number of edges of polygon cell $T$ . For a polygonal mesh, $j=n+k-1$ [27], and in particular, $j=k+1$ when the domain is partitioned into triangles [1].

Similarly, for a scalar-valued function $q\in W_{h}$ , the discrete weak gradient $\tilde{\nabla}_{w}q$ on each cell $T$ is a unique polynomial function in $[P_{k}(T)]^{d}$ satisfying

\displaystyle(\tilde{\nabla}_{w}q,\boldsymbol{\phi})_{T}=-(q,\nabla\cdot\boldsymbol{\phi})_{T}+\left<\left\{q\right\},\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T},\quad\forall\boldsymbol{\phi}\in[P_{k}(T)]^{d}.

(2.4)

For simplicity of notations, we introduce three bilinear forms as follows:

$\displaystyle a({\mathbf{v}},{\mathbf{w}})$	$\displaystyle=(\nabla_{w}{\mathbf{v}},\nabla_{w}{\mathbf{w}})+(\kappa^{-1}{\mathbf{v}},{\mathbf{w}}),$	(2.5)
$\displaystyle b({\mathbf{v}},q)$	$\displaystyle=({\mathbf{v}},\tilde{\nabla}_{w}q),$	(2.6)
$\displaystyle s(p,q)$	$\displaystyle=\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![p]\!],[\![q]\!]\right>_{e}.$	(2.7)

Now we have the following conforming discontinuous Galerkin finite element scheme for Brinkman equations (1.1)-(1.3).

Conforming Discontinuous Galerkin Algorithm 2.1.

Find ${\mathbf{u}}_{h}\in V_{h}^{0}$ and $p_{h}\in W_{h}$ such that

	$\displaystyle a({\mathbf{u}}_{h},{\mathbf{v}})+b({\mathbf{v}},p_{h})$	$\displaystyle=({\mathbf{f}},{\mathbf{v}}),$	$\displaystyle\forall{\mathbf{v}}\in V_{h}^{0},$		(2.8)
	$\displaystyle b({\mathbf{u}}_{h},q)-s(p_{h},q)$	$\displaystyle=0,$	$\displaystyle\forall q\in W_{h}.$		(2.9)

3 The Well-Posedness of CDG Scheme

In this section, we show that the CDG numerical scheme (2.8)-(2.9) has a unique solution.

In order to analyse the well-posedness of the CDG scheme, we first define the following tri-bar norm. For a vector-valued function ${\mathbf{v}}\in V_{h}+[H^{1}(\Omega)]^{d}$ ,

\displaystyle|\!|\!|{\mathbf{v}}|\!|\!|^{2}=a({\mathbf{v}},{\mathbf{v}})=\|\nabla_{w}{\mathbf{v}}\|^{2}+\|\kappa^{-\frac{1}{2}}{\mathbf{v}}\|^{2}.

(3.1)

It is obvious that $|\!|\!|\cdot|\!|\!|$ indeed provides a norm in $V_{h}$ .

For a scalar-valued function $q\in W_{h}$ , we use the following norms in the rest of this paper,

	$\displaystyle\|\!\|\!\|q\|\!\|\!\|^{2}_{1}$	$\displaystyle=\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}+\sum_{e\in{\mathcal{E}}_{h}^{0}}h^{-1}\\|[\![q]\!]\\|_{e}^{2},$		(3.2)
	$\displaystyle\\|q\\|_{h}^{2}$	$\displaystyle=\sum_{e\in{\mathcal{E}}_{h}^{0}}h\\|[\![q]\!]\\|_{e}^{2}.$		(3.3)

By the definition of $|\!|\!|\cdot|\!|\!|$ and the Cauchy-Schwarz inequality, the following the boundedness for the bilinear form $a(\cdot,\cdot)$ holds.

Lemma 3.1.

For any ${\mathbf{v}},{\mathbf{w}}\in V_{h}+[H^{1}(\Omega)]^{d}$ , we have

\displaystyle|a({\mathbf{v}},{\mathbf{w}})|\leq|\!|\!|{\mathbf{v}}|\!|\!||\!|\!|{\mathbf{w}}|\!|\!|.

Next, we present the inf-sup condition of $b(\cdot,\cdot)$ .

Lemma 3.2.

For any $q\in W_{h}$ , there exist constants $C_{1}$ and $C_{2}$ independent of $h$ such that

\displaystyle\underset{{\mathbf{v}}\in V_{h}}{sup}\frac{|b({\mathbf{v}},q)|}{|\!|\!|{\mathbf{v}}|\!|\!|}\geq C_{1}h|\!|\!|q|\!|\!|_{1}-C_{2}\|q\|_{h}.

(3.4)

Proof.

For any ${\mathbf{v}}\in V_{h}+[H^{1}(\Omega)]^{d}$ , from the definition of $\nabla_{w}$ , we get

\displaystyle\|\nabla_{w}{\mathbf{v}}\|^{2}=\sum_{T\in\mathcal{T}_{h}}((\nabla{\mathbf{v}},\nabla{\mathbf{v}})_{T}-\left<{\mathbf{v}}-\left\{{\mathbf{v}}\right\},(\nabla_{w}{\mathbf{v}}+\nabla{\mathbf{v}})\cdot{\mathbf{n}}\right>_{\partial T}),

Using the definition of $\left\{\cdot\right\}$ , the following derivation holds true

	$\displaystyle\sum_{T\in\mathcal{T}_{h}}\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}^{2}$
	$\displaystyle\leq\sum_{e\in{\mathcal{E}}_{h}}\left(\left\\|{\mathbf{v}}\|_{T_{1}}-\frac{{\mathbf{v}}\|_{T_{1}}+{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{\partial T_{1}\cap e}^{2}+\left\\|{\mathbf{v}}\|_{T_{2}}-\frac{{\mathbf{v}}\|_{T_{1}}+{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{\partial T_{2}\cap e}^{2}\right)$
	$\displaystyle\leq C\sum_{e\in{\mathcal{E}}_{h}}\left\\|\frac{{\mathbf{v}}\|_{T_{1}}-{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{e}^{2}$
	$\displaystyle\leq C\sum_{e\in{\mathcal{E}}_{h}}\left(\\|{\mathbf{v}}\\|_{\partial T_{1}\cap e}^{2}+\\|{\mathbf{v}}\\|_{\partial T_{2}\cap e}^{2}\right)$
	$\displaystyle\leq C\sum_{T\in\mathcal{T}_{h}}\\|{\mathbf{v}}\\|_{\partial T}^{2}.$

Using the above inequality, the trace inequality (A.4) and the inverse inequality (A.5), we obtain

	$\displaystyle\\|\nabla_{w}{\mathbf{v}}\\|^{2}$	$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}\\|\nabla_{w}{\mathbf{v}}+\nabla{\mathbf{v}}\\|_{\partial T})$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+h^{-\frac{1}{2}}\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}(\\|\nabla_{w}{\mathbf{v}}\\|_{T}+\\|\nabla{\mathbf{v}}\\|_{T}))$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+Ch^{-\frac{1}{2}}\\|{\mathbf{v}}\\|_{\partial T}(\\|\nabla_{w}{\mathbf{v}}\\|_{T}+\\|\nabla{\mathbf{v}}\\|_{T}))$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+\frac{1}{2}\\|\nabla_{w}{\mathbf{v}}\\|_{T}^{2}+\frac{1}{2}\\|\nabla{\mathbf{v}}\\|_{T}^{2}+Ch^{-1}\\|{\mathbf{v}}\\|_{\partial T}^{2})$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\left(\frac{1}{2}\\|\nabla_{w}{\mathbf{v}}\\|_{T}^{2}+Ch^{-2}\\|{\mathbf{v}}\\|_{T}^{2}\right).$

Taking ${\mathbf{v}}=\kappa\tilde{\nabla}_{w}q\in V_{h}$ and using (1.4), we get the following estimate

	$\displaystyle\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|^{2}$	$\displaystyle=\\|\nabla_{w}{\mathbf{v}}\\|^{2}+\\|\kappa^{-\frac{1}{2}}{\mathbf{v}}\\|^{2}$
		$\displaystyle\leq Ch^{-2}\\|{\mathbf{v}}\\|^{2}+\\|\kappa^{-\frac{1}{2}}{\mathbf{v}}\\|^{2}$
		$\displaystyle\leq Ch^{-2}\\|\kappa\tilde{\nabla}_{w}q\\|^{2}+\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}$
		$\displaystyle\leq Ch^{-2}\lambda_{1}^{-1}\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}+\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}$
		$\displaystyle\leq Ch^{-2}\|\!\|\!\|q\|\!\|\!\|_{1}^{2}.$

Therefore, we have

	$\displaystyle\frac{b({\mathbf{v}},q)}{\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|}$	$\displaystyle=\frac{\|\!\|\!\|q\|\!\|\!\|_{1}^{2}-h^{-2}\\|q\\|_{h}^{2}}{Ch^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}}$
		$\displaystyle\geq\frac{\|\!\|\!\|q\|\!\|\!\|_{1}^{2}-h^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}\\|q\\|_{h}}{Ch^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}}$
		$\displaystyle\geq C_{1}h\|\!\|\!\|q\|\!\|\!\|_{1}-C_{2}\\|q\\|_{h}.$

∎

Lemma 3.3.

The conforming discontinuous Galerkin finite element scheme (2.8)-(2.9) has a unique solution.

Proof.

Consider the corresponding homogeneous equation ${\mathbf{f}}=\boldsymbol{0}$ , let ${\mathbf{v}}={\mathbf{u}}_{h}$ in (2.8) and $q=p_{h}$ in (2.9). Then subtracting (2.9) form (2.8), we have

\displaystyle|\!|\!|{\mathbf{u}}_{h}|\!|\!|^{2}+\|p_{h}\|_{h}^{2}=a({\mathbf{u}}_{h},{\mathbf{u}}_{h})+s(p_{h},p_{h})=0,

which implies ${\mathbf{u}}_{h}=\boldsymbol{0}$ and $\|p_{h}\|_{h}=0$ .

Let ${\mathbf{u}}_{h}=\boldsymbol{0}$ and ${\mathbf{f}}=\boldsymbol{0}$ in (2.8), we have $b({\mathbf{v}},p_{h})=0$ for any ${\mathbf{v}}\in V_{h}^{0}$ . According to the definition of $b(\cdot,\cdot)$ and $\|p_{h}\|_{h}=0$ , let ${\mathbf{v}}=\nabla p_{h}$ on each cell $T$ , it holds that

	$\displaystyle 0=b({\mathbf{v}},p_{h})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}({\mathbf{v}},\tilde{\nabla}_{w}p_{h})_{T}$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-(p_{h},\nabla\cdot{\mathbf{v}})_{T}+\left<\{p_{h}\},{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-(p_{h},\nabla\cdot{\mathbf{v}})_{T}+\left<p_{h},{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T}-\left<p_{h}-\{p_{h}\},{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}({\mathbf{v}},\nabla p_{h})_{T}-\sum_{e\in{\mathcal{E}}_{h}^{0}}\left<[\![p_{h}]\!],\{{\mathbf{v}}\}\right>_{\partial T}$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}({\mathbf{v}},\nabla p_{h})_{T}=\sum_{T\in\mathcal{T}_{h}}(\nabla p_{h},\nabla p_{h})_{T}.$

Then, we have $\nabla p_{h}=0$ on each cell $T$ . Since $p_{h}$ is continuous and $p_{h}\in L^{2}_{0}(\Omega)$ , we have $p_{h}=0$ and complete the proof of the lemma.

∎

4 Error Equations

In this section, we establish the error equations between the numerical solution and the exact solution.

Let $Q_{h}$ , ${\mathbf{Q}}_{h}$ and ${\mathbb{Q}}_{h}$ be the standard $L^{2}$ projection operators onto $[P_{k}(T)]^{d}$ , $[P_{j}(T)]^{d\times d}$ and $P_{k-1}(T)$ , respectively.

First, we recall some properties about the projection operators.

Lemma 4.1.

For the projection operators $Q_{h}$ , ${\mathbf{Q}}_{h}$ and ${\mathbb{Q}}_{h}$ , the following properties hold

$\displaystyle\nabla_{w}{\mathbf{v}}$	$\displaystyle={\mathbf{Q}}_{h}(\nabla{\mathbf{v}}),\quad\forall{\mathbf{v}}\in[H^{1}(\Omega)]^{d},$	(4.1)
$\displaystyle(\tilde{\nabla}_{w}({\mathbb{Q}}_{h}q),\boldsymbol{\phi})_{T}$	$\displaystyle=(Q_{h}(\nabla q),\boldsymbol{\phi})_{T}-\left<q-{\mathbb{Q}}_{h}q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T},$	(4.2)
	$\displaystyle\qquad\qquad\forall q\in H^{1}(\Omega),~{}\forall\boldsymbol{\phi}\in[P_{k}(T)]^{d}.$

Proof.

For any $\tau\in[P_{j}(T)]^{d\times d}$ , we have

	$\displaystyle(\nabla_{w}{\mathbf{v}},\tau)_{T}$	$\displaystyle=-({\mathbf{v}},\nabla\cdot\tau)_{T}+\left<\left\{{\mathbf{v}}\right\},\tau\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=-({\mathbf{v}},\nabla\cdot\tau)_{T}+\left<{\mathbf{v}},\tau\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=(\nabla{\mathbf{v}},\tau)_{T}$
		$\displaystyle=({\mathbf{Q}}_{h}(\nabla{\mathbf{v}}),\tau)_{T}.$

Similarly, for any $\boldsymbol{\phi}\in[P_{k}(T)]^{d}$ , we have

	$\displaystyle(\tilde{\nabla}_{w}({\mathbb{Q}}_{h}q),\boldsymbol{\phi})_{T}$	$\displaystyle=-({\mathbb{Q}}_{h}q,\nabla\cdot\boldsymbol{\phi})_{T}+\left<\left\{{\mathbb{Q}}_{h}q\right\},\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=-(q,\nabla\cdot\boldsymbol{\phi})_{T}+\left<{\mathbb{Q}}_{h}q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=-(q,\nabla\cdot\boldsymbol{\phi})_{T}+\left<q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T}-\left<q-{\mathbb{Q}}_{h}q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=(\nabla q,\boldsymbol{\phi})_{T}-\left<q-{\mathbb{Q}}_{h}q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=(Q_{h}(\nabla q),\boldsymbol{\phi})_{T}-\left<q-{\mathbb{Q}}_{h}q,\boldsymbol{\phi}\cdot{\mathbf{n}}\right>_{\partial T},$

which completes the proof of the lemma. ∎

Let ${\mathbf{e}}_{h}=Q_{h}{\mathbf{u}}-{\mathbf{u}}_{h}$ and $\varepsilon_{h}={\mathbb{Q}}_{h}p-p_{h}$ be the error functions, where $({\mathbf{u}};p)$ be the solution of (1.1)-(1.3) and $({\mathbf{u}}_{h};p_{h})$ be the solution of (2.8)-(2.9). We shall derive the error equations that ${\mathbf{e}}_{h}$ and $\varepsilon_{h}$ satisfy.

Lemma 4.2.

For any ${\mathbf{v}}\in V_{h}^{0}$ and $q\in W_{h}$ , the following equations hold true

	$\displaystyle a({\mathbf{e}}_{h},{\mathbf{v}})+b({\mathbf{v}},\varepsilon_{h})$	$\displaystyle=-l_{1}({\mathbf{u}},{\mathbf{v}})+l_{2}({\mathbf{u}},{\mathbf{v}})-l_{3}(p,{\mathbf{v}}),$		(4.3)
	$\displaystyle b({\mathbf{e}}_{h},q)-s(\varepsilon_{h},q)$	$\displaystyle=l_{4}({\mathbf{u}},q)-s({\mathbb{Q}}_{h}p,q),$		(4.4)

where

	$\displaystyle l_{1}({\mathbf{u}},{\mathbf{v}})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}{\mathbf{v}})_{T},$
	$\displaystyle l_{2}({\mathbf{u}},{\mathbf{v}})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\left<(\nabla{\mathbf{u}}-{\mathbf{Q}}_{h}\nabla{\mathbf{u}})\cdot{\mathbf{n}},{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\right>_{\partial T},$
	$\displaystyle l_{3}(p,{\mathbf{v}})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\left<p-{\mathbb{Q}}_{h}p,{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T},$
	$\displaystyle l_{4}({\mathbf{u}},q)$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\left<({\mathbf{u}}-Q_{h}{\mathbf{u}})\cdot{\mathbf{n}},q-\left\{q\right\}\right>_{\partial T},$
	$\displaystyle s({\mathbb{Q}}_{h}p,q)$	$\displaystyle=\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}p]\!],[\![q]\!]\right>.$

Proof.

Testing (1.1) by ${\mathbf{v}}\in V_{h}^{0}$ yields

\displaystyle-(\Delta{\mathbf{u}},{\mathbf{v}})+(\kappa^{-1}{\mathbf{u}},{\mathbf{v}})+(\nabla p,{\mathbf{v}})=({\mathbf{f}},{\mathbf{v}}).

Applying the definition of projection operators $Q_{h}$ , ${\mathbf{Q}}_{h}$ and ${\mathbb{Q}}_{h}$ , the definition of the weak gradient $\nabla_{w}$ and $\tilde{\nabla}_{w}$ and (4.1), we get

	$\displaystyle-(\Delta{\mathbf{u}},{\mathbf{v}})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}((\nabla{\mathbf{u}},\nabla{\mathbf{v}})_{T}-\left<\nabla{\mathbf{u}}\cdot{\mathbf{n}},{\mathbf{v}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(({\mathbf{Q}}_{h}\nabla{\mathbf{u}},\nabla{\mathbf{v}})_{T}-\left<\nabla{\mathbf{u}}\cdot{\mathbf{n}},{\mathbf{v}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-({\mathbf{v}},\nabla\cdot({\mathbf{Q}}_{h}\nabla{\mathbf{u}}))_{T}+\left<{\mathbf{v}},({\mathbf{Q}}_{h}\nabla{\mathbf{u}}-\nabla{\mathbf{u}})\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}((\nabla_{w}{\mathbf{v}},{\mathbf{Q}}_{h}\nabla{\mathbf{u}})_{T}+\left<{\mathbf{v}}-\left\{{\mathbf{v}}\right\},({\mathbf{Q}}_{h}\nabla{\mathbf{u}}-\nabla{\mathbf{u}})\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}((\nabla_{w}{\mathbf{u}},\nabla_{w}{\mathbf{v}})_{T}+\left<{\mathbf{v}}-\left\{{\mathbf{v}}\right\},({\mathbf{Q}}_{h}\nabla{\mathbf{u}}-\nabla{\mathbf{u}})\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}(Q_{h}{\mathbf{u}}),\nabla_{w}{\mathbf{v}})_{T}+l_{1}({\mathbf{u}},{\mathbf{v}})-l_{2}({\mathbf{u}},{\mathbf{v}}).$

According to the definition of $Q_{h}$ and (4.2), we have

	$\displaystyle(\nabla p,{\mathbf{v}})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(\nabla p,{\mathbf{v}})_{T}$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(Q_{h}(\nabla p),{\mathbf{v}})_{T}$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(\tilde{\nabla}_{w}({\mathbb{Q}}_{h}p),{\mathbf{v}})_{T}+\sum_{T\in\mathcal{T}_{h}}\left<p-{\mathbb{Q}}_{h}p,{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T}$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(\tilde{\nabla}_{w}({\mathbb{Q}}_{h}p),{\mathbf{v}})_{T}+l_{3}(p,{\mathbf{v}}),$

and

\displaystyle(\kappa^{-1}{\mathbf{u}},{\mathbf{v}})

\displaystyle=({\mathbf{u}},\kappa^{-1}{\mathbf{v}})=(Q_{h}{\mathbf{u}},\kappa^{-1}{\mathbf{v}})=(\kappa^{-1}Q_{h}{\mathbf{u}},{\mathbf{v}}).

Using the definition of $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ and the above equations, we have

\displaystyle a(Q_{h}{\mathbf{u}},{\mathbf{v}})+b({\mathbf{v}},{\mathbb{Q}}_{h}p)+l_{1}({\mathbf{u}},{\mathbf{v}})-l_{2}({\mathbf{u}},{\mathbf{v}})+l_{3}(p,{\mathbf{v}})=({\mathbf{f}},{\mathbf{v}}).

(4.5)

Subtracting (2.8) from (4.5), we arrive at

\displaystyle a({\mathbf{e}}_{h},{\mathbf{v}})+b({\mathbf{v}},\varepsilon_{h})+l_{1}({\mathbf{u}},{\mathbf{v}})-l_{2}({\mathbf{u}},{\mathbf{v}})+l_{3}(p,{\mathbf{v}})=0,

which completes the proof of (4.3).

Similarly, testing (1.2) by $q\in W_{h}$ yields

	$\displaystyle 0=(\nabla\cdot{\mathbf{u}},q)$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-(\nabla q,{\mathbf{u}})_{T}+\left<q,{\mathbf{u}}\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-(\nabla q,Q_{h}{\mathbf{u}})_{T}+\left<q,{\mathbf{u}}\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}((q,\nabla\cdot(Q_{h}{\mathbf{u}}))_{T}+\left<q,({\mathbf{u}}-Q_{h}{\mathbf{u}})\cdot{\mathbf{n}}\right>_{\partial T})$
		$\displaystyle=\sum_{T\in\mathcal{T}_{h}}(-(Q_{h}{\mathbf{u}},\tilde{\nabla}_{w}q)_{T},\left<q-\left\{q\right\},({\mathbf{u}}-Q_{h}{\mathbf{u}})\cdot{\mathbf{n}}\right>_{\partial T}).$

Using (2.9), we have

\displaystyle\sum_{T\in\mathcal{T}_{h}}(Q_{h}{\mathbf{u}}-{\mathbf{u}}_{h},\tilde{\nabla}_{w}q)_{T}-\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![p_{h}]\!],[\![q]\!]\right>_{e}=l_{4}({\mathbf{u}},q).

Adding $-s({\mathbb{Q}}_{h}p,q)$ to both sides of this equation, it follows that

\displaystyle\sum_{T\in\mathcal{T}_{h}}(Q_{h}{\mathbf{u}}-{\mathbf{u}}_{h},\tilde{\nabla}_{w}q)_{T}-\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}p-p_{h}]\!],[\![q]\!]\right>_{e}=l_{4}({\mathbf{u}},q)-s({\mathbb{Q}}_{h}p,q),

which implies (4.4) and we complete the proof of this lemma. ∎

5 Error Estimates

In this section, we aim to deriving the error estimates to the numerical scheme (2.8)-(2.9).

Theorem 5.1.

Let $({\mathbf{u}};p)\in[H^{1}_{0}(\Omega)\cap H^{k+1}(\Omega)]^{d}\times(L^{2}_{0}(\Omega)\cap H^{k}(\Omega))$ be the exact solution of (1.1)-(1.3) and $({\mathbf{u}}_{h};p_{h})\in V_{h}^{0}\times W_{h}$ be the solution of (2.8)-(2.9), respectively. Then there exists a constant $C$ such that

	$\displaystyle\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|+\\|\varepsilon_{h}\\|_{h}$	$\displaystyle\leq Ch^{k}(\\|{\mathbf{u}}\\|_{k+1}+\\|p\\|_{k}),$		(5.1)
	$\displaystyle\|\!\|\!\|\varepsilon_{h}\|\!\|\!\|_{1}$	$\displaystyle\leq Ch^{k-1}(\\|{\mathbf{u}}\\|_{k+1}+\\|p\\|_{k}).$		(5.2)

Proof.

Let ${\mathbf{v}}={\mathbf{e}}_{h}$ in (4.3) and $q=\varepsilon_{h}$ in (4.4), we have

\displaystyle|\!|\!|{\mathbf{e}}_{h}|\!|\!|^{2}+\|\varepsilon_{h}\|_{h}^{2}=-l_{1}({\mathbf{u}},{\mathbf{e}}_{h})+l_{2}({\mathbf{u}},{\mathbf{e}}_{h})-l_{3}(p,{\mathbf{e}}_{h})-l_{4}({\mathbf{u}},\varepsilon_{h})+s({\mathbb{Q}}_{h}p,\varepsilon_{h}).

From the estimates (A.7)-(A.11), we obtain

\displaystyle|\!|\!|{\mathbf{e}}_{h}|\!|\!|^{2}+\|\varepsilon_{h}\|_{h}^{2}\leq Ch^{k}(\|{\mathbf{u}}\|_{k+1}+\|p\|_{k})(|\!|\!|{\mathbf{e}}_{h}|\!|\!|+\|\varepsilon_{h}\|_{h}),

which implies (5.1).

Let $p=\varepsilon_{h}$ in (4.3), from the inf-sup condition (3.4), Lemma (3.1) and (A.7)-(A.9), it follows that

	$\displaystyle(C_{1}h\|\!\|\!\|\varepsilon_{h}\|\!\|\!\|_{1}-C_{2}\\|\varepsilon_{h}\\|_{h})\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|$	$\displaystyle\leq\|b({\mathbf{v}},\varepsilon_{h})\|$
		$\displaystyle\leq a({\mathbf{e}}_{h},{\mathbf{v}})+l_{1}({\mathbf{u}},{\mathbf{v}})+l_{2}({\mathbf{u}},{\mathbf{v}})+l_{3}(p,{\mathbf{v}})$
		$\displaystyle\leq\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|+Ch^{k}(\\|{\mathbf{u}}\\|_{k+1}+\\|p\\|_{k})\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|.$

Combining the estimates above, we arrive at

	$\displaystyle h\|\!\|\!\|\varepsilon_{h}\|\!\|\!\|_{1}$	$\displaystyle\leq\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|+\\|\varepsilon_{h}\\|_{h}+Ch^{k}(\\|{\mathbf{u}}\\|_{k+1}+\\|p\\|_{k})$
		$\displaystyle\leq Ch^{k}(\\|{\mathbf{u}}\\|_{k+1}+\\|p\\|_{k}),$

which completes the proof. ∎

In order to obtain $L^{2}$ error estimate, we consider the following dual problem: seek $(\boldsymbol{\psi};\xi)\in[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$ satisfying

$\displaystyle-\Delta\boldsymbol{\psi}+\kappa^{-1}\boldsymbol{\psi}+\nabla\xi$	$\displaystyle={\mathbf{e}}_{h}\quad\text{in}~{}\Omega,$	(5.3)
$\displaystyle\nabla\cdot\boldsymbol{\psi}$	$\displaystyle=0\quad\text{in}~{}\Omega,$	(5.4)
$\displaystyle\boldsymbol{\psi}$	$\displaystyle=\boldsymbol{0}\quad\text{on}~{}\partial\Omega.$	(5.5)

Assume that the following regularity condition holds,

\displaystyle\|\boldsymbol{\psi}\|_{2}+\|\xi\|_{1}\leq C\|{\mathbf{e}}_{h}\|.

(5.6)

Theorem 5.2.

\displaystyle\|{\mathbf{e}}_{h}\|\leq Ch^{k+1}(\|{\mathbf{u}}\|_{k+1}+\|p\|_{k}).

(5.7)

Proof.

Testing (5.3) by ${\mathbf{e}}_{h}$ yields

\displaystyle\|{\mathbf{e}}_{h}\|^{2}=({\mathbf{e}}_{h},{\mathbf{e}}_{h})=-(\Delta\boldsymbol{\psi},{\mathbf{e}}_{h})+(\kappa^{-1}\boldsymbol{\psi},{\mathbf{e}}_{h})+(\nabla\xi,{\mathbf{e}}_{h}).

Similar to the proof of Lemma 4.2, we have

\displaystyle-(\Delta\boldsymbol{\psi},{\mathbf{e}}_{h})

\displaystyle=(\nabla_{w}(Q_{h}\boldsymbol{\psi}),\nabla_{w}{\mathbf{e}}_{h})+l_{1}(\boldsymbol{\psi},{\mathbf{e}}_{h})-l_{2}(\boldsymbol{\psi},{\mathbf{e}}_{h}),

and

\displaystyle(\nabla\xi,{\mathbf{e}}_{h})

\displaystyle=(\tilde{\nabla}_{w}({\mathbb{Q}}_{h}\xi),{\mathbf{e}}_{h})+l_{3}(\xi,{\mathbf{e}}_{h}).

Combining the definition of $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ gives

\displaystyle\|{\mathbf{e}}_{h}\|^{2}=a(Q_{h}\boldsymbol{\psi},{\mathbf{e}}_{h})+b({\mathbf{e}}_{h},{\mathbb{Q}}_{h}\boldsymbol{\psi})+l_{1}(\boldsymbol{\psi},{\mathbf{e}}_{h})-l_{2}(\boldsymbol{\psi},{\mathbf{e}}_{h})+l_{3}(\xi,{\mathbf{e}}_{h}).

(5.8)

Similarly, testing (5.4) by $\varepsilon_{h}$ yields

\displaystyle b(Q_{h}\boldsymbol{\psi},\varepsilon_{h})=l_{4}(\boldsymbol{\psi},\varepsilon_{h}).

(5.9)

Let ${\mathbf{v}}=Q_{h}\boldsymbol{\psi}$ and $q={\mathbb{Q}}_{h}\xi$ in (4.3) and (4.4), we have

	$\displaystyle a({\mathbf{e}}_{h},Q_{h}\boldsymbol{\psi})+b(Q_{h}\boldsymbol{\psi},\varepsilon_{h})$	$\displaystyle=-l_{1}({\mathbf{u}},Q_{h}\boldsymbol{\psi})+l_{2}({\mathbf{u}},Q_{h}\boldsymbol{\psi})-l_{3}(p,Q_{h}\boldsymbol{\psi}),$		(5.10)
	$\displaystyle b({\mathbf{e}}_{h},{\mathbb{Q}}_{h}\xi)-s(\varepsilon_{h},{\mathbb{Q}}_{h}\xi)$	$\displaystyle=l_{4}({\mathbf{u}},{\mathbb{Q}}_{h}\xi)-s({\mathbb{Q}}_{h}p,{\mathbb{Q}}_{h}\xi).$		(5.11)

With (5.8)-(5.11), we obtain

	$\displaystyle\\|{\mathbf{e}}_{h}\\|^{2}=$	$\displaystyle-l_{1}(\boldsymbol{\psi},{\mathbf{e}}_{h})+l_{2}(\boldsymbol{\psi},{\mathbf{e}}_{h})-l_{3}(\xi,{\mathbf{e}}_{h})+l_{4}(\boldsymbol{\psi},\varepsilon_{h})-s({\mathbb{Q}}_{h}\boldsymbol{\psi},\varepsilon_{h})$
		$\displaystyle+(l_{1}({\mathbf{u}},Q_{h}\boldsymbol{\psi})-l_{2}({\mathbf{u}},Q_{h}\boldsymbol{\psi})+l_{3}(p,Q_{h}\boldsymbol{\psi})-l_{4}({\mathbf{u}},{\mathbb{Q}}_{h}\xi)+s({\mathbb{Q}}_{h}p,{\mathbb{Q}}_{h}\xi)).$

Let $\hat{{\mathbf{Q}}}_{h}$ be the projection operator from $[L^{2}(T)]^{d\times d}$ onto $[P_{1}(T)]^{d\times d}$ . For any $q\in P_{1}(T)$ , we have

\displaystyle(\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi},q)_{T}=(\nabla\boldsymbol{\psi},q)_{T}=-(\boldsymbol{\psi},\nabla\cdot q)_{T}+\left<\boldsymbol{\psi},q{\mathbf{n}}\right>_{\partial T}=(\nabla_{w}\boldsymbol{\psi},q)_{T}=(\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi},q)_{T},

which implies $\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi}$ is equal to $\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi}$ on each cell $T$ .

According to the definition of $\nabla_{w}$ , the following equation holds true

	$\displaystyle(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi})_{T}=$	$\displaystyle-({\mathbf{u}}-Q_{h}{\mathbf{u}},\nabla\cdot\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi})_{T}$
		$\displaystyle+\langle\left\{{\mathbf{u}}-Q_{h}{\mathbf{u}}\right\},\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi}\cdot{\mathbf{n}}\rangle_{\partial T}.$

Notice that $k$ is an integer not less than one, from the definition of the projection operator $Q_{h}$ , we obtain

\displaystyle(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi})_{T}=(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\hat{{\mathbf{Q}}}_{h}\nabla_{w}\boldsymbol{\psi})_{T}=0.

(5.12)

Then, using the projection inequalities (A.1)-(A.2), and summing over all the cells, we arrive at

		$\displaystyle\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}\boldsymbol{\psi})_{T}$
	$\displaystyle=$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}\boldsymbol{\psi}-\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi})_{T}$
	$\displaystyle=$	$\displaystyle\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla\boldsymbol{\psi}-\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi})_{T}$
	$\displaystyle\leq$	$\displaystyle\left(\sum_{T\in\mathcal{T}_{h}}\\|\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}})\\|^{2}_{T}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h}}\\|\nabla\boldsymbol{\psi}-\hat{{\mathbf{Q}}}_{h}\nabla\boldsymbol{\psi}\\|^{2}_{T}\right)^{\frac{1}{2}}$
	$\displaystyle\leq$	$\displaystyle Ch^{k+1}\\|{\mathbf{u}}\\|_{k+1}\\|\boldsymbol{\psi}\\|_{2}.$

For $l_{1}({\mathbf{u}},Q_{h}\boldsymbol{\psi})$ , we get

		$\displaystyle\|l_{1}({\mathbf{u}},Q_{h}\boldsymbol{\psi})\|$
	$\displaystyle=$	$\displaystyle\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}Q_{h}\boldsymbol{\psi})_{T}\Bigg{\|}$
	$\displaystyle\leq$	$\displaystyle\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}\boldsymbol{\psi})_{T}\Bigg{\|}+\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{u}}-Q_{h}{\mathbf{u}}),\nabla_{w}(Q_{h}\boldsymbol{\psi}-\boldsymbol{\psi}))_{T}\Bigg{\|}$
	$\displaystyle\leq$	$\displaystyle Ch^{k+1}\\|{\mathbf{u}}\\|_{k+1}\\|\boldsymbol{\psi}\\|_{2}.$

Similarly, we have

	$\displaystyle\|l_{2}({\mathbf{u}},Q_{h}\boldsymbol{\psi})\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<(\nabla{\mathbf{u}}-{\mathbf{Q}}_{h}\nabla{\mathbf{u}})\cdot{\mathbf{n}},Q_{h}\boldsymbol{\psi}-\left\{Q_{h}\boldsymbol{\psi}\right\}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<(\nabla{\mathbf{u}}-{\mathbf{Q}}_{h}\nabla{\mathbf{u}})\cdot{\mathbf{n}},Q_{h}\boldsymbol{\psi}-\boldsymbol{\psi}+\left\{\boldsymbol{\psi}-Q_{h}\boldsymbol{\psi}\right\}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|\nabla{\mathbf{u}}-{\mathbf{Q}}_{h}\nabla{\mathbf{u}}\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h}}h^{-1}\\|Q_{h}\boldsymbol{\psi}-\boldsymbol{\psi}\\|_{\partial T}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k+1}\\|{\mathbf{u}}\\|_{k+1}\\|\boldsymbol{\psi}\\|_{2}.$

According the projection inequality (A.3), we have

	$\displaystyle\|l_{3}(p,Q_{h}\boldsymbol{\psi})\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<p-{\mathbb{Q}}_{h}p,Q_{h}\boldsymbol{\psi}\cdot{\mathbf{n}}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<p-{\mathbb{Q}}_{h}p,(Q_{h}\boldsymbol{\psi}-\boldsymbol{\psi})\cdot{\mathbf{n}}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|p-{\mathbb{Q}}_{h}p\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h}}h^{-1}\\|Q_{h}\boldsymbol{\psi}-\boldsymbol{\psi}\\|_{e}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k+1}\\|p\\|_{k}\\|\boldsymbol{\psi}\\|_{2}.$

In the above derivation, we have used the following fact

\displaystyle\sum_{T\in\mathcal{T}_{h}}\left<p-{\mathbb{Q}}_{h}p,\boldsymbol{\psi}\cdot{\mathbf{n}}\right>_{\partial T}=0.

Similar to $l_{2}({\mathbf{u}},Q_{h}\boldsymbol{\psi})$ , we get

\displaystyle|l_{4}({\mathbf{u}},{\mathbb{Q}}_{h}\xi)|\leq Ch^{k+1}\|{\mathbf{u}}\|_{k+1}\|\xi\|_{1}.

Using the definition of $s(\cdot,\cdot)$ and the projection inequality (A.3), we have

	$\displaystyle\|s({\mathbb{Q}}_{h}p,{\mathbb{Q}}_{h}\xi)\|$	$\displaystyle=\Bigg{\|}\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}p]\!],[\![{\mathbb{Q}}_{h}\xi]\!]\right>_{e}\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}p-p]\!],[\![{\mathbb{Q}}_{h}\xi-\xi]\!]\right>_{e}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|p-{\mathbb{Q}}_{h}p\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h}}h\\|\xi-{\mathbb{Q}}_{h}\xi\\|_{\partial T}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k+1}\\|p\\|_{k}\\|\xi\\|_{1}.$

According to (A.7)-(A.11) and the above five estimates, it follows that

\displaystyle\|{\mathbf{e}}_{h}\|^{2}\leq(Ch(|\!|\!|{\mathbf{e}}_{h}|\!|\!|+\|\varepsilon_{h}\|_{h})+Ch^{k+1}(\|{\mathbf{u}}\|_{k+1}+\|p\|_{k}))(\|\boldsymbol{\psi}\|_{2}+\|\xi\|_{1}).

With the regularity assumption (5.6) and $H^{1}$ error estimates (5.1), we have

\displaystyle\|{\mathbf{e}}_{h}\|^{2}\leq Ch^{k+1}(\|{\mathbf{u}}\|_{k+1}+\|p\|_{k})\|{\mathbf{e}}_{h}\|,

which implies (5.7). ∎

6 Numerical Experiments

In this section, we present several examples in two dimensional domains to verify the stability and order of convergence established in Section 5.

As before, let $({\mathbf{u}};p)$ be the exact solution of (1.1)-(1.3), $({\mathbf{u}}_{h};p_{h})$ be the solution of (2.8)-(2.9). Denote ${\mathbf{e}}_{h}=Q_{h}{\mathbf{u}}-{\mathbf{u}}_{h}$ and $\varepsilon_{h}={\mathbb{Q}}_{h}p-p_{h}$ .

6.1 Example 1

Taking $\Omega=(0,1)\times(0,1)$ , the exact solution is given as follows:

\displaystyle{\mathbf{u}}=\left(\begin{matrix}sin(2\pi x)cos(2\pi y)\\ -cos(2\pi x)sin(2\pi y)\end{matrix}\right),\quad p=x^{2}y^{2}-\frac{1}{9}.

Consider the following permeability

\displaystyle\kappa^{-1}=a(sin({2\pi x})+1.1),

where $a$ is a given positive constant. According to the above parameters, the momentum source term ${\mathbf{f}}$ and the boundary value ${\mathbf{g}}={\mathbf{u}}|_{\partial\Omega}$ can be calculated.

For uniform triangular partition, we choose $k=1$ , $\mu=1$ , $0.01$ and $a=1$ , $10^{4}$ . Table 1-4 show the errors and orders of convergence accordingly. For uniform rectangular partition and polygonal partition, we choose $k=2$ , $3$ , $\mu=10^{4}$ and $a=1$ . Table 5-8 show the errors and orders of convergence accordingly.

As can be seen from the data in the tables, the numerical experiment results are consistent with the theoretical analysis, and both reach the optimal order of convergence. At the same time, the accuracy and stability of the numerical scheme is verified when the permeability $\kappa$ is highly varying.

Table 1: Errors and orders of convergence on triangular partition as

k=1,~{}j=2,~{}\mu=1,~{}a=1

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	1.5213e+00		1.02630e-01		4.3130e-01
$1/8$	6.5903e-01	1.2069	3.7485e-02	1.4531	3.1397e-01	0.4581
$1/16$	2.5867e-01	1.3493	1.0962e-02	1.7738	1.3347e-01	1.2341
$1/32$	1.0300e-01	1.3285	2.8870e-03	1.9249	4.7939e-02	1.4773
$1/64$	4.3307e-02	1.2499	7.3057e-04	1.9825	1.8019e-02	1.4116
$1/128$	1.9323e-02	1.1643	1.8291e-04	1.9979	7.5708e-03	1.2510

Table 2: Errors and orders of convergence on triangular partition as

k=1,~{}j=2,~{}\mu=0.01,~{}a=1

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	1.2093e+00		2.2360e-01		4.1153e-02
$1/8$	4.8033e-01	1.3321	7.8492-02	1.5103	1.8900e-02	1.1226
$1/16$	1.8275e-01	1.3941	2.2513e-02	1.8018	8.5580e-03	1.1430
$1/32$	7.1280e-02	1.3583	5.8547e-03	1.9431	3.9645e-03	1.1101
$1/64$	2.9236e-02	1.2858	1.4752e-03	1.9886	1.8795e-03	1.0768
$1/128$	1.2709e-02	1.2019	3.6902e-04	1.9992	9.1017e-04	1.0462

Table 3: Errors and orders of convergence on triangular partition as

k=1,~{}j=2,~{}\mu=1,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	2.2860e+00		2.9740e-02		1.6643e-01
$1/8$	3.6162e-01	2.6603	3.9858e-03	2.8995	1.9585e-01	-0.2349
$1/16$	1.5733e-01	1.2007	1.4717e-03	1.4374	1.8216e-01	0.1046
$1/32$	8.4025e-02	0.9049	5.4615e-04	1.4301	1.3765e-01	0.4042
$1/64$	4.1075e-02	1.0326	1.7638e-04	1.6306	7.6546e-02	0.8466
$1/128$	1.9215e-02	1.0960	5.0498e-05	1.8044	3.0815e-02	1.3127

Table 4: Errors and orders of convergence on triangular partition as

k=1,~{}j=2,~{}\mu=0.01,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	1.0369e+00		6.7825e-02		1.2495e-01
$1/8$	4.4656e-01	1.2154	2.0320e-02	1.7389	8.4347e-02	0.5670
$1/16$	1.7634e-01	1.3405	7.3711e-03	1.4630	4.7227e-02	0.8367
$1/32$	6.9966e-02	1.3336	2.4262e-03	1.6032	1.9975e-02	1.2414
$1/64$	2.9008e-02	1.2702	6.8638e-04	1.8216	6.4779e-03	1.6246
$1/128$	1.2675e-02	1.1945	1.7914e-04	1.9379	1.9103e-03	1.7618

Table 5: Errors and orders of convergence on rectangular partition as

k=2,~{}j=5,~{}\mu=1,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	1.4777e+00		1.7715e-02		5.4279e-01
$1/8$	2.5250e-01	2.5490	1.7706e-03	3.3227	1.2107e-01	2.1645
$1/16$	5.5071e-02	2.1969	2.4002e-04	2.8830	3.2960e-02	1.8771
$1/32$	1.2294e-02	2.1633	3.2057e-05	2.9045	6.2549e-03	2.3976
$1/64$	2.7958e-03	2.1367	3.8025e-06	3.0756	1.0068e-03	2.6352

Table 6: Errors and orders of convergence on rectangular partition as

k=3,~{}j=6,~{}\mu=1,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	6.4024e-01		5.3190e-03		4.8987e-01
$1/8$	7.5994e-02	3.0746	3.6500e-04	3.8652	4.0092e-02	3.6110
$1/16$	7.9936e-03	3.2490	2.2508e-05	4.0194	2.7254e-03	3.8788
$1/32$	8.2601e-04	3.2746	1.3259e-06	4.0854	1.8280e-04	3.8981
$1/64$	8.8561e-05	3.2214	7.8119e-08	4.0852	1.1620e-05	3.9756

Table 7: Errors and orders of convergence on polygonal partition as

k=2,~{}j=8,~{}\mu=1,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	1.3518e+00		1.5787e-02		5.9253e-01
$1/8$	3.4537e-01	1.9687	2.3238e-03	2.7642	1.8185e-01	1.7041
$1/16$	1.0534e-01	1.7130	3.8524e-04	2.5927	3.8610e-02	2.2357
$1/32$	2.7395e-02	1.9431	5.2055e-05	2.8876	4.9117e-03	2.9747
$1/64$	7.0084e-03	1.9668	6.6314e-06	2.9727	7.1932e-04	2.7715

Table 8: Errors and orders of convergence on polygonal partition as

k=3,~{}j=9,~{}\mu=1,~{}a=10^{4}

$h$	$\|\!\|\!\|{\mathbf{e}}_{h}\|\!\|\!\|$	order	$\\|{\mathbf{e}}_{h}\\|$	order	$\\|\varepsilon_{h}\\|$	order
$1/4$	6.2524e-01		5.2453e-03		4.0361e-01
$1/8$	7.9629e-02	2.9730	3.7527e-04	3.8050	3.5666e-02	3.5003
$1/16$	8.0684e-03	3.3029	2.4399e-05	3.9431	2.1461e-03	4.0548
$1/32$	8.2383e-04	3.2919	1.5631e-06	3.9644	1.5572e-04	3.7847
$1/64$	8.9162e-05	3.2078	9.5552e-08	4.0319	1.1355e-05	3.7775

The rest of examples in this section have the following setting:

\displaystyle k=1,\quad\Omega=(0,1)\times(0,1),\quad\mu=0.01,\quad{\mathbf{f}}=\left(\begin{matrix}0\\ 0\end{matrix}\right),\quad{\mathbf{g}}=\left(\begin{matrix}1\\ 0\end{matrix}\right).

6.2 Example 2

In this example, the permeability coefficient $\kappa$ is selected as the piecewise constant function with highly varying. The profile of the permeability inverse is shown in Fig. 1(a). As we know, this example has no analytic solutions.

In Fig. 1, a 128 $\times$ 128 rectangular partition is used to solve this example. The profiles of the pressure and the two components of the velocity for CDG are plotted in Fig. 1(b)-1(d)

Refer to caption — (a) Profile of $\kappa^{-1}$

6.3 Example 3

In this example, we choose a vuggy medium with the permeability coefficient $\kappa$ highly varying. The profile of the permeability inverse is plotted in Fig. 2(a). For solving this example, a 128 $\times$ 128 rectangular partition is used. And the pressure obtained by CDG method is present in Fig. 2(b). The velocity profiles are shown in Fig. 2(c)-2(d).

6.4 Example 4

The fluid flowing in a fibrous material can also be described by Brinkman equations. Fig. 3(a) shows the inverse of permeability in a common fibrous material. The parameters are the same as in the previous example. We can get the corresponding pressure in Fig. 3(b) and velocity in Fig. 3(c)-3(d) by CDG method.

Appendix A Some Inequality Estimates

In this Appendix, we provide inequalities for projection operators $Q_{h}$ , ${\mathbf{Q}}_{h}$ , ${\mathbb{Q}}_{h}$ and inequality estimates used in the previous paper.

Lemma A.1.

[22] Let ${\mathcal{T}}_{h}$ be a shape regular partition of $\Omega$ , ${\mathbf{w}}\in[H^{k+1}(\Omega)]^{d}$ and $\rho\in H^{k}(\Omega)$ . Then we have the following projection inequalities

$\displaystyle\sum_{T\in\mathcal{T}_{h}}\\|{\mathbf{w}}-Q_{h}{\mathbf{w}}\\|^{2}_{T}$	$\displaystyle\leq Ch^{2(k+1)}\\|{\mathbf{w}}\\|^{2}_{k+1},$	(A.1)
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\\|\nabla{\mathbf{w}}-{\mathbf{Q}}_{h}(\nabla{\mathbf{w}})\\|^{2}_{T}$	$\displaystyle\leq Ch^{2k}\\|{\mathbf{w}}\\|^{2}_{k+1},$	(A.2)
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\\|\rho-{\mathbb{Q}}_{h}\rho\\|^{2}_{T}$	$\displaystyle\leq Ch^{2k}\\|\rho\\|_{k}^{2}.$	(A.3)

where $C$ is a constant independent of the size of mesh $h$ and the functions ${\mathbf{w}}$ and $\rho$ .

Let $T$ be a cell with $e$ as an edge/face. For any function $\rho\in H^{1}(T)$ , the following trace inequality has been proved to be valid in [22]:

\displaystyle\|\rho\|^{2}_{e}\leq C(h^{-1}_{T}\|\rho\|^{2}_{T}+h_{T}\|\nabla\rho\|^{2}_{T}).

(A.4)

Furthermore, if $\rho\in P_{k}(T)$ , we have the following inverse inequality

\displaystyle\|\nabla\rho\|_{T}\leq Ch_{T}^{-1}\|\rho\|_{T}.

(A.5)

Lemma A.2.

For any ${\mathbf{v}}\in V_{h}$ , the following inequality holds true

\displaystyle\sum_{e\in{\mathcal{E}}_{h}}h^{-1}\|[{\mathbf{v}}]\|^{2}_{e}\leq C|\!|\!|{\mathbf{v}}|\!|\!|^{2},

(A.6)

where $C$ is a positive constant.

The proof of this lemma is given in Lemma 3.2 in [27].

Lemma A.3.

For any ${\mathbf{w}}\in[H^{k+1}(\Omega)]^{d}$ , $r\in H^{k}(\Omega)$ ${\mathbf{v}}\in V_{h}$ and $q\in W_{h}$ . Then we have

$\displaystyle\|l_{1}({\mathbf{w}},{\mathbf{v}})\|$	$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|,$	(A.7)
$\displaystyle\|l_{2}({\mathbf{w}},{\mathbf{v}})\|$	$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|,$	(A.8)
$\displaystyle\|l_{3}(r,{\mathbf{v}})\|$	$\displaystyle\leq Ch^{k}\\|r\\|_{k}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|,$	(A.9)
$\displaystyle\|l_{4}({\mathbf{w}},q)\|$	$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\\|q\\|_{h},$	(A.10)
$\displaystyle\|s({\mathbb{Q}}_{h}r,q)\|$	$\displaystyle\leq Ch^{k}\\|r\\|_{k}\\|q\\|_{h}.$	(A.11)

Proof.

Using the definition of $\nabla_{w}$ , the Cauchy-Schwarz inequality, the trace inequality (A.4) and the projection inequality (A.1), we obtain

	$\displaystyle\|l_{1}({\mathbf{w}},{\mathbf{v}})\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}(\nabla_{w}({\mathbf{w}}-Q_{h}{\mathbf{w}}),\nabla_{w}{\mathbf{v}})_{T}\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}(-({\mathbf{w}}-Q_{w}{\mathbf{w}},\nabla\cdot\nabla_{w}{\mathbf{v}})_{T}+\left<\left\{{\mathbf{w}}-Q_{w}{\mathbf{w}}\right\},\nabla_{w}{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T})\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}((\nabla({\mathbf{w}}-Q_{w}{\mathbf{w}}),\nabla_{w}{\mathbf{v}})_{T}+\left<Q_{w}{\mathbf{w}}-\left\{Q_{w}{\mathbf{w}}\right\},\nabla_{w}{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T})\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}((\nabla({\mathbf{w}}-Q_{w}{\mathbf{w}}),\nabla_{w}{\mathbf{v}})_{T}+\left<Q_{w}{\mathbf{w}}-{\mathbf{w}}+\left\{{\mathbf{w}}-Q_{w}{\mathbf{w}}\right\},\nabla_{w}{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T})\Bigg{\|}$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla({\mathbf{w}}-Q_{w}{\mathbf{w}})\\|_{T}\\|\nabla_{w}{\mathbf{v}}\\|_{T}+Ch^{-\frac{1}{2}}\\|{\mathbf{w}}-Q_{w}{\mathbf{w}}\\|_{\partial T}\\|\nabla_{w}{\mathbf{v}}\\|_{T})$
		$\displaystyle\leq\left(\sum_{T\in\mathcal{T}_{h}}(\\|\nabla({\mathbf{w}}-Q_{w}{\mathbf{w}})\\|_{T}+Ch^{-1}\\|{\mathbf{w}}-Q_{w}{\mathbf{w}}\\|_{T})^{2}\right)^{\frac{1}{2}}\left(\sum_{T\in\mathcal{T}_{h}}\\|\nabla_{w}{\mathbf{v}}\\|_{T}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|.$

Combining (A.6), (2.1) and the projection inequality (A.2) gives

	$\displaystyle\|l_{2}({\mathbf{w}},{\mathbf{v}})\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<(\nabla{\mathbf{w}}-{\mathbf{Q}}_{h}\nabla{\mathbf{w}})\cdot{\mathbf{n}},{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|\nabla{\mathbf{w}}-{\mathbf{Q}}_{h}\nabla{\mathbf{w}}\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{e\in{\mathcal{E}}_{h}}h^{-1}\\|[{\mathbf{v}}]\\|_{e}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|.$

According to the projection inequality (A.3), we get

	$\displaystyle\|l_{3}(r,{\mathbf{v}})\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<r-{\mathbb{Q}}_{h}r,{\mathbf{v}}\cdot{\mathbf{n}}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|r-{\mathbb{Q}}_{h}r\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{e\in{\mathcal{E}}_{h}}h^{-1}\\|[{\mathbf{v}}]\\|_{e}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k}\\|r\\|_{k}\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|.$

Similarly, with (2.2), then

	$\displaystyle\|l_{4}({\mathbf{w}},q)\|$	$\displaystyle=\Bigg{\|}\sum_{T\in\mathcal{T}_{h}}\left<({\mathbf{u}}-Q_{h}{\mathbf{u}})\cdot{\mathbf{n}},q-\left\{q\right\}\right>_{\partial T}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h^{-1}\\|{\mathbf{w}}-Q_{h}{\mathbf{w}}\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{e\in{\mathcal{E}}_{h}^{0}}h\\|[\![q]\!]\\|_{e}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k}\\|{\mathbf{w}}\\|_{k+1}\\|q\\|_{h}.$

Using the definition of $s(\cdot,\cdot)$ , we get

	$\displaystyle\|s({\mathbb{Q}}_{h}r,q)\|$	$\displaystyle=\Bigg{\|}\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}r]\!],[\![q]\!]\right>_{e}\Bigg{\|}$
		$\displaystyle=\Bigg{\|}\sum_{e\in{\mathcal{E}}_{h}^{0}}h\left<[\![{\mathbb{Q}}_{h}r-r]\!],[\![q]\!]\right>_{e}\Bigg{\|}$
		$\displaystyle\leq C\left(\sum_{T\in\mathcal{T}_{h}}h\\|r-{\mathbb{Q}}_{h}r\\|_{\partial T}^{2}\right)^{\frac{1}{2}}\left(\sum_{e\in{\mathcal{E}}_{h}^{0}}h\\|[\![q]\!]\\|_{e}^{2}\right)^{\frac{1}{2}}$
		$\displaystyle\leq Ch^{k}\\|r\\|_{k}\\|q\\|_{h},$

which completes the proof. ∎

Acknowledgments

This work was supported by National Natural Science Foundation of China (Grant No. 1901015, 12271208), and Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education, Jilin University. We sincerely thank the anonymous reviewers for their insightful comments, which have helped improve the quality of this paper.

References

[1] A. Al-Taweel, and X. Wang, A note on the optimal degree of the weak gradient of the stabilizer free weak Galerkin finite element method, Appl. Numer. Math., 150 (2020), 444-451.
[2] V. Anaya, G. N. Gatica, D. Mora and R. Ruiz-Baier, An augmented velocity-vorticity-pressure formulation for the Brinkman equations, Internat. J. Numer. Methods Fluids, 79 (2015), 109-137.
[3] A. A. Avramenko, I. V. Shevchuk, M. M. Kovetskaya, and Y. Y. Kovetska, Darcy-Brinkman-forchheimer model for film boiling in porous media, Transp. Porous Media, 134 (2020), 503-536.
[4] S. Badia, and R. Codina, Unified stabilized finite element formulations for the Stokes and the Darcy problems, SIAM J. Numer. Anal., 47 (2009), 1971-2000.
[5] H. Brinkman, A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, Flow Turbul. Combust., 1 (1949), 27-34.
[6] E. Cáceres, G. N. Gatica, and F. A. Sequeira, A mixed virtual element method for the Brinkman problem, Math. Models Methods Appl. Sci., 27 (2017), 707-743.
[7] B. Cockburn, G. Kanschat and D. Schötzau, A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations, J. Sci. Comput., 31 (2007), 61–73.
[8] C. D’Angelo and P. Zunino, A finite element method based on weighted interior penalties for heterogeneous incompressible flows, SIAM J. Numer. Anal., 47 (2009), 3990-4020.
[9] L. de Paulo Ferreira, T. D. S. de Oliveira, R. Surmas, M. A. P. da Silva, and R. P. Peçanha, Brinkman equation in reactive flow: Contribution of each term in carbonate acidification simulations, Adv. Water. Resour., 144 (2020), 103696.
[10] F. Golfier, D. Lasseux, and M. Quintard, Investigation of the effective permeability of vuggy or fractured porous media from a Darcy-Brinkman approach, Comput. Geosci., 19 (2015), 63-78.
[11] J. S. Howell, M. Neilan and N. J. Walkington, A dual-mixed finite element method for the Brinkman problem, SMAI J. Comput. Math., 2 (2016), 1-17.
[12] P. Huang and J. Chen, Some low order nonconforming mixed finite elements combined with Raviart-Thomas elements for a coupled Stokes-Darcy model, Jpn. J. Ind. Appl. Math., 30 (2013), 565–584.
[13] J. Könnö and R. Stenberg, H(div)-conforming finite elements for the Brinkman problem, Math. Models Methods Appl. Sci. , 21 (2011), 2227-2248.
[14] I. Ligaarden, M. Krotkiewski, K. A. Lie, M. Pal and D. W. Schmid, On the Stokes-Brinkman equations for modeling flow in carbonate reservoirs, Oxford: Proceedings of the the ECMOR XII - 12th European Conferenceon the mathematics of Oil Recovery, 2010.
[15] K. A. Mardal, X. C. Tai and R. Winther, A Robust Finite Element Method for Darcy - Stokes Flow, SIAM J. Numer. Anal., 40 (2002), 1605-1631.
[16] L. Mu, J. Wang, and X. Ye, A stable numerical algorithm for the Brinkman equations by weak Galerkin finite element methods, J. Comput. Phys., 273 (2014), 327-342.
[17] Y. Qian, S. Wu, and F. Wang, A mixed discontinuous Galerkin method with symmetric stress for Brinkman problem based on the velocity-pseudostress formulation, Comput. Methods Appl. Mech. Engrg., 368 (2020), 113177.
[18] K. Vafai, Porous media: applications in biological systems and biotechnology. Springer Series in Computational Mathematics, 2010, CRC Press, USA.
[19] G. Wang, Y. Wang, and Y. He, Two robust virtual element methods for the Brinkman equations, Calcolo, 58 (2021), 49.
[20] J. Wang and X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math., 241 (2013), 103-115.
[21] J. Wang and X. Ye, A weak Galerkin finite element method for the stokes equations, Adv. Comput. Math., 42 (2016), 155–174.
[22] J. Wang and X. Ye, A weak Galerkin mixed finite element method for second order elliptic problems, Math. Comp., 83 (2014), 2101-2126.
[23] Y. Wang, F. Gao, and J. Cui, A conforming discontinuous Galerkin finite element method for elliptic interface problems, J. Comput. Appl. Math., 412 (2022), 114304.
[24] Y. Wang, F. Gao, and J. Cui, A conforming discontinuous Galerkin finite element method for linear elasticity interface problems, J. Sci. Comput., 92 (2022), 20.
[25] X. Xie, J. Xu and G. Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, J. Comput. Math., 26 (2008), 437-455.
[26] X. Ye and S. Zhang, A conforming discontinuous Galerkin finite element method, Int. J. Numer. Anal. Model., 17 (2020), 110-117.
[27] X. Ye and S. Zhang, A conforming discontinuous Galerkin finite element method: Part II, Int. J. Numer. Anal. Model., 17 (2020), 281-296.
[28] X. Ye and S. Zhang, A conforming discontinuous Galerkin finite element method for the Stokes problem on polytopal meshes, Internat. J. Numer. Methods Fluids, 93 (2021), 1913-1928.
[29] Q. Zhai, R. Zhang, and L. Mu, A new weak Galerkin finite element scheme for the Brinkman model, Commun. Comput. Phys., 19 (2016), 1409-1434.
[30] X. Zhang and M. Feng, A projection-based stabilized virtual element method for the unsteady incompressible Brinkman equations, Appl. Math. Comp., 408 (2021), 126325.

	$\displaystyle\|\!\|\!\|q\|\!\|\!\|^{2}_{1}$	$\displaystyle=\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}+\sum_{e\in{\mathcal{E}}_{h}^{0}}h^{-1}\\|[\![q]\!]\\|_{e}^{2},$		(3.2)
	$\displaystyle\\|q\\|_{h}^{2}$	$\displaystyle=\sum_{e\in{\mathcal{E}}_{h}^{0}}h\\|[\![q]\!]\\|_{e}^{2}.$		(3.3)

	$\displaystyle\sum_{T\in\mathcal{T}_{h}}\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}^{2}$
	$\displaystyle\leq\sum_{e\in{\mathcal{E}}_{h}}\left(\left\\|{\mathbf{v}}\|_{T_{1}}-\frac{{\mathbf{v}}\|_{T_{1}}+{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{\partial T_{1}\cap e}^{2}+\left\\|{\mathbf{v}}\|_{T_{2}}-\frac{{\mathbf{v}}\|_{T_{1}}+{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{\partial T_{2}\cap e}^{2}\right)$
	$\displaystyle\leq C\sum_{e\in{\mathcal{E}}_{h}}\left\\|\frac{{\mathbf{v}}\|_{T_{1}}-{\mathbf{v}}\|_{T_{2}}}{2}\right\\|_{e}^{2}$
	$\displaystyle\leq C\sum_{e\in{\mathcal{E}}_{h}}\left(\\|{\mathbf{v}}\\|_{\partial T_{1}\cap e}^{2}+\\|{\mathbf{v}}\\|_{\partial T_{2}\cap e}^{2}\right)$
	$\displaystyle\leq C\sum_{T\in\mathcal{T}_{h}}\\|{\mathbf{v}}\\|_{\partial T}^{2}.$

	$\displaystyle\\|\nabla_{w}{\mathbf{v}}\\|^{2}$	$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}\\|\nabla_{w}{\mathbf{v}}+\nabla{\mathbf{v}}\\|_{\partial T})$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+h^{-\frac{1}{2}}\\|{\mathbf{v}}-\left\{{\mathbf{v}}\right\}\\|_{\partial T}(\\|\nabla_{w}{\mathbf{v}}\\|_{T}+\\|\nabla{\mathbf{v}}\\|_{T}))$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+Ch^{-\frac{1}{2}}\\|{\mathbf{v}}\\|_{\partial T}(\\|\nabla_{w}{\mathbf{v}}\\|_{T}+\\|\nabla{\mathbf{v}}\\|_{T}))$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}(\\|\nabla{\mathbf{v}}\\|_{T}^{2}+\frac{1}{2}\\|\nabla_{w}{\mathbf{v}}\\|_{T}^{2}+\frac{1}{2}\\|\nabla{\mathbf{v}}\\|_{T}^{2}+Ch^{-1}\\|{\mathbf{v}}\\|_{\partial T}^{2})$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\left(\frac{1}{2}\\|\nabla_{w}{\mathbf{v}}\\|_{T}^{2}+Ch^{-2}\\|{\mathbf{v}}\\|_{T}^{2}\right).$

	$\displaystyle\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|^{2}$	$\displaystyle=\\|\nabla_{w}{\mathbf{v}}\\|^{2}+\\|\kappa^{-\frac{1}{2}}{\mathbf{v}}\\|^{2}$
		$\displaystyle\leq Ch^{-2}\\|{\mathbf{v}}\\|^{2}+\\|\kappa^{-\frac{1}{2}}{\mathbf{v}}\\|^{2}$
		$\displaystyle\leq Ch^{-2}\\|\kappa\tilde{\nabla}_{w}q\\|^{2}+\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}$
		$\displaystyle\leq Ch^{-2}\lambda_{1}^{-1}\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}+\\|\kappa^{\frac{1}{2}}\tilde{\nabla}_{w}q\\|^{2}$
		$\displaystyle\leq Ch^{-2}\|\!\|\!\|q\|\!\|\!\|_{1}^{2}.$

	$\displaystyle\frac{b({\mathbf{v}},q)}{\|\!\|\!\|{\mathbf{v}}\|\!\|\!\|}$	$\displaystyle=\frac{\|\!\|\!\|q\|\!\|\!\|_{1}^{2}-h^{-2}\\|q\\|_{h}^{2}}{Ch^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}}$
		$\displaystyle\geq\frac{\|\!\|\!\|q\|\!\|\!\|_{1}^{2}-h^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}\\|q\\|_{h}}{Ch^{-1}\|\!\|\!\|q\|\!\|\!\|_{1}}$
		$\displaystyle\geq C_{1}h\|\!\|\!\|q\|\!\|\!\|_{1}-C_{2}\\|q\\|_{h}.$