Hybridized weak Galerkin finite element methods for Brinkman equations

1.   Introduction

The Brinkman equation describes the problem of fluid motion in porous media and is an appropriate model for fluid motion in higher-order non-uniform media. The model can also be seen as a generalization of the Stokes equation, that is, an effective approximation of the Navier-Stokes equation at low Reynolds numbers. Simulating fluid flow in a composite medium with multiphysics effects has significant impacts on many industrial and environmental problems, such as drilling, channels and fluid flow near faults. The permeability with high contrast determines that the flow rate through porous media can vary greatly. Mathematically, the Brinkman equation can be regarded as the combination of the Stokes equation and the Darcy equation, either of which dominantly appear in different area of the domain depending on its characteristic. Due to the change of type, the numerical algorithm [7] for solving the Brinkman equation must be able to handle both the Stokes and the Darcy equation. The numerical experiments in [5,8] show that when a fixed Stokes element is selected, the Brinkman equation is controlled by Darcy and the convergence rate is reduced. Similarly, when a fixed Darcy element is selected, the Brinkman equation becomes controlled by Stokes and the rate of convergence will also be reduced accordingly. That is, the usual Stokes stable elements are not suitable for Darcy fluids and vice versa. At present, the method for solving the Brinkman equation [6,25,28] has a finite volume discrete method, non-coordinated finite element method, weak Galerkin finite element method, etc. This article will introduce a stable and accurate calculation method for the Stokes and Darcy fluid regions, namely the hybrid weak finite element method.

In 2011, weak Galerkin finite element methods [18,24,26,32], referred to as WG method, was proposed by Junping Wang and Xiu Ye. It is a common finite element method for solving partial differential equations. It has played an important role in many fields, such as physics, biology and geosciences [2]. At the same time, the basic theory of mathematics has been improved and the method has become the research project of many researchers and engineers of computational mathematics. Its main characteristics are: (1) differential operators are approximated by discrete weak form; (2) the weak continuity of numerical solution is achieved by introducing stabilizer. The subdivision element of WG can be any polyhedron, and its approximation function space is composed of discontinuous piecewise polynomials. The flexibility of WG in the selection of approximation polynomials makes it as an ideal choice for the stable numerical scheme of partial differential equations with multiple physical properties. In addition, WG has been widely used to solve a variety of partial differential equations, such as the second-order elliptic equation [9,17,22,21,30], Maxwell equation [14], Stokes equation [20,23,27], Brinkman equation [11], and biharmonic equation [10,13,31].

In order to reduce the requirement of the continuity of numerical solution, hybrid technique [4,12,29] was introduced. It has been used as an effective way to solve partial differential equations. For example, HWG reduces the requirement of the continuity of piecewise polynomials in the whole region in the weak finite element method by introducing Lagrange multipliers at the boundary of each subdivision element. As such, it makes its construction simple, highly flexible and efficient.

The aim of this paper is to apply HWG to solve the Brinkman equation, and uses the Schur complement technique to reduce its degree of freedom, so as to improve the calculation efficiency. We shall show that the Schur complement formulation is well-posed. More specifically, we shall apply HWG to solve the Brinkman equation with the following three different boundary conditions:

$(1)$ Brinkman equation under Dirichlet boundary condition

$(2)$ Brinkman equation under Neumann boundary condition

$(3)$ Brinkman equation under Robin boundary condition

where $\mu$ is the viscosity of the fluid, $\kappa$ represents the permeability tensor of a polygon or polyhedron region, $\Omega \in \mathbb{R}^d$ with $(d = 2,3)$, ${\mathit{\boldsymbol{u}}}$ and $p$ represent the velocity and pressure of the fluid respectively, ${\mathit{\boldsymbol{f}}}, \mathit{\boldsymbol{\gamma}}$ and $\mathit{\boldsymbol{\theta}}$ are the source terms, $\alpha > 0$ is a parameter, and ${\mathit{\boldsymbol{n}}}$ is the unit outward normal vector to $\partial \Omega$.

The rest of this paper is organized as follows. In Section 2, we introduce notation for the Sobolev or broken Sobolev spaces, some inequalities, and the concepts of weak gradient and weak divergence. In Section 3, we introduce the HWG finite element method to solve the Brinkman equation under the Dirichlet boundary condition and establish the well-posedness and stability of the numerical solution. We also present error estimates in $H^1$ and $L^2$ norms. The Schur complement technique is then introduced to improve the algorithm. Section 4 describes the numerical algorithm and theoretical analysis of the HWG method for Brinkman equation with Neumann boundary condition. The Robin boundary case is discussed in Section 5. Numerical experiments are then presented to confirm the theoretical analysis in Section 6.

2.   Notation

We let $\Omega \subset \mathbb{R}^d$ be polygonal for $d = 2$ or polyhedral domain. Let $\mathcal{T}_h$ be a finite element partition, which satisfies the shape regular assumption [21]. We then denote all the edges of $\mathcal{T}_h$ by $\mathcal{E}_h$ and all the interior edges by $\mathcal{E}_h^0 = \mathcal{E}_h\backslash\partial\Omega$. We let $h = \max_{T \in \mathcal{T}_h} h_T$, where $h_T$ denotes the diameter of $T$.

On each $T\in \mathcal{T}_h$, we define the weak function spaces $\mathcal{V}(T)$, $V(T)$ by

where ${\mathit{\boldsymbol{n}}}$ is the outward normal direction to $\partial \Omega$. We then define the function space on $\mathcal{T}_h$ and $\mathcal{E}_h$, denoted by $\mathcal{V}$ and $\Lambda$, respectively as follows

For any $e \in \mathcal{E}_h$, we define the jump of both $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}$ and $q$ as follows

For any $e \in \mathcal{E}_h$, we now define the similarity of ${\mathit{\boldsymbol{\lambda}}} \in \Lambda$ as follows

Let $K$ be either $T \in \mathcal{T}_h$ or $e \in \mathcal{E}_h$ and denote the space of polynomial of degree less than or equal to $\ell$ by $P_\ell(K)$. For $T \in \mathcal{T}_h$, we define the discrete analogue of weak function spaces of $\mathcal{V}(T)$ and $V(T)$, denoted by $V_k(T)$ and $V_{k,N}(T)$, respectively as follows:

where $k\geq 1$ is a constant. For $T \in \mathcal{T}_h$, we also define $W_k(T)$ and $\Lambda_k(\partial T)$, respectively, by

We then define the weak finite element function spaces $V_h$, $\Lambda_h$ and $W_h$ as follows:

We shall also consider the subspaces of $V_h$ and $\Lambda_h$. First, we define $V_h^0, \mathcal{V}_h, \mathcal{V}_h^0 \subset V_h$, respectively by

Secondly, we define $\Xi_h \subset \Lambda_h$ as follows:

The space $\Xi_h$ will be taken as Lagrange multiplier approximation space for HWG. For $T \in \mathcal{T}_h$, we shall let $(\cdot,\cdot)_T$ and $\langle \cdot, \cdot \rangle_{\partial T}$ denote the standard $L^2$ inner product on $T$ and $\partial T$, respectively. We are now in a position to introduce a couple of bilinear forms for any given $T \in \mathcal{T}_h$: for $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}, \mathit{\boldsymbol{w}} = \{\mathit{\boldsymbol{w}}_0,\mathit{\boldsymbol{w}}_b\} \in V_k(T)$, $q \in W_k(T)$, $\mathit{\boldsymbol{\lambda}}\in\Lambda_{k}(\partial T)$,

where $\nabla_w \mathit{\boldsymbol{v}}$ and $\nabla_w \cdot \mathit{\boldsymbol{v}}$ are weakly defined gradient and divergence operator in Definition 2.5 and 2.6.

We then define the bilinear forms under different boundary conditions by summing bilinear forms defined locally above, by the following:

We introduce a couple of norms for the space $\mathcal{V}_h$, $\Xi_h$, and $V_h^0$ as follows

Definition 2.1. ([29]) For any $\mathit{\boldsymbol{v}}\in \mathcal{V}_{h}$, we let

where $\|\cdot\|$ is the standard $L^2$ norm on $\Omega$ and $\|\cdot\|_{\partial T}$ is the $L^2$ norm on $\partial T$.

Definition 2.2. ([29]) For $\mathit{\boldsymbol{\lambda}}\in\Xi_{h}$, let

where $h_e$ is the diameter of the edge/face $e \in \mathcal{E}_h$ and $\|\cdot\|_e$ is the $L^2$ norm on $e$.

Definition 2.3. ([29]) For $\mathit{\boldsymbol{v}}\in V_{h}^0$, let

Definition 2.4. For $\mathit{\boldsymbol{v}}\in V_{h,N}$, let

Definition 2.5. For $\mathit{\boldsymbol{v}}\in V_{h}$, let

For any given element $T \in \mathcal {T}_h$ and each edge/face $e\in\mathcal{E}_h$, let $Q_{0}$ and $Q_{b}$ be the $L^{2}$ projection operator from $[L^{2}(T)]^d$ to $[P_{k}(T)]^d$ and from $[L^{2}(e)]^d$ to $[P_{k}(e)]^d$, respectively. Let $\mathit{\boldsymbol{Q}}_{h}$ and $\mathbb{Q}_h$ be the orthogonal $L^{2}$ projection operator from $[L^{2}(T)]^{d\times d}$ to $[P_{k-1}(T)]^{d\times d}$ and from $L^{2}(T)$ to $P_{k-1}(T)$, respectively.

Lastly, following [17], we shall introduce discrete weak gradient and divergence. We begin with the definition of discrete weak gradient as follows:

Definition 2.6. (Discrete weak gradient) For any $\mathit{\boldsymbol{v}} \in\mathcal{V}(T)$, denote the discrete weak gradient operator $\nabla_{w,r,T}{\mathit{\boldsymbol{v}}}$ of $\mathit{\boldsymbol{v}}$ as the unique polynomial in $[P_r(T)]^{d\times d}$ such that for any $\tau\in[P_r(T)]^{d\times d}$, it satisfies

Definition 2.7. ([17]) (Discrete weak divergence) For any $\mathit{\boldsymbol{v}}\in V(T)$, denote the discrete weak divergence operator $\nabla_{w,r,T}\cdot{\mathit{\boldsymbol{v}}}$ of $\mathit{\boldsymbol{v}}$ as the unique polynomial in $P_r(T)$, such that for any $\varphi\in P_r(T)$, it satisfies

From the definition, we notice that the following identities hold: $\forall \mathit{\boldsymbol{v}} \in \mathcal{V}(T)$ and $\tau \in [P_r(T)]^{d\times d}$,

and $\forall \mathit{\boldsymbol{v}} \in V(T)$ and $\varphi \in P_r(T)$,

Denote by $\nabla_{w,k}$ and $\nabla_{w,k-1}\cdot$ the discrete weak gradient operator and the discrete weak divergence operator on the finite element space, which can be computed by using (4) and (5) on each element T, respectively; i.e.,

For simplicity of notation, we shall drop the subscript $k$ and $k-1$ in the notation of $\nabla_{w,k}$ and $(\nabla_{w,k-1}\cdot)$, respectively.

3.   HWG for Brinkman equation with Dirichlet boundary condition (1)

In this section, we present HWG algorithm to solve Brinkman equation with Dirichlet boundary condition (1).

3.1.   Algorithm

The following is the weak Galerkin (WG) finite element numerical scheme of Brinkman first variational formulation [17],

Algorithm 3.1. We seek $(\bar{\mathit{\boldsymbol{u}}}_h;\bar{p}_h)\in V_h\times W_h$ with $\bar{\mathit{\boldsymbol{u}}}_b = Q_b \mathit{\boldsymbol{g}}$ on $\partial\Omega$, such that

for all $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}\in \mathcal V_h$ and $q\in W_h$.

We now present the HWG method for (1). HWG method is attained by introducing the Lagrange multiplier to relax on the boundary of each inner element. Namely, it can be formulated as follows (see [29] for Stokes equation):

Algorithm 3.2. We seek $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_{h}$, with ${\mathit{\boldsymbol{u}}}_b = Q_b \mathit{\boldsymbol{g}}$ on $\partial\Omega$, such that

for all $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}\in V_h^0$, $q\in W_h$, and ${\mathit{\boldsymbol{\mu}}}\in\Xi_{h}$.

We shall establish that the problem (9) is well-posed.

Lemma 3.1. There exists a unique solution to (9).

Proof. Since (9) is linear, we only need to consider the uniqueness of homogeneous equation, let ${\mathit{\boldsymbol{f}}} = \mathit{\boldsymbol{0}}$, $\mathit{\boldsymbol{v}} = {\mathit{\boldsymbol{u}}}_h$, $q = p_h$, ${\mathit{\boldsymbol{\mu}}} = \mathit{\boldsymbol{\lambda}}_h$, then

With the definition of $a(\cdot,\cdot)$ , for any $T\in\mathcal{T}_h$, we have $\nabla_w {\mathit{\boldsymbol{u}}}_h = 0$, ${\mathit{\boldsymbol{u}}}_0 = 0$, ${\mathit{\boldsymbol{u}}}_0 = {\mathit{\boldsymbol{u}}}_b$ on $\partial T$.

Take any $\tau \in [P_{k-1}(T)]^{d\times d}$, according to (6), we have

Then for any $T\in\mathcal{T}_h,\ \nabla{\mathit{\boldsymbol{u}}}_0 = 0$. That is, for any $\partial T,\ {\mathit{\boldsymbol{u}}}_0 = {\mathit{\boldsymbol{u}}}_b = \mathit{\boldsymbol{0}}$. Let $\mathit{\boldsymbol{v}}_b = \mathit{\boldsymbol{0}}$, according to ${\mathit{\boldsymbol{u}}}_h = \{\mathit{\boldsymbol{0}}, \mathit{\boldsymbol{0}}\}$ we have

That is, for all $T\in\mathcal{T}_h,\ \nabla p_h = \mathit{\boldsymbol{0}}$.

For any two adjacent elements $T_1$ and $T_2$ with the common edge $e$, take $\mathit{\boldsymbol{v}}_b|_{e,T_1} = [\![p_h]\!]_e$, $\mathit{\boldsymbol{v}}_b|_{e,T_2} = -[\![p_h]\!]_e$; the same, take $\mathit{\boldsymbol{v}}_b|_e = [\![{\mathit{\boldsymbol{u}}}_h] \!]_e$, and in $\Omega$ , $\mathit{\boldsymbol{v}}_0 = 0$, we have

Since $0 = b(\mathit{\boldsymbol{v}},p_h) = \sum_{e\in\mathcal{\varepsilon}_h}\|[\![p_h]\!]\|_e^2$, we notice that $p_h$ is a constant. Furthermore, since $p_h \in L^2_0(\Omega)$, $p_h = 0$ in $\Omega$. Lastly, let $\mathit{\boldsymbol{v}}_b = \mathit{\boldsymbol{\lambda}}_h$, then

and therefore $\mathit{\boldsymbol{\lambda}}_h = \mathit{\boldsymbol{0}}$. This completes the proof.

Theorem 3.2. We assume that ${\mathit{\boldsymbol{u}}}_h\in V_h$ is the solution to HWG algorithm (9), then ${\mathit{\boldsymbol{u}}}_h$ is the solution of WG algorithm (8).

Proof. For $e \in \mathcal{E}_h^0$ with $\partial T_1 \cap \partial T_2 = e$, let ${\mathit{\boldsymbol{\mu}}} = [\![{\mathit{\boldsymbol{u}}}_h]\!]_{e}$ on $\partial T_1\cap e$ , ${\mathit{\boldsymbol{\mu}}} = -[\![{\mathit{\boldsymbol{u}}}_h]\!]_{e}$ on $\partial T_2\cap e$, and ${\mathit{\boldsymbol{\mu}}} = 0$ elsewhere. According to (9), we have

This leads that $[\![{\mathit{\boldsymbol{u}}}_h]\!]_{e} = \mathit{\boldsymbol{0}}, \forall e \in \mathcal{E}_h^0$. Now, by taking ${\mathit{\boldsymbol{\mu}}} = \mathit{\boldsymbol{0}}$, we have $b({\mathit{\boldsymbol{u}}}_h,q) = 0$. For all $\mathit{\boldsymbol{v}}\in V_h$, take $[\![\mathit{\boldsymbol{v}}]\!]_e = \mathit{\boldsymbol{0}},\ \forall e\in\mathcal{\varepsilon}_h^0$ and $\mathit{\boldsymbol{v}}|_{\partial\Omega} = \mathit{\boldsymbol{0}}$, we derive

This completes the proof.

3.2.   Stability analysis

Lemma 3.3. (Boundedness) There exists a constant $C>0$ such that

Proof. For (10), according to the definition of $a(\cdot,\cdot)$ and Cauchy-Schwarz inequality, we can have

For (11), according to the definition of $b(\cdot,\cdot)$, (7), Cauchy-Schwarz inequality, and trace inequality, we can have

For (12), we invoke the definition of $c(\cdot,\cdot)$ and Cauchy-Schwarz inequality to obtain

This completes the proof. We now establish the coercivity:

Lemma 3.4. We have that

Proof. Since $\forall \mathit{\boldsymbol{v}}\in \mathcal V_h$, it holds that $\|\mathit{\boldsymbol{v}}\|^2_{V_h^0} = |\!|\!|\mathit{\boldsymbol{v}}|\!|\!|^2$. This completes the proof. We shall now establish total three inf-sup conditions.

Lemma 3.5. (inf-sup condition 1) There is a constant $\beta>0$ independent of $h$ such that for any $\rho \in W_h$, we have

Proof. $\forall\rho\in W_h\subset L^2_0(\Omega)$, there is $\tilde{\mathit{\boldsymbol{v}}}\in [H^1_0(\Omega)]^d$ and $C>0$, such that

For $\mathit{\boldsymbol{v}} = Q_h\tilde{\mathit{\boldsymbol{v}}}\in V_h$, $|\!|\!|\mathit{\boldsymbol{v}}|\!|\!|\leq C_0\|\tilde{\mathit{\boldsymbol{v}}}\|_1$. According to the definition of norm and trace inequality, we have

Now due to the identity:

We complete the proof.

Lemma 3.6. (inf-sup condition 1') For any $\rho\in W_h$, there is a constant $\beta>0$ independent of $h$ and $\mathit{\boldsymbol{v}}\in V_h^0$ such that

Proof. $\forall \rho\in W_h\subset L_0^2(\Omega)$, there exists $\tilde{\mathit{\boldsymbol{v}}}\in[H_0^1(\Omega)]^d$ and $C>0$ making

We want to prove $\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert v \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert\leq C_0\Vert\tilde{\mathit{\boldsymbol{v}}}\Vert_1$ with $\mathit{\boldsymbol{v}} = Q_h\tilde{\mathit{\boldsymbol{v}}}\in V_h$. It follows from the definition of norm, trace inequality, and inverse inequality that

Then

This completes the proof.

Lemma 3.7. (inf-sup condition 2) There is a constant $C>0$, for any given $\mathit{\boldsymbol{\tau}}\in \Xi_{h}$, there is $\mathit{\boldsymbol{v}}\in V_h^0,\mathit{\boldsymbol{v}}_0 = \mathit{\boldsymbol{0}}$, so that

Proof. $\forall\mathit{\boldsymbol{\tau}}\in\Xi_{h}, \langle\!\langle\mathit{\boldsymbol{\tau}}\rangle\!\rangle_{e} = \mathit{\boldsymbol{0}}$. Let $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{0}},h_e\mathit{\boldsymbol{\tau}}\}\in V_{h}^0$, according to the definition of bilinear form, we have

where $\mathit{\boldsymbol{v}}_b^i$ and $\mathit{\boldsymbol{\tau}}^i \; (i = 1,2)$ represent the value of $\mathit{\boldsymbol{v}}_b|_{T_i}$ and $\mathit{\boldsymbol{\tau}}|_{T_i}$, respectively. Using Cauchy-Schwarz inequality, trace inequality, and inverse inequality,

where $\mathit{\boldsymbol{v}}_b^{*}\,$ can be selected as $\mathit{\boldsymbol{v}}_b^1\,$ or $\,\mathit{\boldsymbol{v}}_b^2$ and $\mathit{\boldsymbol{\tau}}^{*}\,$ can be selected as $\mathit{\boldsymbol{\tau}}^1\,$ or $\,\mathit{\boldsymbol{\tau}}^2$, it depends on the sectioning unit $\mathit{\boldsymbol{v}}_b\,$ and $\tau$ is in. As a result of $\|\nabla_w \mathit{\boldsymbol{v}}\|_T\leq C\sum_{e\in\partial T}h_e^{\frac{1}{2}}\|\mathit{\boldsymbol{\tau}}^{*}\|_e$ , we can get

This completes the proof.

3.3.   Error equation

The purpose of this section is to construct the error equation [3,19] between the numerical solution and the true solution of HWG according to the numerical algorithm (9).

Now, we shall present the properties of projection operators without proofs: (see [17] for proofs).

Lemma 3.8. The projection operators $Q_{h}$, $\mathbb{Q}_{h}$ and $\mathit{\boldsymbol{Q}}_{h}$ satisfy the following properties:

Lemma 3.9. Assume that $({\mathit{\boldsymbol{u}}};p) \in[H_0^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution of (1), $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_h$ is the solution of (9). Let ${\mathit{\boldsymbol{\lambda}}} = \nabla {\mathit{\boldsymbol{u}}} \cdot {\mathit{\boldsymbol{n}}}-p {\mathit{\boldsymbol{n}}}$, take $\mathit{\boldsymbol{e}}_{h} = \{Q_0 {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_0, Q_b {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_b\},\ \varepsilon_{h} = \mathbb{Q}_h p-p_{h},\ {\mathit{\boldsymbol{\delta}}}_{h} = Q_{b}\mathit{\boldsymbol{\lambda}}-{\mathit{\boldsymbol{\lambda}}}_{h}$. Then the error function $\mathit{\boldsymbol{e}}_{h}$, $\varepsilon_{h}$, and ${\mathit{\boldsymbol{\delta}}}_{h}\,$ satisfy the following equations

where

Proof. First, we invoke the definition of discrete weak gradient (4) and partial integral,

By adding these for all $T \in \mathcal{T}_h$, we obtain

Similarly, from (5) and partial integral, we have that

Hence,

Testing $\mathit{\boldsymbol{v}}_0\,$ for both sides of (1), we obtain

Now, from the identity:

we can have

By combining these with (9), we have that

which then results in

From Theorem 3.2, $[\![\mathit{\boldsymbol{e}}_{h}]\!]_{e} = \mathit{\boldsymbol{0}}$, we can get $c(\mathit{\boldsymbol{e}}_h, {\mathit{\boldsymbol{\mu}}}) = 0,\ \forall {\mathit{\boldsymbol{\mu}}}\in\Xi_{h}.$

Now, for any $q \in W_h$, $b(\mathit{\boldsymbol{e}}_h,q) = 0\,$ and we can get

This completes the proof.

3.4.   Error estimation

In this section, we establish the $H^1$ and $L^2$ norm error estimates using the error equations (13)-(14). To do so, we first provide some simple, but useful lemmas.

Lemma 3.10. If $({\mathit{\boldsymbol{u}}};p) \in[H_0^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution to the problem (1), there is a constant $C\,$ such that

Proof. Using Cauchy-Schwarz inequality, trace inequality, and inverse inequality, we have

Same as the proof of (17), according to Cauchy-Schwarz inequality and trace inequality, we have

By the nature of $Q_b\,$, Cauchy-Schwarz inequality, trace inequality, and inverse inequality, we can get

The theorem is proved.

Theorem 3.11. Assume that $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution satisfying (1), $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h \times \Xi_h$ is the solution of (9), then

Proof. In the error equation (13)-(14), taking $\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{e}}_h, {\mathit{\boldsymbol{\mu}}} = {\mathit{\boldsymbol{\delta}}}_{h}, q = \varepsilon_h\,$, we have

In (16), let $\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{e}}_h$, we have

According to $|\mathit{\boldsymbol{e}}_h|_h\leq C|\!|\!| \mathit{\boldsymbol{e}}_h |\!|\!|$, we can further obtain

There are the following facts $|\!|\!| \mathit{\boldsymbol{e}}_h |\!|\!| = \|\mathit{\boldsymbol{e}}_h\|_{V_{h}^0}$, so

According to Lemma 3.5, by taking $\mathit{\boldsymbol{v}}^* = \{0,\mathit{\boldsymbol{v}}_b\}$, we have

According to (16), error equation (13), and boundedness of bilinear form (10)-(11), we can get

And because $b(\mathit{\boldsymbol{v}}^*,\varepsilon_h)\geq\beta\|\varepsilon_h\|_{\Xi_{h}}|\!|\!|\mathit{\boldsymbol{v}}^*|\!|\!|\,$, we have

Taking $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{0}},\mathit{\boldsymbol{v}}_b\}\,$, same as the proof of (18), we can get

And because $c(\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{\delta}}}_h)\geq C\|{\mathit{\boldsymbol{\delta}}}_h\|_{\Xi_{h}}|\!|\!|\mathit{\boldsymbol{v}}|\!|\!|\,$, we have

This completes the proof.

Finally, the dual technique is used to derive the optimal order error estimates of the WG scheme under $L^2$ norm. Consider the following dual problems

with $(\psi;\xi)\in [H^2(\Omega)]^d\times H^1(\Omega)$. Assume that the dual problem is $H^2$ -regular, that is, the constant $C$ makes

Theorem 3.12. Suppose $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution to the problem (1), $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times \Xi_h$ is the solution of (9), then

Proof. Multiplying $\mathit{\boldsymbol{e}}_0\,$ to both sides of (19) gives

Take ${\mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{\psi}},\,\mathit{\boldsymbol{v}}_0 = \mathit{\boldsymbol{e}}_0,\,p = \xi\,$ in the above formula, from the error equation

The following estimates the above formula item by item

Because of the following fact

we can get $|\phi_{{\mathit{\boldsymbol{u}}},p}(Q_h\psi)|\leq Ch^{k+1}(\|{\mathit{\boldsymbol{u}}}\|_{k+1}+\|p\|_k)\|\psi\|_2\,$.

In (16), by taking ${\mathit{\boldsymbol{u}}} = \psi, \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{e}}_h, p = \xi$, we can get

Then

From regularity (20), we have

This completes the proof.

Note that Under the condition of Dirichlet boundary value, change the space of Lagrange multiplier and redefine it as

Denote $\widetilde{Q}_{b}$ the $L^{2}\,$ projection operator from $[L^{2}(e)]^d$ to $[P_{k-1}(e)]^d\,$. Then from (15), we can get

The error equation is

and $c(\mathit{\boldsymbol{v}},\widetilde{Q}_b \mathit{\boldsymbol{\lambda}}-\mathit{\boldsymbol{\lambda}}) = \langle \mathit{\boldsymbol{v}}_b,(\nabla{\mathit{\boldsymbol{u}}}\cdot {\mathit{\boldsymbol{n}}}-p{\mathit{\boldsymbol{n}}})-\widetilde{Q}_b(\nabla{\mathit{\boldsymbol{u}}}\cdot {\mathit{\boldsymbol{n}}}-p{\mathit{\boldsymbol{n}}})\rangle_{\partial T} = 0\,$, so the error equation is the same as Theorem 3.9, we can get the same error estimates.

3.5.   Theoretical analysis of Schur complement method

Due to the introduction of Lagrange multipliers, the number of unknowns to be solved is increased in HWG method. The purpose of this section is to apply $Schur$ complement technique [24,29] to reduce degrees of freedom, based on the numerical scheme constructed by HWG method. That is, boundary function ${\mathit{\boldsymbol{u}}}_b$ is used to express internal function ${\mathit{\boldsymbol{u}}}_0$ and Lagrange multiplier $\mathit{\boldsymbol{\lambda}}_h$.

First, we define the boundary finite element space $B_h$ as follows

For Hilbert space $B_h$, we define inner product as follows

$B_h^0$ is a subspace of $B_h$, consisting of functions in $B_h$, with zero boundary value. Obviously, $B_h$ is isomorphic to $\Xi_{h}$. In order to eliminate Lagrange multiplier $\mathit{\boldsymbol{\lambda}}_h$ and interior unknowns ${\mathit{\boldsymbol{u}}}_0$ by Schur complement technique, we introduce mapping $S_{{\mathit{\boldsymbol{f}}}}:B_h\rightarrow B_h^0$.

For a fixed function $p_h$ and any given function $\mathit{\boldsymbol{\omega}}_b\in B_h$, we shall define $S_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{\omega}}_b; p_h)$ by the following three steps:

Step 1: On each element $T\in\mathcal T_h$, $\mathit{\boldsymbol{\omega}}_0$ is represented by $\mathit{\boldsymbol{\omega}}_b$ and $p_h$ through the following equation:

where $\mathit{\boldsymbol{\omega}}_h = \{\mathit{\boldsymbol{\omega}}_0,\mathit{\boldsymbol{\omega}}_b\}\in V_k(T)$, $p_h\in W_k(T)$. Then we can work out $\mathit{\boldsymbol{\omega}}_0 = D_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{\omega}}_b;p_h)$ from (21).

Step 2: On each element $T\in \mathcal T_h$, we represent $\mathit{\boldsymbol{\zeta}}_{h,T}\in \Lambda_k(\partial T)$ by $\mathit{\boldsymbol{\omega}}_h = \{\mathit{\boldsymbol{\omega}}_0,\mathit{\boldsymbol{\omega}}_b\}\in V_k(T)$ and $p_h$

Then we can work out $\mathit{\boldsymbol{\zeta}}_{h,T}\in \Lambda_h,\ \mathit{\boldsymbol{\zeta}}_{h,T} = L_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{\omega}}_b;p_h)$ from (22).

Step 3: We then define $S_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{\omega}}_b;p_h)$ by the following: the similarity of $\mathit{\boldsymbol{\zeta}}_h$ on the inner boundary and $\mathit{\boldsymbol{0}}$ on the outer boundary, that is

We observe that by (23), $S_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{\omega}}_b;p_h) \in B_h^0$. Furthermore, the operator $S_{{\mathit{\boldsymbol{f}}}}$ has the following properties:

(1) Summing (21) and (22), we obtain that

(2) From the superposition principle, we have that

where $S_{\mathit{\boldsymbol{0}}}$ corresponds to the operator of ${\mathit{\boldsymbol{f}}} = \mathit{\boldsymbol{0}}$.

Lemma 3.13. For operator $S_{\mathit{\boldsymbol{0}}}$, the following equation holds true

where $\mathit{\boldsymbol{\omega}}_h = \{D_0(\mathit{\boldsymbol{\omega}}_b;p_h),\mathit{\boldsymbol{\omega}}_b\},\ \mathit{\boldsymbol{q}}_h = \{D_0(\mathit{\boldsymbol{q}}_b;p_h),\mathit{\boldsymbol{q}}_b\}$, $D_{\mathit{\boldsymbol{0}}}$ and $L_{\mathit{\boldsymbol{0}}}$ correspond to the operator of ${\mathit{\boldsymbol{f}}} = \mathit{\boldsymbol{0}}$.

Proof. For any $\mathit{\boldsymbol{\omega}}_b, \mathit{\boldsymbol{q}}_b\in B_h^0$. From the definition of operator $S_{{\mathit{\boldsymbol{f}}}}$, we obtain

Let ${\mathit{\boldsymbol{f}}} = \mathit{\boldsymbol{0}}$ in (23), we have

We complete the proof.

Lemma 3.14. Assume that $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_h$ is the only solution of HWG algorithm (9), we have that ${\mathit{\boldsymbol{u}}}_h\in\mathcal V_h$ and ${\mathit{\boldsymbol{u}}}_b\in B_h$ are well defined functions and $S_{{\mathit{\boldsymbol{f}}}}({\mathit{\boldsymbol{u}}}_b;p_h) = \langle\!\langle\mathit{\boldsymbol{\zeta}}_h\rangle\!\rangle_e = \mathit{\boldsymbol{0}}$.

Proof. Because $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_h$ is the only solution of HWG algorithm (9). Then by Theorem 3.2, we obtain that for any $e \in \mathcal{E}_h^0$, $[\![{\mathit{\boldsymbol{u}}}_h]\!]_{e} = 0$, and ${\mathit{\boldsymbol{u}}}_b = Q_b\mathit{\boldsymbol{g}}$ on $\partial\Omega\,$, so ${\mathit{\boldsymbol{u}}}_h\in \mathcal V_h$ and ${\mathit{\boldsymbol{u}}}_b\in B_h$ are well defined. $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{0}}\}\in V_k(T)$ on $T$, and $\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{0}}$ elsewhere. Substituting (9), we obtain

where $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{0}},\mathit{\boldsymbol{v}}_b\}\in V_k(T)$ on the element $T$, and elsewhere $\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{0}}$. Substituting (9), then $\mathit{\boldsymbol{\lambda}}_h$ satisfies the following equation

where $\mathit{\boldsymbol{\lambda}}_{h,T}$ is the limit of $\mathit{\boldsymbol{\lambda}}_{h}$ on $\partial T$. From the definition of operator $S_{{\mathit{\boldsymbol{f}}}}$, we obtain

So $\langle\!\langle\mathit{\boldsymbol{\lambda}}_h\rangle\!\rangle_e = \mathit{\boldsymbol{0}}$, that is $S_{{\mathit{\boldsymbol{f}}}}({\mathit{\boldsymbol{u}}}_b;p_h) = \langle\!\langle \mathit{\boldsymbol{\zeta}}_h\rangle\!\rangle_e = \mathit{\boldsymbol{0}}$.

Lemma 3.15. Assume that $\bar{\mathit{\boldsymbol{u}}}_b\in B_h$ is a function satisfying $\bar{\mathit{\boldsymbol{u}}}_b = Q_b\mathit{\boldsymbol{g}}\,$ on $\partial\Omega\,$, and $\bar{\mathit{\boldsymbol{u}}}_b$ and $p_h$ satisfy the following operator equations:

Then $(\bar{\mathit{\boldsymbol{u}}}_h;p_h)\in V_h\times W_h$ is the solution of the WG algorithm (8), where $\bar{\mathit{\boldsymbol{u}}}_0$ is the solution of the following problem on each element $T\in \mathcal T_h$.

Proof. For each element $T\in \mathcal T_h$, $\bar{\mathit{\boldsymbol{\lambda}}}_{h,T}\in \Lambda_k(\partial T)$ can be solved from the following equation

Define $\bar{\mathit{\boldsymbol{\lambda}}}_h\in\Lambda_h$ as $\bar{\mathit{\boldsymbol{\lambda}}}_h|_{\partial T} = \bar{\mathit{\boldsymbol{\lambda}}}_{h,T}$. Because $(\bar{\mathit{\boldsymbol{u}}}_b;p_h)\in B_h\times W_h$ satisfies Lemma 3.15, $\bar{\mathit{\boldsymbol{u}}}_b$ satisfies the boundary condition and $\bar{\mathit{\boldsymbol{u}}}_0$ satisfies (25), then

Using (25) subtract (26), we have

Adding up all $T$ on $\mathcal{T}_{h}$ so that

Limiting $\mathit{\boldsymbol{v}}$ in weak function space $V_h^0$, and using (27), it is easy to get

Then

According to the assumption $\bar{\mathit{\boldsymbol{u}}}_b|_{\partial \Omega} = Q_b \mathit{\boldsymbol{g}}$ and Theorem 3.2, $\bar{\mathit{\boldsymbol{u}}}_h\in V_h$ is the solution of WG method (8).

From the above lemma, it is not difficult to prove the following theorem.

Theorem 3.16. Assume that $\bar{\mathit{\boldsymbol{u}}}_b \in B_h$ is a function satisfying that ${\mathit{\boldsymbol{u}}}_b = Q_b \mathit{\boldsymbol{g}}\,$ on $\partial\Omega\,$, $\bar{\mathit{\boldsymbol{u}}}_0$ is the solution to (25).Then $(\bar{\mathit{\boldsymbol{u}}}_h;p_h)\in V_h\times W_h$ is the solution of WG problem (8) if and only if $\bar{\mathit{\boldsymbol{u}}}_b$ satisfies the following operator equation

By (24) and (28), we have

Seeking the finite element $\mathit{\boldsymbol{G}}_h\in B_h$ satisfying: $\mathit{\boldsymbol{G}}_b = Q_b \mathit{\boldsymbol{g}}$ on $\partial\Omega\,$, and $\mathit{\boldsymbol{0}}$ elsewhere. Since $S_{\mathit{\boldsymbol{0}}}$ is a linear operator, we obtain

Substituting the above equation in (29) gives

We define $\mathit{\boldsymbol{H}}_b = \bar{\mathit{\boldsymbol{u}}}_b-\mathit{\boldsymbol{G}}_b$ such that $\mathit{\boldsymbol{H}}_b = \mathit{\boldsymbol{0}}$ on $\partial\Omega\,$. Let $\mathit{\boldsymbol{r}}_b = -S_{{\mathit{\boldsymbol{f}}}}(\mathit{\boldsymbol{0}};0)-S_{\mathit{\boldsymbol{0}}}(\mathit{\boldsymbol{G}}_b;p_h)$, then

Subtraction algorithm 1 The solution $({\mathit{\boldsymbol{u}}}_h;p_h)$ of the WG algorithm (8) can be obtained by the following steps

Step 1: On each element $T\in \mathcal T_h$, $\mathit{\boldsymbol{r}}_b$ can be solved by the following equation

Step 2: Solving $\{\mathit{\boldsymbol{H}}_b;p_h\}$ through (30).

Step 3: Calculating ${\mathit{\boldsymbol{u}}}_b = \mathit{\boldsymbol{G}}_b+\mathit{\boldsymbol{H}}_b$, we get the solution on the element boundary, and then on each element $T\in \mathcal T_h$, we calculate ${\mathit{\boldsymbol{u}}}_0 = D_{{\mathit{\boldsymbol{f}}}}({\mathit{\boldsymbol{u}}}_b,p_h)$ through (21).

4.   HWG for Brinkman equation with Neumann boundary condition (2)

In this section, we present HWG algorithm for Brinkman equation with Neumann boundary condition.

4.1.   Algorithm

First we present the WG numerical scheme of Brinkman first variational formulation.

Algorithm 4.1. ([17]) Find $(\bar{\mathit{\boldsymbol{u}}}_h;\bar{p}_h)\in V_h\times W_h$ such that $\nabla\bar{\mathit{\boldsymbol{u}}}_b \cdot {\mathit{\boldsymbol{n}}} = Q_b\mathit{\boldsymbol{\theta}}$ on $\partial\Omega$ and the following equation holds true:

for all $\mathit{\boldsymbol{v}} = \{{\mathit{\boldsymbol{v}}}_0,{\mathit{\boldsymbol{v}}}_b\}\in {V}_{h,N}$ and $q\in W_h$.

Similarly to the case of Dirichlet boundary condition, by introducing Lagrange multipliers, we introduce HWG method:

Algorithm 4.2. ([29]) Find $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_{h}$ such that $\nabla{{\mathit{\boldsymbol{u}}}_b}\cdot{\mathit{\boldsymbol{n}}} = Q_b\mathit{\boldsymbol{\theta}}$ on $\partial\Omega$ and satisfying the following equations:

for any $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}\in {V}_{h,N}$, $q\in W_h$ and ${\mathit{\boldsymbol{\mu}}}\in\Xi_{h}$.

Lemma 4.1. The problem (31) is well-posed.

Proof. The argument is similar to that for the Lemma 3.1. This completes the proof.

4.2.   Stability analysis

The proofs of the following lemmas are the same as Lemma 3.3 - 3.7.

Lemma 4.2. (Boundedness) There exists a constant $C>0$ such that

Lemma 4.3. (Positivity) For any $\mathit{\boldsymbol{v}}\in \mathcal V_h$, we have $\|\mathit{\boldsymbol{v}}\|^2_{V_{h,N}} = |\!|\!|\mathit{\boldsymbol{v}}|\!|\!|^2$, then

Lemma 4.4. (inf-sup condition 1) There exists a constant $\beta>0$ independent of $h$, for any $\rho\in W_h$, we have

Lemma 4.5. (inf-sup condition 1') For any $\rho\in W_h$, there exists a constant $\beta>0$ independent of $h$ and $\mathit{\boldsymbol{v}}\in V_{h,N}$, we have

Lemma 4.6. (inf-sup condition 2) For any $\mathit{\boldsymbol{\tau}}\in \Xi_{h}$, there exists $\mathit{\boldsymbol{v}}\in V_{h,N}$ satisfying $\mathit{\boldsymbol{v}}_0 = \mathit{\boldsymbol{0}}$, such that

where $C>0$ is a constant independent of $h$.

4.3.   Error equation

The purpose of this section is to construct the error equation between the numerical solution and the true solution for HWG according to the numerical solution algorithm (31).

Lemma 4.7. Assume that $({\mathit{\boldsymbol{u}}};p) \in[H_0^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution of (2), $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_h$ is the solution of (31). Let ${\mathit{\boldsymbol{\lambda}}} = \nabla {\mathit{\boldsymbol{u}}} \cdot {\mathit{\boldsymbol{n}}}-p {\mathit{\boldsymbol{n}}}$, take $\mathit{\boldsymbol{e}}_{h} = \{Q_0 {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_0, Q_b {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_b\},\ \varepsilon_{h} = \mathbb{Q}_h p-p_{h},\ {\mathit{\boldsymbol{\delta}}}_{h} = Q_{b}\mathit{\boldsymbol{\lambda}}-{\mathit{\boldsymbol{\lambda}}}_{h}$. Error function $\mathit{\boldsymbol{e}}_{h}$, $\varepsilon_{h}$, and ${\mathit{\boldsymbol{\delta}}}_{h}\,$ satisfy the following equation

where

Proof. By Lemma 3.9, we obtain

By combining with (31), we obtain

Now, from Theorem 3.2 and Lemma 3.9, we have

We complete the proof.

4.4.   Error estimation

In this section, we first give the following lemmas to make the corresponding $H^1$ and $L^2$ error estimates for the error equation.

Lemma 4.8. Assume that $({\mathit{\boldsymbol{u}}};p)\in [H^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution to the problem (2). Then, there exists a constant $C$ such that

Proof. By Lemma 3.10, we obtain

This completes the proof. The proofs of the following Theorems are similar to those of Theorem 3.11 and Theorem 3.12. Therefore, we only state the the conclusion without proofs.

Theorem 4.9. Assume that $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution to the problem (2), and $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_{h,N}\times W_h \times \Xi_h$ is the solution of (31). We have

Finally, the dual technique is used to derive the optimal order error estimates of the weak finite element scheme under $L^2$ norm. We consider the following dual problems

where $(\psi;\xi)\in [H^2(\Omega)]^d\times H^1(\Omega)$. Assume that the above dual problem has $H^2$ - regularity, that is, there is a constant $C$, which makes

Theorem 4.10. Assume that $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution to the problem (2), and $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_{h,N}\times W_h\times \Xi_h$ is the solution of (31). We have

5.   HWG for Brinkman problem with Robin boundary condition (3)

In this section, we present HWG for Brinkman equation with Robin boundary condition.

5.1.   Algorithm

We, first give a weak Galerkin finite element numerical scheme of Brinkman first variational formulation under Robin boundary condition as follows:

Algorithm 5.1. ([17]) Find $(\bar{\mathit{\boldsymbol{u}}}_h;\bar{p}_h)\in V_h\times W_h$ such that $\nabla\bar{\mathit{\boldsymbol{u}}}_b\cdot\mathit{\boldsymbol{n}}+\alpha\bar {\mathit{\boldsymbol{u}}}_b = \gamma$ on $\partial\Omega$ and the following equations hold true

for any $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}\in \mathcal V_h$ and $q\in W_h$.

Similar to the other two boundary cases, by introducing Lagrange multipliers, we introduce the HWG method for Robin boundary case as follows:

Algorithm 5.2. ([29]) Find $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_{h}$ such that $\nabla{\mathit{\boldsymbol{u}}}_b\cdot\mathit{\boldsymbol{n}}+\alpha\mathit{\boldsymbol{u}}_b = \gamma$ on $\partial\Omega$ and the following equations hold true:

for any $\mathit{\boldsymbol{v}} = \{\mathit{\boldsymbol{v}}_0,\mathit{\boldsymbol{v}}_b\}\in V_h$, $q\in W_h\,$ and ${\mathit{\boldsymbol{\mu}}}\in\Xi_{h}$.

We can easily establish the following well-posedness of the problem (32).

Lemma 5.1. The problem (32) is well-posed.

5.2.   Stability analysis

Lemma 5.2. (Boundedness) There exists a constant $C>0$ such that

Lemma 5.3. (Positivity) For any $\mathit{\boldsymbol{v}}\in \mathcal V_h$, we have $\|\mathit{\boldsymbol{v}}\|^2_{V_{h}} = |\!|\!|\mathit{\boldsymbol{v}}|\!|\!|^2$, then

Lemma 5.4. (inf-sup condition 1) There exists a constant $\beta>0$ independent of $h$, for any $\rho\in W_h$, we have

Lemma 5.5. (inf-sup condition 1') For any $\rho\in W_h$, there exists a constant $\beta>0$ and $\mathit{\boldsymbol{v}}\in V_{h}$ independent of $h$, we have

Lemma 5.6. (inf-sup condition 2) For any $\mathit{\boldsymbol{\tau}}\in \Xi_{h}$, there exists $\mathit{\boldsymbol{v}}\in V_{h,N}$ satisfying $\mathit{\boldsymbol{v}}_0 = \mathit{\boldsymbol{0}}$, such that

where $C>0$ is a constant independent of $h$.

5.3.   Error equation

The purpose of this section is to construct the error equations between the numerical solution and the true solution for HWG, according to the numerical solution algorithm (32).

Lemma 5.7. Assume that $({\mathit{\boldsymbol{u}}};p) \in[H_0^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution to the problem (3), and $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h\times\Xi_h$ is the solution of (32). Let ${\mathit{\boldsymbol{\lambda}}} = \nabla {\mathit{\boldsymbol{u}}} \cdot {\mathit{\boldsymbol{n}}}-p {\mathit{\boldsymbol{n}}}$, take $\mathit{\boldsymbol{e}}_{h} = \{Q_0 {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_0, Q_b {\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_b\},\ \varepsilon_{h} = \mathbb{Q}_h p-p_{h},\ {\mathit{\boldsymbol{\delta}}}_{h} = Q_{b}\mathit{\boldsymbol{\lambda}}-{\mathit{\boldsymbol{\lambda}}}_{h}$. Error function $\mathit{\boldsymbol{e}}_{h}$, $\varepsilon_{h}$, and ${\mathit{\boldsymbol{\delta}}}_{h}\,$ satisfy the following equation $\varepsilon_{h}$ satisfy the following equation

where

Proof. Similar to the Lemma 3.9, we have that

By combining it with (32), we obtain that

which gives

From Theorem 3.2 and Lemma 3.9, we obtain

This completes the proof.

5.4.   Error estimation

In this section, we want to give the $H^1$ and $L^2$ error estimates, so the following lemmas are given first.

Lemma 5.8. Assume that $({\mathit{\boldsymbol{u}}};p) \in[H_0^1(\Omega)]^d\times L_0^2(\Omega)$ is the true solution to the problem (3). There exists constant $C$ satisfying

Proof. From Lemma 3.10, we have

From Cauchy-Schwarz inequality and trace inequality, we obtain

we complete the proof.

The proofs of the following Theorems are similar to that of Theorem 3.11.

Theorem 5.9. Assume that $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution which satisfies the problem (3), and $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_h)\in V_h\times W_h \times \Xi_h$ is the solution of (32). We have

Finally, the dual technique is used to derive the optimal order error estimates of the weak Galerkin finite element scheme under $L^2$ norm. Consider the following dual problem:

with $(\psi;\xi)\in [H^2(\Omega)]^d\times H^1(\Omega)$. Assume that the above dual problem has $H^2$ - regularity, that is, there is a constant $C$, which makes

Theorem 5.10. Assume that $({\mathit{\boldsymbol{u}}};p) \in \{[H^{k+1}(\Omega)]^d\cap [H^1_0(\Omega)]^d\}\times L_0^2(\Omega)$ is the true solution to the problem (3), $({\mathit{\boldsymbol{u}}}_h;p_h;\mathit{\boldsymbol{\lambda}}_{h})\in V_{h}\times W_h\times \Xi_{h}$ is the solution of (32). When $k\geq 2$, we have

Proof. Using $\mathit{\boldsymbol{e}}_0\,$ to act on both ends of (33), we obtain

Take ${\mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{\psi}},\,\mathit{\boldsymbol{v}}_0 = \mathit{\boldsymbol{e}}_0,\,p = \xi\,$ in the above formula. From the error equation, we obtain

Using Theorem 3.12, we have

Since we have

we obtain

This completes the proof.

6.   Numerical experiments

In this section, we consider Brinkman problem (1) on the partition region $\Omega = (0,1)^2$, where we consider different $\mu$ and $\kappa$ given as follows:

where $a$ is a constant and a number of different values of $a$ have been tested. Our results show that the proposed method produce robust numerical solutions for varying parameters $\mu$ and $a$.

We shall take the following analytical solution:

According to (1), we can get the exact $\mathit{\boldsymbol{f}}$ and let $h$ denote the grid size. For simplicity, we choose the polynomial degree $k = 1$. We shall set $\mathit{\boldsymbol{e}}_h = Q_h{\mathit{\boldsymbol{u}}}-{\mathit{\boldsymbol{u}}}_h$, $\varepsilon_h = \mathbb Q_h p-p_h$, and ${\mathit{\boldsymbol{\delta}}}_h = Q_b \mathit{\boldsymbol{\lambda}}-\mathit{\boldsymbol{\lambda}}_h$.