#  Analysis of weak Galerkin mixed FEM based on the velocity–pseudostress
formulation for Navier–Stokes equation on polygonal meshes

Zeinab Gharibi Mehdi Dehghan CI2MA and GIMNAP-Departamento de Matemática,
Universidad del Bío-Bío, Casilla 5-C, Concepción, Chile, email:
<EMAIL_ADDRESS>of Applied Mathematics, Faculty of Mathematics
and Computer Sciences, Amirkabir University of Technology (Tehran
Polytechnic), No. 424, Hafez Ave., Tehran, 15914, Iran, email:
<EMAIL_ADDRESS>
###### Abstract

The present article introduces, mathematically analyzes, and numerically
validates a new weak Galerkin (WG) mixed-FEM based on Banach spaces for the
stationary Navier–Stokes equation in pseudostress-velocity formulation. More
precisely, a modified pseudostress tensor, called ${\boldsymbol{\sigma}}$,
depending on the pressure, and the diffusive and convective terms has been
introduced in the proposed technique, and a dual-mixed variational formulation
has been derived where the aforementioned pseudostress tensor and the
velocity, are the main unknowns of the system, whereas the pressure is
computed via a post-processing formula. Thus, it is sufficient to provide a WG
space for the tensor variable and a space of piecewise polynomial vectors of
total degree at most ’k’ for the velocity. Moreover, in order to define the
weak discrete bilinear form, whose continuous version involves the classical
divergence operator, the weak divergence operator as a well-known alternative
for the classical divergence operator in a suitable discrete subspace is
proposed. The well-posedness of the numerical solution is proven using a
fixed-point approach and the discrete versions of the Babuška–Brezzi theory
and the Banach–Nečas–Babuška theorem. Additionally, an a priori error estimate
is derived for the proposed method. Finally, several numerical results
illustrating the method’s good performance and confirming the theoretical
rates of convergence are presented.

Keywords: Weak Galerkin. Pseudostress-velocity formulation. Mixed finite
element methods. Navier–Stokes equation. Well-posedness. Error analysis.
Mathematics Subject Classification: 34B15, 65L60, 65M70.

##  1 Introduction

###  1.1 Scope

The weak Galerkin method (WG), introduced in [31] for second-order elliptic
equations, approximates the differential operators in the variational
formulation through a framework that emulates the theory of distributions for
piecewise polynomials. Additionally, the usual regularity of the approximating
functions in this method is compensated by carefully-designed stabilizers. In
turn, the WG methods have been studied for solving various models (cf. (cf.
[28, 36, 36, 29, 31, 32, 11, 25, 27, 12, 21, 13])), demonstrating their
substantial potential as a formidable numerical method in scientific
computing. The key distinction between WG methods and other existing finite
element techniques lies in their utilization of weak derivatives and weak
continuities in formulating numerical schemes based on conventional weak forms
for the underlying PDE problems. Because of their significant structural
pliability, WG methods are aptly suited for a broad class of PDEs, ensuring
the requisite stability and accuracy in approximations.

On the other hand, the development of certain mixed finite element methods has
become a new research field aimed at solving linear and nonlinear problems
using velocity-pressure, stress-velocity-pressure, and pseudostress-velocity-
based formulations. A principal benefit of stress-based and pseudostress-based
formulations is the immediate calculation of physical quantities, such as
stress, without the need for velocity derivatives, preventing the degradation
of precision. Despite the extensive research on the weak Galerkin FEM, there
are few publications on the weak Galerkin mixed-FEM in the literature.
Specifically, and up to our knowledge, [32] and [21] are the only two works in
which the weak Galerkin mixed-FEM for linear and nonlinear PDEs based on the
dual formulation are proposed and analyzed, including corresponding optimal
errors estimates. Although many research works have been studied and analyzed
on the weak Galerkin method for the Navier–Stokes problem (cf. [25, 35, 33,
34, 24]), most existing schemes are extracted by the Stokes-type (or primal-
mixed) formulation, where the velocity and pressure are principle unknowns.

According to the above discussion, the objective of the present paper is to
develop and design a new mixed FEM based on the nonstandard Banach spaces
drawn by Gatica et al. [18] within the framework of the weak Galerkin method
to solve the stationary Navier–Stokes equation. More precisely, we are
particularly interested in the development of mixed formulations not involving
any augmentation procedure (as done, e.g. in [7, 6]). To this end, we now
extend the applicability of the approach employed in [8] for the fluid model,
in the framework weak Galerkin method. In fact, instead of using the primal or
dual-mixed method (see e.g. [32, 21]), we now employ a modified mixed
formulation and adapt the approach from [22] to propose, up to our knowledge
by the first time, a weak Galerkin mixed-FEM for Navier–Stokes, which consists
of introducing the gradient of velocity and a vector version of the Bernoulli
tensor as further unknowns. In this way, and besides eliminating the pressure,
which can be approximated later on via postprocessing, the resulting mixed
variational formulation does not need to incorporate any augmented term, and
it yields basically the same Banach saddle-point structure.

###  1.2 Outline and notations

The remainder of the paper is organized as follows. In Sec. 2, we refer to
Ref. [8] to introduce the model of interest, recapitulate the associated dual-
mixed variational formulation with the unknowns ${\boldsymbol{\sigma}}$ and
${\mathbf{u}}$ living in suitable Banach spaces, and state the main result
establishing its well-posedness. We then introduce the WG discretization in
Sec. 3, following Refs. [32] and [21]. This section comprises four main parts:
first, we state basic assumptions on the mesh; second, we define the local WG
space, projections, and weak differential operators; third, we discuss their
approximation properties; and fourth, we derive the global WG subspace and the
disceret scheme. In Sec. 4, we analyze the solvability of our discrete scheme
using a fixed-point strategy. To accomplish this, we derive common estimates
for the bilinear and trilinear forms, as well as the discrete inf-sup
condition. The classical Banach fixed-point theorem then allows us to obtain
the principal result. Sec. 5 is dedicated to deriving a priori error estimates
for the numerical solution under a small data assumption. Finally, in Sec. 6,
we present some numerical experiments to confirm the theoretical correctness
and effectiveness of the discrete schemes.

For any vector fields $\mathbf{v}=(v_{1},v_{2})^{\tt t}$ and
$\mathbf{w}=(w_{1},w_{2})^{\tt t}$, we set the gradient, divergence and tensor
product operators as

$\boldsymbol{\nabla}\mathbf{v}:=(\nabla v_{1},\nabla
v_{2}),\quad\operatorname{div}(\mathbf{v}):=\partial_{x}v_{1}+\partial_{y}v_{2},{\quad\hbox{and}\quad}\mathbf{v}\otimes\mathbf{w}\,:=\,(v_{i}\,w_{j})_{i,j=1,2}\,,$

respectively. In addition, denoting by $\mathbb{I}$ the identity matrix of
$\mathrm{R}^{2}$, for any tensor fields
${\boldsymbol{\tau}}=(\tau_{ij}),~{}{\boldsymbol{\zeta}}=(\zeta_{ij})\in\mathrm{R}^{2\times
2}$, we write as usual

${\boldsymbol{\tau}}^{\mathtt{t}}:=(\tau_{ji}),\quad\operatorname{tr}({\boldsymbol{\tau}}):=\tau_{11}+\tau_{22},\quad{\boldsymbol{\tau}}^{\mathtt{d}}:={\boldsymbol{\tau}}-\dfrac{1}{2}\operatorname{tr}({\boldsymbol{\tau}})\mathbb{I}\,,{\quad\hbox{and}\quad}{\boldsymbol{\tau}}:{\boldsymbol{\zeta}}:=\sum_{i,j=1}^{2}\tau_{ij}\zeta_{ij}\,,$

which corresponds, respectively, to the transpose, the trace, and the deviator
tensor of ${\boldsymbol{\tau}}$, and to the tensorial product between
${\boldsymbol{\tau}}$ and ${\boldsymbol{\zeta}}$.

Throughout the paper, given a bounded domain $\Omega$, we let $\mathcal{O}$ be
any given open subset of $\Omega$. By $(\cdot,\cdot)_{0,\mathcal{O}}$ and
$\|\cdot\|_{0,\mathcal{O}}$ we denote the usual integral inner product and the
corresponding norm of $\mathrm{L}^{2}(\mathcal{O})$, respectively. For
positive integers $m$ and $r$, we shall use the common notation for the
Sobolev spaces $\mathrm{W}^{m,r}(\mathcal{O})$ with the corresponding norm and
semi-norm $\|\cdot\|_{m,r,\mathcal{O}}$ and $|\cdot|_{m,r,\mathcal{O}}$,
respectively; and if $r=2$, we set
$\mathrm{H}^{m}(\mathcal{O}):=\mathrm{W}^{m,2}(\mathcal{O})$,
$\|\cdot\|_{m,\mathcal{O}}:=\|\cdot\|_{m,2,\mathcal{O}}$ and
$|\cdot|_{m,\mathcal{O}}:=|\cdot|_{m,2,\mathcal{O}}$. If $\mathcal{O}=\Omega$,
the subscript will be omitted. Furthermore, M and $\mathbb{M}$ represent
corresponding vectorial and tensorial counterparts of the scalar functional
space $\mathrm{M}$. On the other hand, given $t\in(1,+\infty)$, letting
$\operatorname{\mathbf{div}}$ be the usual divergence operator
$\operatorname{div}$ acting along the rows of a given tensor, we introduce the
standard Banach space

$\mathbb{H}({\mathbf{div}}_{t};\Omega)\,:=\,\Big{\\{}{\boldsymbol{\tau}}\in\mathbb{L}^{2}(\Omega):\quad{\mathbf{div}}({\boldsymbol{\tau}})\in\mathbf{L}^{t}(\Omega)\Big{\\}}\,,$

equipped with the usual norm

$\|{\boldsymbol{\tau}}\|_{{\mathbf{div}}_{t};\Omega}\,:=\,\|{\boldsymbol{\tau}}\|_{0,\Omega}+\|{\mathbf{div}}({\boldsymbol{\tau}})\|_{0,t;\Omega},\quad\quad\forall\,{\boldsymbol{\tau}}\in\mathbb{H}({\mathbf{div}}_{t};\Omega),$

##  2 The model problem and its continuous formulation

Consider a spatial bounded domain $\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) with
a Lipschitz continuous boundary $\partial\Omega$ with outward-pointing unit
normal $\mathbf{n}$. We focus on solving the Navier–Stokes equation with
viscosity $\nu$, where, given the body force term
$\textbf{f}\in\textbf{L}^{2}(\Omega)$ and suitable boundary data
$\mathbf{g}\in\textbf{H}^{1/2}(\partial\Omega)$, the objective is to find a
velocity field $\mathbf{u}$ and a pressure field $p$ such that

$\begin{array}[]{rcll}-\nu\,\boldsymbol{\Delta}{\mathbf{u}}+{\mathbf{u}}\cdot\boldsymbol{\nabla}{\mathbf{u}}+\nabla
p&=&\,\textbf{f}&{\quad\hbox{in}\quad}\Omega\,,\\\\[8.61108pt]
{\mathrm{div}}({\mathbf{u}})&=&\,0&{\quad\hbox{in}\quad}\Omega\,,\\\\[8.61108pt]
{\mathbf{u}}&=&\,\mathbf{g}&{\quad\hbox{on}\quad}\partial\Omega\,,\end{array}$
(2.1)

In addition, in order to guarantee uniqueness of the pressure, this unknown
will be sought in the space

$L_{0}^{2}(\Omega)\,:=\,\Big{\\{}q\in
L^{2}(\Omega):\quad\int_{\Omega}p\,=\,0\Big{\\}}\,.$ (2.2)

Note that, due to the incompressibility of the fluid (cf. second row of
(2.1)), $\mathbf{g}$ must satisfy

$\displaystyle\int_{\Omega}\mathbf{g}\cdot{\mathbf{n}}\,=\,0\,.$

For the subsequent analysis we assume that coefficient $\nu$ is piecewise
constant and positive.

Next, to obtain a velocity–pseudostress formulation, the first step is to
rewrite equation (2.1) so that stress and velocity are only unknowns of
equation. To achieve this, we introduce a tensor field denoted by
${\boldsymbol{\sigma}}$, represented as

$\boldsymbol{\sigma}\,:=\,\nu\boldsymbol{\nabla}\mathbf{u}-\mathbf{u}\otimes\mathbf{u}-(p+r_{\mathbf{u}})\mathbb{I}\quad{\quad\hbox{in}\quad}\Omega\,,$
(2.3)

where

$r_{\mathbf{u}}\,:=\,-c_{r}(\operatorname{tr}(\mathbf{u}\otimes\mathbf{u}),1)_{0,\Omega}\,=\,-c_{r}(\mathbf{u},\mathbf{u})_{0,\Omega}\qquad\text{with}~{}~{}c_{r}\,=\,\dfrac{1}{2|\Omega|}\,.$

In this way, applying the trace operator to both sides of (2.3), and utilizing
the incompressibility condition $\operatorname{div}(\mathbf{u})=0$, one
arrives at

$p\,=\,-\dfrac{1}{2}\left(\operatorname{tr}(\boldsymbol{\sigma})+\operatorname{tr}(\mathbf{u}\otimes\mathbf{u})\right)-r_{\mathbf{u}}\quad{\quad\hbox{in}\quad}\Omega\,.$
(2.4)

which allows us to eliminate the pressure variable from the formulation. In
turn, according to (2.4), the assumption (2.2) becomes

$\displaystyle\int_{\Omega}{\mathrm{tr}}({\boldsymbol{\sigma}})\,=\,0\,.$
(2.5)

Hence, after replacing back (2.3) in (2.1), gathering the resulting equation
and (2.5), we have the following problem, which contains unknowns
${\boldsymbol{\sigma}}$ and ${\mathbf{u}}$.

######  Problem 1 (Model problem).

Find $\boldsymbol{\sigma},\mathbf{u}$ such that

$\begin{cases}\hskip
2.84544pt\boldsymbol{\sigma}^{\mathtt{d}}+(\mathbf{u}\otimes\mathbf{u})^{\mathtt{d}}\,=\,\nu\,\boldsymbol{\nabla}\mathbf{u}&{\quad\hbox{in}\quad}\Omega\,,\\\\[2.84526pt]
\hskip
36.98866pt\operatorname{\mathbf{div}}(\boldsymbol{\sigma})\,=\,-\mathbf{f}&{\quad\hbox{in}\quad}\Omega\,,\\\\[2.84526pt]
\hskip
31.2982pt\displaystyle\int_{\Omega}{\mathrm{tr}}({\boldsymbol{\sigma}})\,=\,0\,.\end{cases}$

supplied with the following boundary condition

$\mathbf{u}\,=\,\mathbf{g}\quad{\quad\hbox{on}\quad}\partial\Omega\,.$

Next, in order to derive a velocity–pseudostress based-mixed formulation for
Problem 1, we let $\mathbb{X}$ and $\mathbf{Y}$ be corresponding test spaces,
and then proceed to multiply the first and second equations of Problem 1 by
${\boldsymbol{\tau}}$ and ${\mathbf{v}}$, respectively, and use the fact that
${\mathrm{tr}}({\boldsymbol{\tau}}^{\mathtt{d}})=0$, to get

$\dfrac{1}{\nu}\int_{\Omega}{\boldsymbol{\sigma}}^{\mathtt{d}}:{\boldsymbol{\tau}}^{\mathtt{d}}+\dfrac{1}{\nu}\int_{\Omega}({\mathbf{u}}\otimes{\mathbf{u}})^{\mathtt{d}}:{\boldsymbol{\tau}}\,=\,\int_{\Omega}\boldsymbol{\nabla}{\mathbf{u}}:{\boldsymbol{\tau}}\qquad\forall\,{\boldsymbol{\tau}}\in\mathbb{X}\,,$
(2.6)

and

$\int_{\Omega}{\mathbf{div}}({\boldsymbol{\sigma}})\cdot{\mathbf{v}}\,=\,-\int_{\Omega}\textbf{f}\cdot{\mathbf{v}}\qquad\forall\,{\mathbf{v}}\in\mathbf{Y}\,,$
(2.7)

it is easy to notice that, thanks to Cauchy–Schwarz’s inequality, the first
term on the left hand side of (2.6) makes sense for
${\boldsymbol{\sigma}},{\boldsymbol{\tau}}\in\mathbb{L}^{2}(\Omega)$. In turn,
regarding the term on the right hand side of (2.6), assuming originally that
${\mathbf{u}}\in\mathbf{H}^{1}(\Omega)$, and given
$t,t^{\prime}\in(1,\infty)$, conjugate to each other, we can integrate by
parts with ${\boldsymbol{\tau}}\in\mathbb{H}({\mathbf{div}}_{t};\Omega)$, so
that using the Dirichlet boundary conditions provided in Problem 1, we obtain

$\int_{\Omega}\boldsymbol{\nabla}{\mathbf{u}}:{\boldsymbol{\tau}}\,=\,-\int_{\Omega}{\mathbf{u}}\cdot{\mathbf{div}}({\boldsymbol{\tau}})+\langle{\boldsymbol{\tau}}{\mathbf{n}},\mathbf{g}\rangle\quad\forall\,{\boldsymbol{\tau}}\in\mathbb{H}({\mathbf{div}}_{t};\Omega)\,,$
(2.8)

where $\langle\cdot,\cdot\rangle_{\Gamma}$ stands for the duality
$(\mathbf{H}^{-1/2}(\Gamma),\mathbf{H}^{1/2}(\Gamma)\big{)}$. Now, from the
first term on the right hand side of the foregoing equation, along with the
Sobolev embedding
$\mathbf{H}^{1}(\Omega)\subset\mathbf{L}^{t^{\prime}}(\Omega)$ , we realize
that it actually suffices to look for
${\mathbf{u}}\in\mathbf{L}^{t^{\prime}}(\Omega)$. However, it is clear from
(2.6) that its second term is well defined if
${\mathbf{u}}\in\mathbf{L}^{4}(\Omega)$, which yields $t^{\prime}=4$ and thus
$t=4/3$.

At this point, in order to deal with the null mean value of
${\mathrm{tr}}({\boldsymbol{\sigma}})$ (cf. third row of Problem 1), we
introduce the subspace of $\mathbb{H}({\mathbf{div}}_{4/3};\Omega)$ given by

$\mathbb{H}_{0}(\operatorname{\mathbf{div}}_{4/3};\Omega)\,:=\,\Big{\\{}{\boldsymbol{\tau}}\in\mathbb{X}:\quad\int_{\Omega}{\mathrm{tr}}({\boldsymbol{\tau}})\,=\,0\Big{\\}}\,.$

Then, testing the new (2.6) against
${\boldsymbol{\tau}}\in\mathbb{H}(\operatorname{\mathbf{div}}_{4/3};\Omega)$
is equivalent to doing it against
${\boldsymbol{\tau}}\in\mathbb{H}_{0}(\operatorname{\mathbf{div}}_{4/3};\Omega)$,
and taking into account the above discussion, we define the testing spaces as

$\mathbb{X}\,:=\,\mathbb{H}_{0}(\operatorname{\textbf{div}}_{4/3};\Omega),\quad\quad\text{with}\quad\quad\|\cdot\|_{\mathbb{X}}\,:=\,\|\cdot\|_{\operatorname{\textbf{div}}_{4/3};\Omega}$

and

$\mathbf{Y}\,:=\,\textbf{L}^{4}(\Omega),\qquad\qquad\text{with}\quad\quad\|\cdot\|_{\mathbf{Y}}\,:=\,\|\cdot\|_{0,4,\Omega}.$

Let us introduce the following bilinear (and trilinear) forms

$\displaystyle\mathcal{A}(\cdot,\cdot):\mathbb{X}\times\mathbb{X}\rightarrow\mathrm{R}$
$\displaystyle\qquad\qquad\mathcal{A}({\boldsymbol{\zeta}},{\boldsymbol{\tau}})\,:=\,\dfrac{1}{\nu}\int_{\Omega}{\boldsymbol{\zeta}}^{\mathtt{d}}:{\boldsymbol{\tau}}^{\mathtt{d}}\,,$
$\displaystyle\mathcal{B}(\cdot,\cdot):\mathbb{X}\times\mathbf{Y}\rightarrow\mathrm{R}$
$\displaystyle\qquad\qquad\mathcal{B}({\boldsymbol{\tau}},{\mathbf{v}})\,:=\,\int_{\Omega}{\mathbf{v}}\cdot\operatorname{\textbf{div}}({\boldsymbol{\tau}})\,,$
$\displaystyle\mathcal{C}(\cdot;\cdot,\cdot):\mathbf{Y}\times\mathbf{Y}\times\mathbb{X}\rightarrow\mathrm{R}$
$\displaystyle\qquad\qquad\mathcal{C}(\textbf{w},{\mathbf{v}};{\boldsymbol{\tau}})\,:=\,\dfrac{1}{\nu}\int_{\Omega}(\textbf{w}\otimes{\mathbf{v}})^{\mathtt{d}}:{\boldsymbol{\tau}}\,,$

and the linear functionals (associated to given data)

$\mathcal{G}({\boldsymbol{\tau}})\,:=\,\int_{\partial\Omega}\textbf{g}\cdot{\boldsymbol{\tau}}\mathbf{n}\quad\forall\,{\boldsymbol{\tau}}\in\mathbb{X}{\quad\hbox{and}\quad}\mathcal{F}({\mathbf{v}})\,:=\,-\int_{\Omega}\textbf{f}\cdot{\mathbf{v}}\quad\forall\,{\mathbf{v}}\in\mathbf{Y}\,.$

Then, with these forms at hand, the variational formulation of Problem 1 reads
as follows:

######  Problem 2 (variational problem).

Find the tensor $\boldsymbol{\sigma}\in\mathbb{X}$ and the velocity
$\mathbf{u}\in\mathbf{Y}$ such that

$\left\\{\begin{array}[]{ll}\mathcal{A}(\boldsymbol{\sigma},{\boldsymbol{\tau}})+\mathcal{C}(\mathbf{u},\mathbf{u};{\boldsymbol{\tau}})+\mathcal{B}({\boldsymbol{\tau}},\mathbf{u})\,=\,\mathcal{G}({\boldsymbol{\tau}})\quad\forall\,{\boldsymbol{\tau}}\in\mathbb{X}\,,\\\\[8.53581pt]
\hskip
119.50148pt\mathcal{B}(\boldsymbol{\sigma},\mathbf{v})\,=\,\mathcal{F}(\mathbf{v})\quad\forall\,{\mathbf{v}}\in\mathbf{Y}\,.\end{array}\right.$

The solvability result concerning Problem (2) is established as follows.

######  Theorem 2.1.

Let $\delta>0$ be constant related to the inf-sup condition of the linear part
of the left-hand side of Problem 2 (cf. Ref. [8, Eq. (3.29)]) and $c_{g}$ be
the upper bound of $\mathcal{G}(\cdot)$, define the ball

$\widehat{\mathbf{Y}}\,:=\,\Big{\\{}\mathbf{z}\in\mathbf{Y}:\quad\quad\|\mathbf{z}\|_{\mathbf{Y}}\,\leq\,\dfrac{\delta\nu}{2}\Big{\\}}\,,$

and assume that the given data satisfy

$\big{(}\dfrac{\nu\delta}{2}\big{)}^{-2}\Big{(}c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}+\|\mathbf{f}\|_{0,4/3}\Big{)}\,\leq\,\dfrac{1}{\nu}\,.$

Then, there exists a unique solution
$({\boldsymbol{\sigma}},\mathbf{u})\in\mathbb{X}\times\widehat{\mathbf{Y}}$
for Problem 2, and there holds the following stability estimate

$\|{\boldsymbol{\sigma}}\|_{\mathbb{X}}+\|\mathbf{u}\|_{\mathbf{Y}}\,\leq\,\dfrac{2}{\delta}\Big{(}c_{g}\|\mathbf{g}\|_{1/2,\partial\Omega}+\|\mathbf{f}\|_{0,4/3}\Big{)}\,.$

###### Proof.

See [8, proof of Theorem 3.8]. ∎

##  3 Weak Galerkin approximation

This section aims to introduce the weak Galerkin spaces and the discrete
bilinear form essential for introducing a weak Galerkin mixed FEM scheme. We
focus on the construction of method in the 2D case for simplicity.

###  3.1 Various tools in weak Galerkin method

A fundamental aspect of the weak Galerkin method lies in employing uniquely
defined weak derivatives instead of traditional derivative operators. Our
emphasis here centers on the weak divergence operator, a crucial step in
introducing our Weak Galerkin technique. To facilitate this discussion, we
begin by offering an overview of the mesh structure.
Mesh notation. Let $\mathcal{K}_{h}=\\{K\\}$ be a shape regular mesh of domain
$\Omega$ that consists of arbitrary polygon elements, where the mesh size
$h=\max\\{h_{K}\\}$, $h_{K}$ is the diameter of element $K$. The interior and
the boundary of any element $K\in\mathcal{K}_{h}$, are represented by $K^{0}$
and $\partial K$, respectively. Denote by $\mathcal{E}_{h}$ the set of all
edges in $\mathcal{K}_{h}$, and let
$\mathcal{E}_{h}^{0}=\mathcal{E}_{h}\backslash\partial\Omega$ be the set of
all interior edges. Here is a set of normal directions on $\mathcal{E}_{h}$:

$\mathrm{D}_{h}\,:=\,\bigg{\\{}\mathbf{n}_{e}:~{}~{}~{}~{}\mathbf{n}_{e}~{}\text{is
unit and outward normal to}~{}e,~{}~{}~{}~{}\text{for
all}~{}e\in\mathcal{E}_{h}\bigg{\\}}\,.$ (3.1)

Weak divergence operator and weak Galerkin space. It is well known that the
weak divergence operator is well-defined for weak matrix-valued functions
${\boldsymbol{\tau}}=\\{{\boldsymbol{\tau}}_{0},{\boldsymbol{\tau}}_{b}\\}$ on
the element $K$ such that ${\boldsymbol{\tau}}_{0}\in\mathbb{L}^{2}(K)$ and
${\boldsymbol{\tau}}_{b}\mathbf{n}_{e}\in\mathbf{H}^{-1/2}(\partial K)$, where
$\mathbf{n}_{e}\in\mathrm{D}_{h}|_{K}$. Components ${\boldsymbol{\tau}}_{0}$
and ${\boldsymbol{\tau}}_{b}$ can be understood as the value of function
${\boldsymbol{\tau}}$ in $K^{0}$ and on $\partial K$, respectively. We follow
[32], and introduce for each $K\in\mathcal{K}_{h}$ the local weak tensor space

$\mathbb{X}_{\mathtt{w}}(K)\,:=\,\Big{\\{}{\boldsymbol{\tau}}=\\{{\boldsymbol{\tau}}_{0},{\boldsymbol{\tau}}_{b}\\}:~{}~{}~{}{\boldsymbol{\tau}}_{0}\in\mathbb{L}^{2}(K){\quad\hbox{and}\quad}{\boldsymbol{\tau}}_{b}\mathbf{n}_{e}\in\mathbf{H}^{-1/2}(\partial
K)\quad\forall\,\mathbf{n}_{e}\in\mathrm{D}_{h}|_{K}\Big{\\}}\,.$ (3.2)

The global space $\mathbb{X}_{\mathtt{w}}$ is defined by gluing together all
local spaces $\mathbb{X}_{\mathtt{w}}(K)$ for any $K\in\mathcal{K}_{h}$. Now,
we define the weak divergence operator for matrix-valued functions as follows.

######  Definition 3.1 ([31]).

For any weak matrix-valued function
${\boldsymbol{\tau}}\in\mathbb{X}_{\mathtt{w}}(K)$ and element
$K\in\mathcal{K}_{h}$, the weak divergence operator, denoted by
$\operatorname{\mathbf{div}}_{\mathtt{w}}$, is defined as the unique vector-
valued function
$\operatorname{\mathbf{div}}_{\mathtt{w}}({\boldsymbol{\tau}})\in\mathbf{H}^{1}(K)$
satisfying

$({\mathbf{div}}_{\mathtt{w}}({\boldsymbol{\tau}}),\,{\boldsymbol{\zeta}})_{0,K}\,:=\,-({\boldsymbol{\tau}}_{0},\,\nabla{\boldsymbol{\zeta}})_{0,K^{0}}\,+\,\langle{\boldsymbol{\tau}}_{b}\mathbf{n}_{e},{\boldsymbol{\zeta}}\rangle_{0,\partial
K}\quad\forall\,{\boldsymbol{\zeta}}\in\mathbf{H}^{1}(K)\,.$ (3.3)

Our focus will be on a subspace of $\mathbb{X}_{\mathtt{w}}$ in which
$({\boldsymbol{\tau}}_{b}|_{e})=({\boldsymbol{\tau}}|_{e}\mathbf{n}_{e})\mathbf{n}_{e}$.
On the other hand, discrete weak divergence operator can be introduced using a
finite-dimensional space $\mathbb{X}_{h}\subset\mathbb{X}_{\mathtt{w}}$, which
will be stated in the next. First, for any mesh object
$\varpi\in\mathcal{K}_{h}\cup\mathcal{E}_{h}$ and for any $r\in\mathrm{N}$ let
us introduce the space $\mathrm{P}_{r}(\varpi)$ to be the space of polynomials
defined on $\varpi$ of degree $\leq r$, with the extended notation
$\mathrm{P}_{-1}(\varpi)=\\{0\\}$. Similarly, we let $\mathbf{P}_{r}(\varpi)$
and $\mathbb{P}_{r}(\varpi)$ be the vectorial and tensorial versions of
$\mathrm{P}_{r}(\varpi)$. Then, given $k\in\mathrm{N}_{0}$, we define for any
$K\in\mathcal{K}_{h}$ the local discrete weak Galerkin space

$\begin{split}\mathbb{X}_{h}(K)\,:=\,\Big{\\{}{\boldsymbol{\tau}}_{h}&=\\{{\boldsymbol{\tau}}_{0h},{\boldsymbol{\tau}}_{bh}\\}\in\mathbb{X}_{\mathtt{w}}(K):\quad{\boldsymbol{\tau}}_{0h}|_{K}\in\mathbb{P}_{k}(K){\quad\hbox{and}\quad}\\\
&~{}~{}~{}~{}~{}~{}{\boldsymbol{\tau}}_{bh}|_{e}=\tau_{b}\otimes\mathbf{n}_{e}\,,\,\,\tau_{b}\in\mathbf{P}_{k}(e)\,,\,\,\forall\,e\subset\partial
K\,,\,\,\forall\,\mathbf{n}_{e}\in\mathrm{D}_{h}\Big{\\}}\,.\end{split}$

In addition, the global finite dimensional space $\mathbb{X}_{h}$, associated
with the partition $\mathcal{K}_{h}$, is defined so that the restriction of
every weak function ${\boldsymbol{\tau}}_{h}$ to the mesh element $K$ belongs
to $\mathbb{X}_{h}(K)$.

######  Definition 3.2 ([31]).

For any ${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}(K)$ and element
$K\in\mathcal{K}_{h}$, the discrete weak divergence operator, denoted by
$\operatorname{\mathbf{div}}_{\mathtt{w,h}}$, is defined as the unique vector-
valued polynomial
$\operatorname{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\in\mathbf{P}_{k+1}(K)$
satisfying

$({\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h}),\,{\boldsymbol{\zeta}}_{h})_{0,K}\,:=\,-({\boldsymbol{\tau}}_{0h},\,\nabla{\boldsymbol{\zeta}}_{h})_{0,K^{0}}\,+\,\langle{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e},{\boldsymbol{\zeta}}_{h}\rangle_{0,\partial
K}\quad\forall\,{\boldsymbol{\zeta}}_{h}\in\mathbf{P}_{k+1}(K)\,.$ (3.4)

On the other hand, for approximating the velocity unknowns we simply consider
the piecewise polynomial space of degree up to $k+1$:

$\mathbf{Y}_{h}\,:=\,\bigg{\\{}\mathbf{v}_{h}\in\mathbf{Y}:\quad\mathbf{v}_{h}|_{K}\in\mathbf{P}_{k+1}(K)\quad\text{for
all}~{}K\in\mathcal{K}_{h}\bigg{\\}}\,.$

L2-orthogonal projections and approximation properties. For any $r\in$ N and
$K\in\mathcal{K}_{h}$, we introduce L2-projection operators
$\mathcal{P}\mathcal{P}_{0,r}^{K}:\mathbb{L}^{2}(K)\rightarrow\mathbb{P}_{r}(K)$
and ${\boldsymbol{\mathcal{P}}}_{b,r}^{K}:\mathbf{L}^{2}(\partial
K)\rightarrow\mathbf{P}_{r}(\partial K)$ which are type of interior and
boundary, and are given by

$\int_{K}\mathcal{P}\mathcal{P}_{0,r}^{K}({\boldsymbol{\tau}}):\widehat{\mathbf{q}}_{r}=\int_{K}{\boldsymbol{\tau}}:\widehat{\mathbf{q}}_{r}\quad{\quad\hbox{and}\quad}\quad\int_{\partial
K}{\boldsymbol{\mathcal{P}}}_{b,r}^{K}(\mathbf{v})\cdot\mathbf{q}_{r}=\int_{\partial
K}\mathbf{v}\cdot\mathbf{q}_{r}\,,$ (3.5)

for all
$({\boldsymbol{\tau}},\mathbf{v})\in\mathbb{L}^{2}(K)\times\mathbf{L}^{2}(\partial
K)$ and
$(\widehat{\mathbf{q}}_{r},\mathbf{q}_{r})\in\mathbb{P}_{r}(K)\times\mathbf{P}_{r}(\partial
K)$.
Now, we introduce projection operator $\mathcal{P}\mathcal{P}_{h}^{K}$ into
the tensorial weak Galerkin space $\mathbb{X}_{h}(K)$ as:

$\mathcal{P}\mathcal{P}_{h}^{K}{\boldsymbol{\tau}}\,:=\,\left\\{\mathcal{P}\mathcal{P}_{0,k}^{K}{\boldsymbol{\tau}}_{0}\,,\,{\boldsymbol{\mathcal{P}}}_{b,k}^{K}(\tau_{b})\otimes\mathbf{n}_{e}\right\\},\qquad\text{for
all}~{}{\boldsymbol{\tau}}\in\mathbb{X}_{h}(K)\,.$

Also, for each element $K\in\mathcal{K}_{h}$ and function
${\boldsymbol{\tau}}\in\mathbb{X}_{h}$, the global projection operator
$\mathcal{P}\mathcal{P}_{h}$ on the space $\mathbb{X}_{h}$ is defined by

$\mathcal{P}\mathcal{P}_{h}({\boldsymbol{\tau}})|_{K}\,=\,\mathcal{P}\mathcal{P}_{h}^{{K}}({\boldsymbol{\tau}}|_{K})\,.$

The approximation properties of $\mathcal{P}\mathcal{P}_{0}$ and
$\boldsymbol{\mathcal{P}}_{h}$ are stated as follows.

######  Lemma 3.3.

Let $\mathcal{K}_{h}$ be a finite element partition of $\Omega$ satisfying the
shape regularity assumptions $\mathbf{A1}$-$\mathbf{A4}$ stated in [32]. Then,
for $k,s,m\in\mathrm{N}_{0}$ such that $m\in\\{0,1\\}$ there exist constants
$C_{1},C_{2}$, independent on the mesh size $h$, such that

$\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\|{\boldsymbol{\tau}}-\mathcal{P}\mathcal{P}_{0}^{K}({\boldsymbol{\tau}})\|_{m,0;K}^{2}$
$\displaystyle\,\leq\,C_{1}h^{2(s-m)}|{\boldsymbol{\tau}}|_{s}^{2}\qquad s\leq
k\,,$ (3.6)
$\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\|\mathbf{v}-\boldsymbol{\mathcal{P}}_{k+1}^{K}(\mathbf{v})\|_{0;K}^{2}$
$\displaystyle\,\leq\,C_{2}h^{2s}|\mathbf{v}|_{s}^{2}\qquad~{}~{}~{}s\leq
k+1\,.$ (3.7)

###  3.2 Weak Galerkin scheme

In order to define our weak Galerkin scheme for Problem 2, we now introduce,
when necessary, the discrete versions of the bilinear forms and functionals
involving the weak spaces. Following the usual procedure in the WG setting,
the construction of them is based on the weak derivatives to ensure
computability for all weak functions.

Notice that, for each
${\boldsymbol{\zeta}},{\boldsymbol{\tau}}\in\mathbb{X}_{h}$ and
${\mathbf{z}},{\mathbf{w}},{\mathbf{v}}\in\mathbf{Y}_{h}$, the quantities

$\mathcal{A}({\boldsymbol{\zeta}},{\boldsymbol{\tau}})\,,\qquad\mathcal{C}({\mathbf{z}},{\mathbf{w}};{\boldsymbol{\tau}})\,,\qquad\mathcal{G}({\boldsymbol{\tau}}){\quad\hbox{and}\quad}\mathcal{F}(\mathbf{v})\,,$

are computable, while the bilinear form
$\mathcal{B}|_{\mathbb{X}_{h}\times\mathbf{Y}_{h}}$ is not computable because
it involve the divergence operator which cannot be evaluated for weak
functions. To overcome this matter, employing Definition 3.2 we define the
discrete bilinear form
$\mathcal{B}_{h}^{K}:\mathbb{X}_{h}(K)\times\mathbf{P}_{k+1}(K)\rightarrow\mathrm{R}$
by

$\mathcal{B}_{h}^{K}({\boldsymbol{\tau}},{\mathbf{v}})\,:=\,\displaystyle\int_{K}{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\cdot{\mathbf{v}}_{h}\qquad\forall\,({\boldsymbol{\tau}},{\mathbf{v}})\in\mathbb{X}_{h}(K)\times\mathbf{P}_{k+1}(K)\,.$
(3.8)

On the other hand, despite the computability of
$\mathcal{A}|_{\mathbb{X}_{h}\times\mathbb{X}_{h}}$, this form needs an
additional stabilizer term to achieve the well-posedness of our weak Galerkin
scheme. More precisely, we define the corresponding discrete bilinear form as
follows:

$\mathcal{A}_{h}^{K}({\boldsymbol{\zeta}},{\boldsymbol{\tau}})\,:=\,\dfrac{1}{\nu}\displaystyle\int_{K}{\boldsymbol{\sigma}}_{0}^{\mathtt{d}}:{\boldsymbol{\tau}}_{0}^{\mathtt{d}}\,+\,\rho\,\mathcal{S}^{K}({\boldsymbol{\zeta}},{\boldsymbol{\tau}})\qquad\forall\,{\boldsymbol{\zeta}},{\boldsymbol{\tau}}\in\mathbb{X}_{h}(K)\,,$
(3.9)

where $\rho$ is the piecewise constant on $\mathcal{K}_{h}$ and the
stabilization form
$\mathcal{S}^{K}(\cdot,\cdot):\mathbb{X}_{h}(K)\times\mathbb{X}_{h}(K)\rightarrow\mathrm{R}$
is given by

$\mathcal{S}^{K}({\boldsymbol{\zeta}},{\boldsymbol{\tau}})\,:=\,h_{K}\left\langle{\boldsymbol{\zeta}}_{0}\mathbf{n}-{\boldsymbol{\zeta}}_{b}\mathbf{n},\,{\boldsymbol{\tau}}_{0}\mathbf{n}-{\boldsymbol{\tau}}_{b}\mathbf{n}\right\rangle_{0,\partial
K}\qquad\forall\,{\boldsymbol{\zeta}},{\boldsymbol{\tau}}\in\mathbb{X}_{h}(K)\,.$
(3.10)

In addition, the global bilinear forms $\mathcal{A}_{h}$ and $\mathcal{B}_{h}$
can be derived by adding the local contributions, that is,

$\mathcal{A}_{h}(\cdot,\cdot)\,:=\,\sum\limits_{K\in\mathcal{K}_{h}}\mathcal{A}_{h}^{K}(\cdot,\cdot)\,{\quad\hbox{and}\quad}\,\mathcal{B}_{h}(\cdot,\cdot)\,:=\,\sum\limits_{K\in\mathcal{K}_{h}}\mathcal{B}_{h}^{K}(\cdot,\cdot)\,.$

Finally, let us introduce the subspace of $\mathbb{X}_{h}$ as

$\mathbb{X}_{0,h}\,:=\,\bigg{\\{}{\boldsymbol{\tau}}_{h}=\\{{\boldsymbol{\tau}}_{0h},{\boldsymbol{\tau}}_{bh}\\}\in\mathbb{X}_{h}:\quad\displaystyle\int_{\Omega}{\mathrm{tr}}({\boldsymbol{\tau}}_{0h})\,=\,0\bigg{\\}}\,.$

Referring to the above space, the discrete bilinear form (3.9) and (3.8), the
discrete weak Galerkin problem reads as follow.

######  Problem 3 (WG problem).

Find ${\boldsymbol{\sigma}}_{h}\in\mathbb{X}_{0,h}$ and
$\mathbf{u}_{h}\in\mathbf{Y}_{h}$ such that

$\left\\{\begin{array}[]{ll}\mathcal{A}_{h}({\boldsymbol{\sigma}}_{h},{\boldsymbol{\tau}}_{h})+\mathcal{C}({\mathbf{u}}_{h},{\mathbf{u}}_{h};{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},{\mathbf{u}}_{h})\,=\,\mathcal{G}({\boldsymbol{\tau}}_{h})\qquad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{0,h}\,,\\\\[8.53581pt]
\hskip
150.79968pt\mathcal{B}_{h}({\boldsymbol{\sigma}}_{h},{\mathbf{v}}_{h})\,=\,\mathcal{F}({\mathbf{v}}_{h})\qquad\forall\,{\mathbf{v}}_{h}\in\mathbf{Y}_{h}\,.\end{array}\right.$

##  4 Solvability analysis

This section goals to examine the solvability of Problem 3. Initially, we
delve into the discussion of the stability properties inherent in the discrete
forms of Problem 3. Subsequently, we employ the Banach fixed-point and
classical Banach–Nečas–Babuška theorems to establish the well-posedness of the
discrete scheme, assuming an appropriate smallness condition on the data.

###  4.1 Stability properties

We begin by equipping the approximate pair spaces $\mathbb{X}_{h}$ and
$\mathbf{Y}_{h}+\mathbf{H}^{1}(\Omega)$ with the following discrete norms,
respectively (see e.g., [32])

$\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}^{2}\,:=\,\sum_{K\in\mathcal{K}_{h}}\big{[}\|{\boldsymbol{\tau}}_{0h}\|_{0,K}^{2}+\mathcal{S}^{K}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\big{]}\qquad\text{for
all}~{}{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,,$

and

$\|\mathbf{v}_{h}\|_{1,h}^{2}\,:=\,\sum_{K\in\mathcal{K}_{h}}|\mathbf{v}_{h}|_{1,K}^{2}+\sum_{e\in\mathcal{E}_{h}}h_{e}^{-1}\|{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket\|_{0,e}^{2}\qquad\quad\text{for
all}~{}\mathbf{v}_{h}\in\mathbf{Y}_{h}+\mathbf{H}^{1}(\Omega)\,,$

Next, we provide the following result which is a counterpart of [8, Lemma 3.1]
for the weak divergence operator.

######  Lemma 4.1.

There exists $\hat{c}_{\Omega}>0$ depends on $\Omega$ but independent of mesh
size, such that

$\hat{c}_{\Omega}\|{\boldsymbol{\tau}}_{0h}\|_{0}^{2}\,\leq\,\|{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}\|_{0}^{2}+\|\operatorname{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\|_{0,4/3}^{2}\quad\quad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,.$
(4.1)

###### Proof.

See [22, Lemma 3.1]. ∎

Then, some properties of $\mathcal{A}_{h}$ (cf. (3.9)) are established as
follows.

######  Lemma 4.2.

The discrete bilinear form $\mathcal{A}_{h}$ defined in (3.9), satisfy the
following properties:

  * •

consistency: for any $\boldsymbol{\xi}\in\mathbb{H}^{1}(\Omega)$, we have that

$\mathcal{A}_{h}({\boldsymbol{\zeta}}_{h},\boldsymbol{\xi})\,=\,\mathcal{A}({\boldsymbol{\zeta}}_{h},\boldsymbol{\xi})\qquad\forall\,{\boldsymbol{\zeta}}_{h}\in\mathbb{X}_{h}\,.$

  * •

stability and boundedness: there exists positive constant $c_{\mathcal{A}}$,
independent of $K$ and $h$, such that:

$\big{|}\mathcal{A}_{h}({\boldsymbol{\zeta}}_{h},{\boldsymbol{\tau}}_{h})\big{|}\,\leq\,c_{\mathcal{A}}\,\|{\boldsymbol{\zeta}}_{h}\|_{\mathbb{H},h}\,\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}\qquad\forall\,{\boldsymbol{\zeta}}_{h},{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,,$
(4.2)

and let $\widetilde{\mathbb{X}}_{h}$ be the discrete kernel of the bilinear
form $\mathcal{B}_{h}$. Then, there exists constant $\alpha>0$, independent of
$K$, such that

$\mathcal{A}_{h}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\,\geq\,\alpha\,\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}^{2}\qquad\forall\,{\boldsymbol{\tau}}_{h}\in\widetilde{\mathbb{X}}_{h}\,.$
(4.3)

###### Proof.

By considering $\boldsymbol{\xi}\in\mathbb{H}^{1}(\Omega)$ and
${\boldsymbol{\zeta}}_{h}\in\mathbb{X}_{h}$ and employing the definition of
the discrete form $\mathcal{A}_{h}$ (cf. (3.9)), along with the observation
that the stabilization term (cf. (3.10)) vanishes when one of the components
is sufficiently regular, we deduce the ”consistency” property. Next, to verify
the boundedness of $\mathcal{A}_{h}$, we use the Cauchy–Schwarz inequality and
virtue of

$\|{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}\|_{0,K}^{2}\,=\,\|{\boldsymbol{\tau}}_{0h}\|_{0,K}^{2}-\dfrac{1}{2}\|\operatorname{tr}({\boldsymbol{\tau}}_{0h})\|_{0,K}^{2}\,\leq\,\|{\boldsymbol{\tau}}_{0h}\|_{0,K}^{2}\,.$
(4.4)

This gives

$\displaystyle\big{|}\mathcal{A}_{h}({\boldsymbol{\zeta}}_{h},{\boldsymbol{\tau}}_{h})\big{|}$
$\displaystyle\,=\,\bigg{|}\sum_{K\in\mathcal{K}_{h}}\Big{(}\dfrac{1}{\nu}\displaystyle\int_{K}{\boldsymbol{\zeta}}_{0h}^{\mathtt{d}}:{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}+\rho\,\mathcal{S}^{K}({\boldsymbol{\zeta}}_{h},{\boldsymbol{\tau}}_{h})\Big{)}\bigg{|}$
$\displaystyle\,\leq\,\sum_{K\in\mathcal{K}_{h}}\Big{(}\dfrac{1}{\nu}\|{\boldsymbol{\zeta}}_{0h}\|_{0,K}\,\|{\boldsymbol{\tau}}_{0h}\|_{0,K}+\rho\left(\mathcal{S}^{K}({\boldsymbol{\zeta}}_{h},{\boldsymbol{\zeta}}_{h})\right)^{1/2}\left(\mathcal{S}^{K}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\right)^{1/2}\Big{)}$
$\displaystyle\,\leq\,\max\left\\{\dfrac{1}{\nu},\rho\right\\}\|{\boldsymbol{\zeta}}_{h}\|_{\mathbb{H},h}\,\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}\,.$

We now aim to establish the ellipticity of $\mathcal{A}_{h}$ on the discrete
kernel of $\mathcal{B}_{h}|_{\mathbb{X}_{h}\times\mathbf{Y}_{h}}$, that is,

$\widetilde{\mathbb{X}}_{h}\,:=\,\bigg{\\{}{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}:\quad\int_{\Omega}\operatorname{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\cdot\mathbf{v}_{h}\,=\,0\quad\forall\,\mathbf{v}_{h}\in\mathbf{Y}_{h}\bigg{\\}}\,,$

which, together with fact
$\operatorname{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\in\mathbf{Y}_{h}$
(cf. Definition 3.2) implies that

$\widetilde{\mathbb{X}}_{h}\,:=\,\Big{\\{}{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}:\quad{\mathbf{div}}_{\mathtt{w,h}}({\boldsymbol{\tau}}_{h})\,=\,0\quad\text{in}\,\,\Omega\,\Big{\\}}\,.$
(4.5)

Therefore, employing Lemma 4.1 and the fact that ${\boldsymbol{\tau}}_{h}$ is
weak divergence-free, we obtain

$\displaystyle\mathcal{A}_{h}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})$
$\displaystyle\,=\,\sum_{K\in\mathcal{K}_{h}}\Big{(}\dfrac{1}{\nu}\|{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}\|_{0,K}^{2}+\rho\,\mathcal{S}^{K}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\Big{)}$
$\displaystyle\,\geq\,\sum_{K\in\mathcal{K}_{h}}\Big{(}\dfrac{\hat{c}_{\Omega}}{\nu}\|{\boldsymbol{\tau}}_{0h}\|_{0,K}^{2}+\rho\,\mathcal{S}^{K}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\Big{)}$
$\displaystyle\,\geq\,\min\left\\{\dfrac{\hat{c}_{\Omega}}{\nu},\rho\right\\}\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}^{2}\,,$

which together with setting
$c_{\mathcal{A}}\,=\,\max\left\\{\dfrac{1}{\nu},\rho\right\\}$ and
$\alpha\,=\,\min\left\\{\dfrac{1}{\nu},\rho\right\\}$ completes the proof. ∎

To establish the discrete inf-sup condition for the bilinear form
$\mathcal{B}_{h}$, we require a preliminary result, as stated in the following
lemma.

######  Lemma 4.3.

For any ${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}$ and
$\mathbf{v}_{h}\in\mathbf{Y}_{h}$, we have

$\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\,=\,\sum_{e\in\mathcal{E}_{h}}\langle{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e},\llbracket\mathbf{v}_{h}\rrbracket\rangle_{e}-\sum_{K\in\mathcal{K}_{h}}({\boldsymbol{\tau}}_{0h},\boldsymbol{\nabla}\mathbf{v}_{h})_{0,K}\,.$

###### Proof.

It straightforwardly follows from the definition of $\mathcal{B}_{h}$, as
given in (3.8), and the application of the discrete divergence operator (cf.
Definition 3.2). ∎

We are now in a position to establish the discrete inf-sup condition for the
bilinear from $\mathcal{B}_{h}$.

######  Lemma 4.4.

There exists a positive constant $\widehat{\beta}$ independent of $h$ such
that

$\sup_{\mathbf{0}\neq{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}}\dfrac{\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})}{\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}}\,\geq\,\widehat{\beta}\,\|\mathbf{v}_{h}\|_{1,h}\qquad\forall\,\mathbf{v}_{h}\in\mathbf{Y}_{h}\,.$
(4.6)

###### Proof.

In what follows, we proceed similarly to the proof of [32, Lemma 3.3]. In
fact, given $\mathbf{v}_{h}\in\mathbf{Y}_{h}$, we set

$\left\\{\begin{array}[]{l}{\boldsymbol{\tau}}_{0h}\,=\,-\boldsymbol{\nabla}\mathbf{v}_{h}\qquad\qquad\qquad~{}\text{in}~{}K\,,\\\\[5.69054pt]
{\boldsymbol{\tau}}_{bh}\,=\,h_{e}^{-1}{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket\otimes\mathbf{n}_{e}\qquad\text{on}~{}e\,,\end{array}\right.$

which satisfies

$\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\,=\,\sum_{e\in\mathcal{E}_{h}}\langle{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket,{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket\rangle_{e}+\sum_{K\in\mathcal{K}_{h}}(\boldsymbol{\nabla}\mathbf{v}_{h},\boldsymbol{\nabla}\mathbf{v}_{h})_{0,K}\,=\,\|\mathbf{v}_{h}\|_{1,h}^{2}\,.$
(4.7)

On the other hand, we use the definition of the discrete norm
$\|\cdot\|_{\mathbb{H},h}$ to the above chosen
${\boldsymbol{\tau}}_{h}\,:=\,\\{{\boldsymbol{\tau}}_{0h},{\boldsymbol{\tau}}_{bh}\\}$
and the trace inequlaity, to obtain

$\displaystyle\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}^{2}$
$\displaystyle\,=\,\sum_{K\in\mathcal{K}_{h}}\Big{(}\|\boldsymbol{\nabla}\mathbf{v}_{h}\|_{0,K}^{2}+\mathcal{S}^{K}({\boldsymbol{\tau}}_{h},{\boldsymbol{\tau}}_{h})\Big{)}$
$\displaystyle\,=\,\sum_{K\in\mathcal{K}_{h}}\Big{(}\|\boldsymbol{\nabla}\mathbf{v}_{h}\|_{0,K}^{2}+h_{K}\|(\boldsymbol{\nabla}\mathbf{v}_{h})\mathbf{n}+h_{e}^{-1}{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket\|_{0,\partial
K}^{2}\Big{)}$
$\displaystyle\,\leq\,(1+\mathcal{C}_{\mathtt{tr}})\sum_{K\in\mathcal{K}_{h}}\|\boldsymbol{\nabla}\mathbf{v}_{h}\|_{0,K}^{2}+\sum_{e\in\mathcal{E}_{h}}h_{e}^{-1}\|{\boldsymbol{\mathcal{P}}}_{b,k}\llbracket\mathbf{v}_{h}\rrbracket\|_{0,e}^{2}\,\leq\,\max\\{1,(1+\mathcal{C}_{\mathtt{tr}})\\}\|\mathbf{v}_{h}\|_{1,h}^{2}\,.$

Combining the above result with (4.7) concludes the desired inequality (4.6)
with setting
$\widehat{\beta}:=\frac{1}{\sqrt{\max\\{1,(1+\mathcal{C}_{\mathtt{tr}})\\}}}$.
∎

The boundedness properties of the discrete bilinear form $\mathcal{B}_{h}$
(cf. (3.8)) is established as follows.

######  Lemma 4.5.

The discrete bilinear form $\mathcal{B}_{h}(\cdot,\cdot)$ is bounded in
$\mathbb{X}_{h}\times\mathbf{Y}_{h}$. In other words, there exists a positive
constant $c_{\mathcal{B}}$ such that

$\big{|}\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\big{|}\,\leq\,c_{\mathcal{B}}\,\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}\,\|\mathbf{v}_{h}\|_{1,h}\qquad\forall\,({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}\,.$
(4.8)

###### Proof.

An application of Lemma 4.3 and Cauchy–Schwarz inequality, yields

$\begin{array}[]{c}\big{|}\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\big{|}\,\leq\,\displaystyle\sum_{e\in\mathcal{E}_{h}}\|{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e}\|_{0,e}\,\|\llbracket\mathbf{v}_{h}\rrbracket\|_{0,e}+\displaystyle\sum_{K\in\mathcal{K}_{h}}\|{\boldsymbol{\tau}}_{0h}\|_{0,K}\,\|\boldsymbol{\nabla}\mathbf{v}_{h}\|_{0,K}\\\\[8.61108pt]
\,\leq\,\left(\displaystyle\sum_{e\in\mathcal{E}_{h}}h_{e}\|{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e}\|_{0,e}^{2}\right)^{1/2}\left(\displaystyle\sum_{e\in\mathcal{E}_{h}}h_{e}^{-1}\|\llbracket\mathbf{v}_{h}\rrbracket\|_{0,e}^{2}\right)^{1/2}+\|{\boldsymbol{\tau}}_{0h}\|_{0,\Omega}\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}\|\boldsymbol{\nabla}\mathbf{v}_{h}\|_{0,K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,\bigg{\\{}\left(\displaystyle\sum_{e\in\mathcal{E}_{h}}h_{e}\|{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e}\|_{0,e}^{2}\right)^{1/2}+\|{\boldsymbol{\tau}}_{0h}\|_{0,\Omega}\bigg{\\}}\,\|\mathbf{v}_{h}\|_{1,h}\,.\end{array}$
(4.9)

Our aim is now to show that

$\left(\sum_{e\in\mathcal{E}_{h}}h_{e}\|{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e}\|_{0,e}^{2}\right)^{1/2}+\|{\boldsymbol{\tau}}_{0h}\|_{0,\Omega}\,\lesssim\,\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}\qquad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,.$

To attain the aforementioned inequality, we employ the trace inequality,
resulting

$\displaystyle\sum_{e\in\mathcal{E}_{h}}h_{e}\|{\boldsymbol{\tau}}_{bh}\mathbf{n}_{e}\|_{0,e}^{2}$
$\displaystyle\,\leq\,2\sum_{e\in\mathcal{E}_{h}}\Big{(}h_{e}\|({\boldsymbol{\tau}}_{bh}-{\boldsymbol{\tau}}_{0h})\mathbf{n}_{e}\|_{0,e}^{2}+h_{e}\|{\boldsymbol{\tau}}_{0h}\mathbf{n}_{e}\|_{0,e}^{2}\Big{)}$
$\displaystyle\,\leq\,2\sum_{e\in\mathcal{E}_{h}}h_{e}\|({\boldsymbol{\tau}}_{bh}-{\boldsymbol{\tau}}_{0h})\mathbf{n}_{e}\|_{0,e}^{2}+2\,\mathcal{C}_{\mathtt{tr}}^{2}\sum_{K\in\mathcal{K}_{h}}\|{\boldsymbol{\tau}}_{0h}\|_{0,K}^{2}\,.$

Substituting the above result back into (4.9) finishes the proof of lemma. ∎

Now, employing the aforementioned stability properties of $\mathcal{A}_{h}$
(cf. 4.2, eqs. (4.2) and (4.3)), the discrete inf-sup condition and the
boundedness of $\mathcal{B}_{h}$ (cf. Lemmas 4.4 and 4.5), and applying [14,
Proposition 2.36] it is not difficult to see that the bilinear form
$\Lambda:(\mathbb{X}_{h}\times\mathbf{Y}_{h})\times(\mathbb{X}_{h}\times\mathbf{Y}_{h})\rightarrow\mathrm{R}$
defined by

$\Lambda_{h}[({\boldsymbol{\zeta}}_{h},\textbf{w}_{h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,:=\,\mathcal{A}_{h}({\boldsymbol{\zeta}}_{h},{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\textbf{w}_{h})-\mathcal{B}_{h}({\boldsymbol{\zeta}}_{h},{\mathbf{v}}_{h})\,,$
(4.10)

satisfies:

$\sup_{\textbf{0}\neq({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\Lambda_{h}[({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h}),({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})]}{\|({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\|_{\mathbb{X}\times\mathbf{Y}}}\,\geq\,\alpha_{\Lambda,\mathtt{d}}\,\|({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h})\|_{\mathbb{X}\times\mathbf{Y}}\,,$
(4.11)

where $\alpha_{\Lambda,\mathtt{d}}$ is the positive constant dependent on
$c_{\mathcal{A}},c_{\mathcal{B}},\alpha,\widehat{\beta}$.

Furthermore, the boundedness of the trilinear form
$\mathcal{C}(\cdot;\cdot,\cdot)$ on
$\mathbf{Y}_{h}\times\mathbf{Y}_{h}\times\mathbb{X}_{h}$ can be readily
inferred by utilizing the Hölder inequality, the Sobolev embedding
$\mathbf{H}^{1}\subset\mathbf{L}^{4}$ and inequality (4.4), as follows:

$\begin{array}[]{c}\big{|}\mathcal{C}(\textbf{w}_{h};{\mathbf{v}}_{h},{\boldsymbol{\tau}}_{h})\big{|}\,=\,\big{|}\dfrac{1}{\nu}\left((\textbf{w}_{h}\otimes{\mathbf{v}}_{h})^{\mathtt{d}},{\boldsymbol{\tau}}_{0h}\right)_{0,\Omega}\big{|}\\\\[8.61108pt]
\,=\,\big{|}\dfrac{1}{\nu}\left((\textbf{w}_{h}\otimes{\mathbf{v}}_{h}),{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}\right)_{0,\Omega}\big{|}\,\leq\,\dfrac{1}{\nu}\|\textbf{w}_{h}\|_{\mathbf{Y}}\|{\mathbf{v}}_{h}\|_{\mathbf{Y}}\|{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}\|_{0,\Omega}\\\\[8.61108pt]
\,\leq\,\dfrac{c_{\mathtt{em}}^{2}}{\nu}\|\textbf{w}_{h}\|_{1,h}\|{\mathbf{v}}_{h}\|_{1,h}\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}\,.\end{array}$
(4.12)

###  4.2 The fixed-point strategy

We begin by introducing the associated fixed-point operator for any
$\textbf{z}_{h}\in\mathbf{Y}_{h}$ as
$\mathbf{S}_{\mathtt{d}}(\textbf{z}_{h})\,=\,{\mathbf{u}}_{\star,h}$, where
$({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h})\in\mathbb{X}_{0,h}\times\mathbf{Y}_{h}$
is the solution of the linearized version of Problem 3, that is,

$\left\\{\begin{array}[]{rcll}\mathcal{A}_{h}({\boldsymbol{\sigma}}_{\star,h},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\textbf{z}_{h},{\mathbf{u}}_{\star,h};{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},{\mathbf{u}}_{\star,h})&=&\,\mathcal{G}({\boldsymbol{\tau}}_{h})&\quad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{0,h}\,,\\\\[8.53581pt]
\mathcal{B}_{h}({\boldsymbol{\sigma}}_{\star,h},{\mathbf{v}}_{h})&=&\,\mathcal{F}({\mathbf{v}}_{h})&\quad\forall\,{\mathbf{v}}_{h}\in\mathbf{Y}_{h}\,.\end{array}\right.$
(4.13)

Equivalently, given ${\mathbf{z}}_{h}\in\mathbf{Y}_{h}$, introducing the
bilinear form
$\Lambda_{h,{\mathbf{z}}_{h}}:(\mathbb{X}_{h}\times\mathbf{Y}_{h})\times(\mathbb{X}_{h}\times\mathbf{Y}_{h})\rightarrow\mathrm{R}$
given by

$\Lambda_{h,{\mathbf{z}}_{h}}[({\boldsymbol{\zeta}}_{h},\textbf{w}_{h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,:=\,\Lambda_{h}[({\boldsymbol{\zeta}}_{h},\textbf{w}_{h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,+\,\mathcal{C}(\textbf{z}_{h};\textbf{w}_{h},{\mathbf{v}}_{h})\,,$
(4.14)

and the linear functional
$\mathrm{F}:\mathbb{X}_{h}\times\mathbf{Y}_{h}\rightarrow\mathrm{R}$ as

$\mathrm{F}({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\,:=\,\mathcal{G}({\boldsymbol{\tau}}_{h})+\mathcal{F}({\mathbf{v}}_{h})\,.$
(4.15)

for all
$({\boldsymbol{\zeta}}_{h},{\mathbf{w}}_{h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}$,
the linearized problem (4.13) can be rewritten as

$\Lambda_{\textbf{z}_{h}}[({\boldsymbol{\sigma}}_{\star,h},\textbf{u}_{\star,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,=\,\mathrm{F}({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\qquad\forall\,({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\in\mathbb{X}_{0,h}\times\mathbf{Y}_{h}\,.$
(4.16)

It can be observed that solving Problem 3 is equivalent to seeking a fixed
point of $\textbf{S}_{\mathtt{d}}$, that is: Find
$\textbf{u}_{h}\in\mathbf{Y}_{h}$ such that

$\textbf{S}_{\mathtt{d}}(\textbf{u}_{h})\,=\,\textbf{u}_{h}\,.$

Next, we utilize the classical Banach–Nečas–Babuška theorem (see, for instance
[14, Theorem 2.6]) to demonstrate that, for any arbitrary
$\textbf{z}_{h}\in\mathbf{Y}_{h}$, the problem (4.13) (or equivalently (4.16))
is well-posed, implying the well-definedness of $\mathbf{S}_{\mathtt{d}}$.

The following lemma demonstrates the global discrete inf-sup condition for
form $\Lambda_{{\mathbf{z}}_{h}}$ on pair-space
$\mathbb{X}_{h}\times\mathbf{Y}_{h}$.

######  Lemma 4.6.

For any $\mathbf{z}_{h}\in\mathbf{Y}_{h}$ such that

$\|\mathbf{z}_{h}\|_{1,h}\,\leq\,\dfrac{\nu\,\alpha_{\Lambda,\mathtt{d}}}{2c_{\mathtt{em}}^{2}}\,,$
(4.17)

there holds

$\sup_{\mathbf{0}\neq({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\Lambda_{h,{\mathbf{z}}_{h}}[({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h}),({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})]}{\|({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\|_{h}}\,\geq\,\dfrac{\alpha_{\Lambda,\mathtt{d}}}{2}\,\|({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h})\|_{h}\qquad\forall\,({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}\,.$
(4.18)

where $\alpha_{\Lambda,\mathtt{d}}$ is the positive constant in (4.11).

###### Proof.

Bearing in mind the definition of $\Lambda_{h,{\mathbf{z}}_{h}}$ (cf. (4.14))
for each ${\mathbf{z}}_{h}\in\mathbf{Y}_{h}$, and combining (4.11) with the
effect of the extra term given by $\mathcal{C}({\mathbf{z}}_{h};\cdot,\cdot)$,
which means invoking the upper bound provided by (4.12), we arrive at

$\sup_{\mathbf{0}\neq({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\Lambda_{h,{\mathbf{z}}_{h}}[({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h}),({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})]}{\|({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\|_{h}}\,\geq\,\Big{\\{}\alpha_{\Lambda,\mathtt{d}}-\dfrac{c_{\mathtt{em}}^{2}}{\nu}\,\|{\mathbf{z}}_{h}\|_{\mathbf{Y}}\Big{\\}}\|(\boldsymbol{\zeta}_{h},{\mathbf{w}}_{h})\|_{h}\quad\forall~{}({\boldsymbol{\zeta}}_{h},{\mathbf{w}}_{h})\in\mathcal{V}_{h},$

from which, under the assumption (4.17) we arrive at (4.18), which ends the
proof. ∎

Now we ready to show that mapp $\mathbf{S}_{\mathtt{d}}$ is well-defined or
equivalently problem (4.13) is uniquely solvable.

######  Lemma 4.7.

Let the assumption of Lemma 4.6 be satisfied. Then, there exists a unique
$({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}$
solution to problem (4.13). In addition, there holds

$\|\mathbf{S}_{\mathtt{d}}(\mathbf{z}_{h})\|_{1,h}\,=\,\|{\mathbf{u}}_{\star,h}\|_{1,h}\,\leq\,\dfrac{2}{\alpha_{\Lambda,\mathtt{d}}}\left(\|\mathbf{f}\|_{0,4/3}+c_{g}\|\mathbf{g}\|_{1/2,\partial\Omega}\right)\,.$
(4.19)

###### Proof.

A straightforward application of the classical Babuška–Brezzi theory and Lemma
4.6 implies that problem (4.13) is well-posed. For the second part of proof,
by combining (4.16) and Lemma 4.6 with considering
$({\boldsymbol{\zeta}}_{h},\mathbf{w}_{h})\,:=\,({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h})$,
we readily obtain

$\displaystyle\dfrac{\alpha_{\Lambda,\mathtt{d}}}{2}\,\|({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h})\|_{h}$
$\displaystyle\,\leq\,\sup_{\mathbf{0}\neq({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\Lambda_{h,{\mathbf{z}}_{h}}[({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h}),({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})]}{\|({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\|_{h}}$
$\displaystyle\,=\,\sup_{\mathbf{0}\neq({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\mathrm{F}({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})}{\|({\boldsymbol{\tau}}_{h},\mathbf{v}_{h})\|_{h}}\,\leq\,\left(c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}+\|\mathbf{f}\|_{0,4/3}\right)\,,$

where the boundness of $\mathcal{G}$ and $\mathcal{F}$ was used in the last
step, and this finishes the proof. ∎

We continue the analysis by establishing sufficient conditions under which
$\mathbf{S}_{\mathtt{d}}$ maps a closed ball of $\mathbf{Y}_{h}$ into itself.
Let us define the set

$\widehat{\mathbf{Y}}_{h}\,:=\,\bigg{\\{}\textbf{z}_{h}\in\mathbf{Y}_{h}:\quad\|\textbf{z}_{h}\|_{1,h}\,\leq\,\dfrac{\nu\,\alpha_{\Lambda,\mathtt{d}}}{2c_{\mathtt{em}}^{2}}\bigg{\\}}\,,$

and state the following result.

######  Lemma 4.8.

Assume that the data are sufficiently small so that

$\left(c_{g}\|\mathbf{g}\|_{1/2,\partial\Omega}+\|\mathbf{f}\|_{0,4/3}\right)\,\leq\,\dfrac{\nu\,\alpha_{\Lambda,\mathtt{d}}^{2}}{4c_{\mathtt{em}}^{2}}\,,$

Then
$\mathbf{S}_{\mathtt{d}}(\widehat{\mathbf{Y}}_{h})\subseteq\widehat{\mathbf{Y}}_{h}$.

###### Proof.

It deduces straightly from a priori estimate stated by (4.19). ∎

We now address the continuity properties of $\mathbf{S}_{\mathtt{d}}$.

######  Lemma 4.9.

For any $\mathbf{z}_{h},\mathbf{y}_{h}\in\mathbf{Y}_{h}$, there holds

$\|\mathbf{S}_{\mathtt{d}}(\mathbf{z}_{h})-\mathbf{S}_{\mathtt{d}}(\mathbf{y}_{h})\|_{1,h}\,\leq\,\dfrac{4c_{\mathtt{em}}^{2}}{\nu\,\alpha_{\Lambda,\mathtt{d}}^{2}}\,\left(\|\mathbf{f}\|_{0,3/4}+c_{g}\|\mathbf{g}\|_{1/2,\partial\Omega}\right)\|\mathbf{z}_{h}-\mathbf{y}_{h}\|_{1,h}\,.$
(4.20)

###### Proof.

Given ${\mathbf{z}}_{h},{\mathbf{y}}_{h}\in\widehat{\mathbf{Y}}_{h}$ we let
$\mathbf{S}_{\mathtt{d}}(\textbf{z}_{h})={\mathbf{u}}_{\star,h}$ and
$\mathbf{S}_{\mathtt{d}}(\textbf{y}_{h})={\mathbf{u}}_{\circ,h}$, where
$({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h})$ and
$({\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\circ,h})$ are the
corresponding solutions of equation (4.16). It follows from (4.16) that

$\Lambda_{h,\textbf{z}_{h}}[({\boldsymbol{\sigma}}_{\star,h},{\mathbf{u}}_{\star,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,=\,\Lambda_{h,\textbf{y}_{h}}[({\boldsymbol{\sigma}}_{\circ,h},{\mathbf{y}}_{\circ,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,,\quad\quad\forall\,({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}\,,$

from which, according to the definitions of $\Lambda_{h,{\mathbf{z}}_{h}}$ and
$\Lambda_{h,{\mathbf{y}}_{h}}$ (cf. (4.14)), we have

$\Lambda_{h}[({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,=\,\mathcal{C}(\textbf{y}_{h},{\mathbf{u}}_{\circ,h},{\boldsymbol{\tau}}_{h})-\mathcal{C}(\textbf{z}_{h},{\mathbf{u}}_{\star,h},{\boldsymbol{\tau}}_{h})\,.$

This result, combined with (4.10) by setting
$({\boldsymbol{\zeta}}_{h},\textbf{w}_{h}):=({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h})$,
yields

$\begin{array}[]{c}\Lambda_{h,\textbf{z}_{h}}[({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,=\,\Lambda_{h}[({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\\\\[8.61108pt]
\,+\,\mathcal{C}(\textbf{z}_{h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h};{\boldsymbol{\tau}}_{h})\\\\[8.61108pt]
\,=\,\mathcal{C}(\textbf{y}_{h},{\mathbf{u}}_{\circ,h};{\boldsymbol{\tau}}_{h})-\mathcal{C}(\textbf{z}_{h},{\mathbf{u}}_{\star,h};{\boldsymbol{\tau}}_{h})+\mathcal{C}(\textbf{z}_{h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h};{\boldsymbol{\tau}}_{h})\\\\[8.61108pt]
\,=\,-\mathcal{C}({\mathbf{z}}_{h}-{\mathbf{y}}_{h},{\mathbf{u}}_{\circ,h};{\boldsymbol{\tau}}_{h})\,.\end{array}$

Now, we apply the discrete global inf–sup condition (4.18) to the left-hand
side of the above equation and utilize the estimates (4.12) and (4.19), to get

$\displaystyle\dfrac{\alpha_{\Lambda,\mathtt{d}}}{2}\|({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h})\|_{h}$
$\displaystyle\,\leq\,\sup_{\textbf{0}\neq({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\Lambda_{\textbf{z}_{h}}[({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]}{\|({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\|_{h}}$
$\displaystyle\,=\,\sup_{\textbf{0}\neq({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}}\dfrac{\mathcal{C}(\textbf{y}_{h}-\textbf{z}_{h},{\mathbf{u}}_{\circ,h};{\boldsymbol{\tau}}_{h})}{\|({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})\|_{h}}$
$\displaystyle\,\leq\,\dfrac{c_{\mathtt{em}}^{2}}{\nu}\,\|\textbf{z}_{h}-\textbf{y}_{h}\|_{1,h}\,\|{\mathbf{u}}_{\circ,h}\|_{1,h}$
$\displaystyle\,\leq\,\dfrac{2c_{\mathtt{em}}^{2}}{\nu\alpha_{\Lambda,\mathtt{d}}}\left(\|\textbf{f}\|_{0,3/4}+c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}\right)\|\textbf{z}_{h}-\textbf{y}_{h}\|_{1,h}\,.$

So, thanks to the above inequality, we arrive at

$\displaystyle\|\mathbf{S}_{\mathtt{d}}(\textbf{z}_{h})-\mathbf{S}_{\mathtt{d}}(\textbf{y}_{h})\|_{1,h}$
$\displaystyle\,=\,\|{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h}\|_{1,h}$
$\displaystyle\,\leq\,\|({\boldsymbol{\sigma}}_{\star,h}-{\boldsymbol{\sigma}}_{\circ,h},{\mathbf{u}}_{\star,h}-{\mathbf{u}}_{\circ,h})\|_{h}\,\leq\,\dfrac{4}{\alpha_{\Lambda,\mathtt{d}}^{2}\nu}\left(\|\textbf{f}\|_{0,3/4}+c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}\right)\|\textbf{z}_{h}-\textbf{y}_{h}\|_{1,h}\,,$

and thus ends the proof. ∎

The main result of this section is summarized in the following theorem.

######  Theorem 4.10.

Assume that the hypothesis of Lemma 4.8 hold. Then, there exists
$({\boldsymbol{\sigma}}_{h},\mathbf{u}_{h})\in\mathbb{X}_{0,h}\times\mathbf{Y}_{h}$
solution for Problem 3. Moreover, there holds

$\|({\boldsymbol{\sigma}}_{h},\mathbf{u}_{h})\|_{h}\,\leq\,\dfrac{2}{\alpha_{\Lambda,\mathtt{d}}}\,\left(\|\mathbf{f}\|_{0,4/3}+c_{g}\|\mathbf{g}\|_{1/2,\partial\Omega}\right)\,.$
(4.21)

###### Proof.

The compactness of $\mathbf{S}_{\mathtt{d}}$ on space
$\widehat{\mathbf{Y}}_{h}$ and the Lipschitz-continuity of
$\mathbf{S}_{\mathtt{d}}$ are guaranteed by Lemmas 4.8 and 4.9, respectivly.
Hence, by applying the Banach fixed-point theorem directly to Problem 3, one
can conclude its existence and uniqueness. Furthermore, the stability result
(4.21) is derived straightly from (4.19) stated in Lemma 4.7. ∎

##  5 Convergence Analysis

In this section, we are interested in deriving an a priori estimate for the
error

$\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|_{\mathbb{H},h}\,,\quad\|\mathbf{u}-\mathbf{u}_{h}\|_{1,h}\,,{\quad\hbox{and}\quad}\|\mathbf{u}-\mathbf{u}_{h}\|_{0,\Omega}.$

where $({\boldsymbol{\sigma}},{\mathbf{u}})\in\mathbb{X}\times\mathbf{Y}$ with
${\mathbf{u}}\in\widehat{\mathbf{Y}}$, is the unique solution of Problem 2,
and
$({\boldsymbol{\sigma}}_{h},{\mathbf{u}}_{h})\in\mathbb{X}_{0,h}\times\mathbf{Y}_{h}$
with ${\mathbf{u}}_{h}\in\widehat{\mathbf{Y}}_{h}$ is the unique solution of
Problem 3. As a byproduct of this, we also derive an a priori estimate for
$\|p-p_{h}\|_{0,\Omega}$, where $p_{h}$ is the discrete pressure computed
according to the postprocessing formula suggested by (2.4), that is

$p_{h}\,:=\,-\dfrac{1}{2}\left(\operatorname{tr}({\boldsymbol{\sigma}}_{0h})+\operatorname{tr}(\mathbf{u}_{h}\otimes\mathbf{u}_{h})\right)-r_{\mathbf{u}_{h}}\,.$
(5.1)

To this end, we divide our result into suboptimal and optimal convergence for
the velocity.

###  5.1 Suboptimal convergence

We begin by introducing the following lemmas, which will play an essential
role in the error analysis.

######  Lemma 5.1 (trace inequality).

Assume that the partition $\mathcal{K}_{h}$ satisfies the shape regularity
assumptions $\mathbf{A1}$-$\mathbf{A4}$ stated in [32]. Then, there exists a
constant $\mathcal{C}_{\mathtt{tr}}$ such that for any $K\in\mathcal{K}_{h}$
and edge/face $e\subset\partial K$, we have

$\|v\|_{0,e}^{2}\,\leq\,\mathcal{C}_{\mathtt{tr}}\left(h_{K}^{-1}\|v\|_{0,K}^{2}+h_{K}\|\nabla
v\|_{0,K}^{2}\right)\,.$ (5.2)

######  Lemma 5.2 (inverse inequality).

Let $K\subset\mathrm{R}^{d}$ be a d-simplex which has diameter $h_{K}$ and is
shape regular. Assume that $p$ and $q$ be non-negative integers such that
$p\leq q$. Then, there exists a constant $\mathcal{C}_{\mathtt{inv}}$ such
that

$\|\varrho\|_{0,q,K}\,\leq\,\mathcal{C}_{\mathtt{inv}}\,h^{d(\frac{1}{q}-\frac{1}{p})}\,\|\varrho\|_{0,p,K}\,,$
(5.3)

for any polynomial $\varrho$ of degree no more than $k$.

######  Lemma 5.3.

For ${\boldsymbol{\sigma}}\in\mathbb{H}^{r+1}(\Omega)$ and $r\in N$ such that
$r\leq k$, there exist positive constants $C_{0}$ and $C_{\mathtt{s}}$ such
that

$\sum\limits_{K\in\mathcal{K}_{h}}h_{K}\|{\boldsymbol{\sigma}}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n})\|_{0,\partial
K}^{2}\leq C_{0}\,h^{2(r+1)}\|{\boldsymbol{\sigma}}\|_{r+1}^{2}\,,$ (5.4)

and

$\big{|}\mathcal{S}^{K}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},{\boldsymbol{\tau}}_{h})\big{|}\leq
C_{\mathtt{s}}h^{r+1}\|{\boldsymbol{\sigma}}\|_{r+1}\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h},$
(5.5)

for all ${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}$

###### Proof.

The proof of vector versions of (5.4) and (5.5) can be found in [21, eqs. (33)
and (34)] and by following the similar arguments can be easily conclude them.
∎

Thanks to the projection error estimate stated in Lemma 3.3, we only analyze
the error functions defined by

${\boldsymbol{\theta}}_{h}:=\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}=\big{\\{}\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{0h},({\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n}_{e})-{\boldsymbol{\sigma}}_{bh}\mathbf{n}_{e})\otimes\mathbf{n}_{e}\big{\\}}{\quad\hbox{and}\quad}\theta_{h}:=\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}_{h}.$

As the primary step in convergence analysis, we derive error equation as
follows:

  * •

For any ${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}$, from the first equation of
Problem 2, we deduce

$\begin{array}[]{c}\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\theta_{h})\,=\,\big{[}\mathcal{A}_{h}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\boldsymbol{\mathcal{P}}_{h}\mathbf{u})\big{]}-\big{[}\mathcal{A}_{h}({\boldsymbol{\sigma}}_{h},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{u}_{h})\big{]}\\\\[8.61108pt]
\,=\,\mathcal{A}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\boldsymbol{\mathcal{P}}_{h}\mathbf{u})+\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})-\mathcal{G}({\boldsymbol{\tau}}_{h})\,.\end{array}$

We test the first equation of Problem 1, i.e.,
${\boldsymbol{\sigma}}^{\mathtt{d}}+(\mathbf{u}\otimes\mathbf{u})^{\mathtt{d}}-\nu\boldsymbol{\nabla}\mathbf{u}=\textbf{0}$
against ${\boldsymbol{\tau}}_{0h}$, where
${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}$, and add it on the right-hand side
of the above equation as well as use the L2-orthogonality property of
projections $\mathcal{P}\mathcal{P}_{0}$, $\boldsymbol{\mathcal{P}}_{h}$,
along with the definition of discrete weak divergence operator (cf. Definition
3.2) and Green formula, to obtain

$\begin{array}[]{c}\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\theta_{h})\,=\,\dfrac{1}{\nu}\displaystyle\int_{\Omega}(\mathcal{P}\mathcal{P}_{0}{\boldsymbol{\sigma}})^{\mathtt{d}}:{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}+\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\boldsymbol{\mathcal{P}}_{h}\mathbf{u})+\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})\\\\[8.61108pt]
\,-\,\Big{(}\mathcal{G}({\boldsymbol{\tau}}_{h})+\dfrac{1}{\nu}\displaystyle\int_{\Omega}{\boldsymbol{\sigma}}^{\mathtt{d}}:{\boldsymbol{\tau}}_{0h}^{\mathtt{d}}-\displaystyle\int_{\Omega}\boldsymbol{\nabla}\mathbf{u}:{\boldsymbol{\tau}}_{0h}+\dfrac{1}{\nu}\displaystyle\int_{\Omega}(\mathbf{u}\otimes\mathbf{u})^{\mathtt{d}}:{\boldsymbol{\tau}}_{0h}\Big{)}\\\\[10.76385pt]
\,=\,\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\Big{(}\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})-\mathcal{C}(\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})\Big{)}\\\\[8.61108pt]
\,+\,\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\langle({\boldsymbol{\tau}}_{0h}-{\boldsymbol{\tau}}_{bh})\mathbf{n},~{}\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\rangle_{\partial
K}\,.\end{array}$

  * •

An application of the second equation of Problems 1, that is,
$\operatorname{\textbf{div}}({\boldsymbol{\sigma}})+\mathbf{f}=0$ with test
function ${\mathbf{v}}_{h}\in\mathbf{Y}_{h}$ and Definition 3.2 and Green
formula, gives

$\mathcal{B}_{h}({\boldsymbol{\theta}}_{h},{\mathbf{v}}_{h})\,=\,\mathcal{B}_{h}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},{\mathbf{v}}_{h})-\mathcal{B}({\boldsymbol{\sigma}},{\mathbf{v}}_{h})=\sum\limits_{K\in\mathcal{K}_{h}}\langle{\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n})-{\boldsymbol{\sigma}}\mathbf{n},{\mathbf{v}}_{h}\rangle_{\partial
K}\,.$ (5.6)

Hence, the error equation reads as follows.

######  Problem 4 (error problem).

Find error functions ${\boldsymbol{\theta}}_{h}\in\mathbb{X}_{h}$ and
$\theta_{h}\in\mathbf{Y}_{h}$ such that

$\left\\{\begin{array}[]{rcll}\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\theta_{h})+\big{[}\mathcal{C}(\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})-\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})\big{]}&=&\mathcal{L}_{1}(({\boldsymbol{\sigma}},\mathbf{u});{\boldsymbol{\tau}}_{h})\,,\\\\[8.53581pt]
\mathcal{B}_{h}({\boldsymbol{\theta}}_{h},\mathbf{v}_{h})&=&\mathcal{L}_{2}({\boldsymbol{\sigma}},\mathbf{v}_{h})\,,\end{array}\right.$
(5.7)

in which

$\displaystyle\mathcal{L}_{1}(({\boldsymbol{\sigma}},\mathbf{u});{\boldsymbol{\tau}}_{h})$
$\displaystyle:=\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}{\boldsymbol{\tau}}_{h})+\sum\limits_{K\in\mathcal{K}_{h}}\langle({\boldsymbol{\tau}}_{0h}-{\boldsymbol{\tau}}_{bh})\mathbf{n},~{}\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\rangle_{\partial
K}\,,$ $\displaystyle\mathcal{L}_{2}({\boldsymbol{\sigma}},\mathbf{v}_{h})$
$\displaystyle:=\sum\limits_{K\in\mathcal{K}_{h}}\langle{\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n})-{\boldsymbol{\sigma}}\mathbf{n},\mathbf{v}_{h}\rangle_{\partial
K}\,.$

for all ${\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}$ and
$\mathbf{v}_{h}\in\mathbf{Y}_{h}$.

The following result states the estimate for error terms appeared in Problem
4.

######  Lemma 5.4.

Let ${\boldsymbol{\sigma}}\in\mathbb{X}\cap\mathbb{H}^{r+1}(\Omega)$,
$\mathbf{u}\in\mathbf{Y}\cap\mathbf{H}^{r+2}(\Omega)$ and $r\leq k$. Then,
there exist the constants $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ such that

$\begin{array}[]{rcll}\big{|}\mathcal{L}_{1}(({\boldsymbol{\sigma}},\mathbf{u});{\boldsymbol{\tau}}_{h})\big{|}&\leq&\,\mathcal{C}_{1}\,h^{r+1}\left(\|{\boldsymbol{\sigma}}\|_{r+1}+\|\mathbf{u}\|_{r+2}\right)\|{\boldsymbol{\tau}}_{h}\|_{\mathbb{H},h}&\quad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,,\\\\[8.61108pt]
\big{|}\mathcal{L}_{2}({\boldsymbol{\sigma}};\mathbf{v}_{h})\big{|}&\leq&\,\mathcal{C}_{2}\,h^{r+1}\|{\boldsymbol{\sigma}}\|_{r+1}\|\mathbf{v}_{h}\|_{1,h}&\quad\forall\,\mathbf{v}_{h}\in\mathbf{Y}_{h}\,.\end{array}$

###### Proof.

See [21, Lemma 4]. ∎

Now we are in position of establishing the main result of this section,
namely, the suboptimal rate of convergence for the weak Galerkin scheme
provided by Problem 3.

######  Theorem 5.5.

Assume that the data satisfy

$\dfrac{4}{\nu\alpha_{\Lambda}\alpha_{\Lambda,\mathtt{d}}}\left(\|\textbf{f}\|_{0,4/3}+c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}\right)\,\leq\,\dfrac{1}{2}\,,$
(5.8)

and let
$({\boldsymbol{\sigma}}_{h},{\mathbf{u}}_{h})\in\mathbb{X}_{h}\times\mathbf{Y}_{h}$
be the solution to Problem 3 and
$({\boldsymbol{\sigma}},{\mathbf{u}})\in\mathbb{X}\times\mathbf{Y}$ be the
solution of Problem 2 satisfying the regularity conditions
${\boldsymbol{\sigma}}\in\mathbb{H}^{k+1}(\Omega)$ and
$\mathbf{u}\in\mathbf{H}^{k+2}(\Omega)$. Then, for all $k\in\mathrm{N}_{0}$,
the following estimate holds

$\|\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|_{\mathbb{H},h}+\|\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}_{h}\|_{1,h}\,\leq\,\mathcal{C}_{\mathtt{sub}}\,h^{k+1}\left(\|{\boldsymbol{\sigma}}\|_{k+1}+\|\mathbf{u}\|_{k+2}\right)\,.$
(5.9)

###### Proof.

Bearing in mind the definition of $\Lambda_{h,{\mathbf{u}}_{h}}$ (cf. (4.14)),
utilizing Problem 4 (cf. (5.7)), and some simple computations we obtain

$\begin{array}[]{c}\Lambda_{h,\textbf{u}_{h}}[({\boldsymbol{\theta}}_{h},\theta_{h}),({\boldsymbol{\tau}}_{h},{\mathbf{v}}_{h})]\,=\,\big{[}\mathcal{C}(\mathbf{u}_{h};\theta_{h},{\boldsymbol{\tau}}_{h})-\mathcal{C}(\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})\big{]}\\\\[8.61108pt]
\,+\,\mathcal{L}_{1}(({\boldsymbol{\sigma}},\mathbf{u});{\boldsymbol{\tau}}_{h})+\mathcal{L}_{2}({\boldsymbol{\sigma}},\mathbf{v}_{h})\,.\end{array}$
(5.10)

Next, adding and subtracting some suitable terms and recalling that
$\theta_{h}=\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}_{h}$ the first
term on the right-hand side of the above equation can be rewritten as

$\begin{array}[]{c}\mathcal{C}(\mathbf{u}_{h};\theta_{h},{\boldsymbol{\tau}}_{h})-\mathcal{C}(\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},{\boldsymbol{\tau}}_{h})\,=\,\mathcal{C}(\mathbf{u}_{h};\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\mathbf{u}_{h}-\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})\\\\[8.61108pt]
\,=\,\big{[}\mathcal{C}(\mathbf{u}_{h};\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u};\mathbf{u},{\boldsymbol{\tau}}_{h})\big{]}-\mathcal{C}(\theta_{h};\textbf{u},{\boldsymbol{\tau}}_{h})\,.\end{array}$
(5.11)

Now, it suffices to substitute (5.11) back into (5.10) and employ the discrete
global inf-sup of $\Lambda_{h,\textbf{u}_{h}}$ (cf. Lemma 4.6, eq. (4.18)),
the continuity property of $\mathcal{C}$ (cf. (4.12)) and the estimats of
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ (cf. Lemma 5.4), to arrive at

$\begin{array}[]{c}\dfrac{\alpha_{\Lambda,\mathtt{d}}}{2}\|({\boldsymbol{\theta}}_{h},\theta_{h})\|_{h}\,\leq\,\dfrac{1}{\nu}\|\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}\|_{\mathbf{Y}}\left(c_{\mathtt{em}}\|\mathbf{u}_{h}\|_{1,h}+\|\mathbf{u}\|_{\mathbf{Y}}\right)+\dfrac{c_{\mathtt{em}}}{\nu}\|\mathbf{u}\|_{\mathbf{Y}}\|\theta_{h}\|_{1,h}\\\\[8.61108pt]
\,+\,\mathcal{C}_{3}\,h^{k+1}\left(\|{\boldsymbol{\sigma}}\|_{k+1}+\|\mathbf{u}\|_{k+2}\right)\,,\end{array}$
(5.12)

where $\mathcal{C}_{3}$ is a constant depending on $\mathcal{C}_{1}$,
$\mathcal{C}_{2}$.

In this way, using the fact that $\mathbf{u}\in\widehat{\mathbf{Y}}$ and
$\mathbf{u}_{h}\in\widehat{\mathbf{Y}}_{h}$, in particularly a priori bounds
provided by (4.19) and (4.21), and the inequality (5.12) we deduce that

$\begin{array}[]{c}\dfrac{\alpha_{\Lambda,\mathtt{d}}}{2}\|({\boldsymbol{\theta}}_{h},\theta_{h})\|_{h}\,\leq\,\left(\dfrac{4}{\nu\min\\{\alpha_{\Lambda},\alpha_{\Lambda,\mathtt{d}}\\}}\|\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}\|_{\mathbf{Y}}+\dfrac{2}{\alpha_{\Lambda}\nu}\|({\boldsymbol{\theta}}_{h},\theta_{h})\|_{h}\right)\left(\|\textbf{f}\|_{0,4/3}+c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}\right)\\\\[10.33327pt]
\,+\,\mathcal{C}_{3}\,h^{k+1}\left(\|{\boldsymbol{\sigma}}\|_{k+1}+\|\mathbf{u}\|_{k+2}\right)\,,\end{array}$

which together with assumption (5.8) implies (5.9) with considering

$C_{\mathtt{sub}}\,:=\,\dfrac{2}{\alpha_{\Lambda,\mathtt{d}}}\max\bigg{\\{}\dfrac{4}{\nu\min\\{\alpha_{\Lambda},\alpha_{\Lambda,\mathtt{d}}\\}}\left(\|\textbf{f}\|_{0,4/3}+c_{g}\|\textbf{g}\|_{1/2,\partial\Omega}\right),\mathcal{C}_{3}\bigg{\\}},$

and concludes the proof. ∎

###  5.2 Optimal convergence

To obtain an optimal order error estimate for the vector component
$\theta_{h}=\mathbf{u}_{h}-\boldsymbol{\mathcal{P}}_{h}\textbf{u}$ in the
usual $\mathbf{L}^{2}$-norm, we consider a dual problem that seeks
$\boldsymbol{\vartheta}$ and $\Phi$ satisfying

$\begin{cases}\hskip
2.84544pt\boldsymbol{\vartheta}^{\mathtt{d}}+(\mathbf{u}\otimes\Phi)^{\mathtt{d}}\,=\,\nu\boldsymbol{\nabla}\Phi&{\quad\hbox{in}\quad}\Omega\,,\\\\[4.30554pt]
\hskip
28.45274pt-\operatorname{\mathbf{div}}(\boldsymbol{\vartheta})\,=\,\theta_{h}&{\quad\hbox{in}\quad}\Omega\,,\\\\[4.30554pt]
\displaystyle\int_{\Omega}{\mathrm{tr}}(\boldsymbol{\vartheta})=0{\quad\hbox{and}\quad}\Phi\,=\,\mathbf{0}&{\quad\hbox{on}\quad}\partial\Omega\,,\end{cases}$
(5.13)

Assume that the usual $\mathrm{H}^{2}$-regularity is satisfied for the dual
problem; i.e., for any $\theta_{h}\in\mathbf{L}^{2}(\Omega)$, there exists a
unique solution
$(\boldsymbol{\vartheta},\Phi)\in\mathbb{H}^{1}(\Omega)\times\mathbf{H}^{2}(\Omega)$
such that

$\|\boldsymbol{\vartheta}\|_{1}+\|\Phi\|_{2}\,\leq\,C_{\mathtt{reg}}\|\theta_{h}\|_{0,\Omega}\,.$
(5.14)

Next, we give an error estimate for $\theta_{h}$ in the L2-norm.

######  Theorem 5.6.

Let the assumption of Theorem 5.5 be satisfied. Also, assume that
$(\boldsymbol{\vartheta},\Phi)\in\mathbb{H}^{1}(\Omega)\times\mathbf{H}^{2}(\Omega)$
be a solution of (5.13). Then, the following error estimation holds

$\|\boldsymbol{\mathcal{P}}_{h}\mathbf{u}-\mathbf{u}_{h}\|_{0}\leq\mathcal{C}_{\mathtt{opt}}\,h^{k+2}\left(\|\mathbf{u}\|_{k+2}+\|{\boldsymbol{\sigma}}\|_{k+1}\right)\,.$
(5.15)

###### Proof.

We divide the proof into three steps.
Step 1: discrete evolution equation for the error. First, for any
$\mathbf{v}\in\mathbf{H}^{1}(\Omega)$ using Definition 3.2 and Green formula,
we obtain the following result

$\begin{array}[]{c}\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\boldsymbol{\mathcal{P}}_{h}\mathbf{v})\,=\,\displaystyle\sum_{K\in\mathcal{K}_{h}}\big{[}-({\boldsymbol{\tau}}_{0h},\boldsymbol{\nabla}(\boldsymbol{\mathcal{P}}_{h}\mathbf{v}))_{K}+\langle{\boldsymbol{\tau}}_{bh}\mathbf{n},\boldsymbol{\mathcal{P}}_{h}\mathbf{v}\rangle_{\partial
K}\big{]}\\\\[8.61108pt]
\,=\,\displaystyle\sum_{K\in\mathcal{K}_{h}}\big{[}(\operatorname{\mathbf{div}}({\boldsymbol{\tau}}_{0h}),\boldsymbol{\mathcal{P}}_{h}\mathbf{v})_{K}+\langle({\boldsymbol{\tau}}_{bh}-{\boldsymbol{\tau}}_{0h})\mathbf{n},\boldsymbol{\mathcal{P}}_{h}\mathbf{v}\rangle_{\partial
K}\big{]}\\\\[8.61108pt]
\,=\,\displaystyle\sum_{K\in\mathcal{K}_{h}}\big{[}(\operatorname{\mathbf{div}}({\boldsymbol{\tau}}_{0h}),\mathbf{v})_{K}+\langle({\boldsymbol{\tau}}_{bh}-{\boldsymbol{\tau}}_{0h})\mathbf{n},\boldsymbol{\mathcal{P}}_{h}\mathbf{v}\rangle_{\partial
K}\big{]}\\\\[8.61108pt]
\,=\,\displaystyle\sum_{K\in\mathcal{K}_{h}}\big{[}-({\boldsymbol{\tau}}_{0h},\boldsymbol{\nabla}\mathbf{v})_{K}+\langle{\boldsymbol{\tau}}_{0h}\mathbf{n},\mathbf{v}\rangle_{\partial
K}+\langle({\boldsymbol{\tau}}_{bh}-{\boldsymbol{\tau}}_{0h})\mathbf{n},\boldsymbol{\mathcal{P}}_{h}\mathbf{v}\rangle_{\partial
K}\big{]}\\\\[8.61108pt]
\,=\,\displaystyle\sum_{K\in\mathcal{K}_{h}}\big{[}-({\boldsymbol{\tau}}_{0h},\boldsymbol{\nabla}\mathbf{v})_{K}+\langle({\boldsymbol{\tau}}_{0h}-{\boldsymbol{\tau}}_{bh})\mathbf{n},\mathbf{v}-\boldsymbol{\mathcal{P}}_{h}\mathbf{v}\rangle_{\partial
K}+\langle{\boldsymbol{\tau}}_{bh}\mathbf{n},\mathbf{v}\rangle_{\partial
K}\big{]}\,.\end{array}$ (5.16)

Now, the testing of the Eq. (5.13) against
$({\boldsymbol{\theta}}_{h},\theta_{h})$ yields

$\displaystyle\mathcal{A}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})+\mathcal{C}(\mathbf{u};\Phi,{\boldsymbol{\theta}}_{h})$
$\displaystyle\,=\,\displaystyle\int_{\Omega}\boldsymbol{\nabla}\Phi:{\boldsymbol{\theta}}_{0h}\,,$
(5.17a) $\displaystyle-\mathcal{B}(\boldsymbol{\vartheta},\theta_{h})$
$\displaystyle\,=\,(\theta_{h},\theta_{h})_{0,\Omega}\,.$ (5.17b)

By employing (5.16) and the fact that $\Phi=\textbf{0}$ on $\partial\Omega$,
the term on the right-hand side of (5.17a) can be rewritten as:

$\displaystyle\int_{\Omega}\boldsymbol{\nabla}\Phi:{\boldsymbol{\theta}}_{0h}\,=\,-\mathcal{B}_{h}({\boldsymbol{\theta}}_{h},\boldsymbol{\mathcal{P}}_{h}\Phi)+\sum_{K\in\mathcal{K}_{h}}\langle({\boldsymbol{\theta}}_{0h}-{\boldsymbol{\theta}}_{bh})\mathbf{n},\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\rangle_{\partial
K}\,,$

from which, replacing the first term on the right hand by the second row of
(5.7), and using the fact that $\Phi\in\mathbf{H}^{1}(\Omega)$ and
$\Phi=\textbf{0}$ on $\partial\Omega$ gives

$\begin{array}[]{c}\displaystyle\int_{\Omega}\boldsymbol{\nabla}\Phi:{\boldsymbol{\theta}}_{0h}\,=\,\mathcal{L}_{2}({\boldsymbol{\sigma}},\boldsymbol{\mathcal{P}}_{h}\Phi)+\displaystyle\sum_{K\in\mathcal{K}_{h}}\langle({\boldsymbol{\theta}}_{0h}-{\boldsymbol{\theta}}_{bh})\mathbf{n},\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\rangle_{\partial
K}\\\\[8.61108pt]
\,=\,\underbrace{\sum\limits_{K\in\mathcal{K}_{h}}\langle{\boldsymbol{\sigma}}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n}),\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\rangle_{\partial
K}}_{=:\mathrm{I}_{1}}+\underbrace{\sum_{K\in\mathcal{K}_{h}}\langle({\boldsymbol{\theta}}_{0h}-{\boldsymbol{\theta}}_{bh})\mathbf{n},\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\rangle_{\partial
K}}_{=:\mathrm{I}_{2}}\,.\end{array}$ (5.18)

On the other hand, by following the similar arguments of (5.6) the term on the
left-hand side of (5.17b) can be rewritten by

$\displaystyle\mathcal{B}(\boldsymbol{\vartheta},\theta_{h})$
$\displaystyle=\mathcal{B}_{h}(\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta},\theta_{h})+\sum\limits_{K\in\mathcal{K}_{h}}\langle\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n}),\theta_{h}\rangle_{\partial
K}$
$\displaystyle=\mathcal{L}_{1}(({\boldsymbol{\sigma}},\mathbf{u});\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\big{[}\mathcal{C}(\mathbf{u};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{]}$
$\displaystyle~{}~{}~{}+\sum\limits_{K\in\mathcal{K}_{h}}\langle\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n}),\theta_{h}\rangle_{\partial
K}$
$\displaystyle=\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\big{[}\mathcal{C}(\mathbf{u};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{]}$
$\displaystyle~{}~{}~{}+\sum\limits_{K\in\mathcal{K}_{h}}\langle(\mathcal{P}\mathcal{P}_{0}\boldsymbol{\vartheta}-\mathcal{P}\mathcal{P}_{b}\boldsymbol{\vartheta})\mathbf{n},~{}\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\rangle_{\partial
K}+\sum\limits_{K\in\mathcal{K}_{h}}\langle\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n}),\theta_{h}\rangle_{\partial
K}$ $\displaystyle=:\sum_{i=3}^{7}\mathrm{I}_{i}\,,$

where the first row of (5.7) and the definition of $\mathcal{L}_{1}$ were
applied in the second and fourth lines, respectively.
Step 2: bounding the error terms $\mathrm{I}_{1}$-$\mathrm{I}_{7}$. For the
term $\mathrm{I}_{1}$, applying the Cauchy-Schwarz inequality, the bound (5.4)
stated in Lemma 5.3 and trace inequality (cf. Lemma 5.1) we estimate

$\begin{array}[]{c}\big{|}\mathrm{I}_{1}\big{|}\,\leq\,\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}h_{K}\|{\boldsymbol{\sigma}}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}({\boldsymbol{\sigma}}\mathbf{n})\|_{0,\partial
K}^{2}\right)^{1/2}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}h_{K}^{-1}\|\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0,\partial
K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,Ch^{k+1}\|{\boldsymbol{\sigma}}\|_{k+1}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\big{[}h_{K}^{-2}\|\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0,K}^{2}+\|\nabla(\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi)\|_{0,K}^{2}\big{]}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,Ch^{k+2}\|{\boldsymbol{\sigma}}\|_{k+1}\|\Phi\|_{2}\,.\end{array}$
(5.19)

As for the term $\mathrm{I}_{2}$, the Cauchy-Schwarz and trace inequalities
and the definition of discrete norm $\|\cdot\|_{\mathbb{H},h}$ (cf. first
paragraph of Sec. 4.1) implies that,

$\begin{array}[]{c}\big{|}\mathrm{I}_{2}\big{|}\,\leq\,\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\|({\boldsymbol{\theta}}_{0h}-{\boldsymbol{\theta}}_{bh})\mathbf{n}\|_{0,\partial
K}^{2}\right)^{1/2}\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}^{-1}\|\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0,\partial
K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,\mathcal{C}_{\mathtt{tr}}\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\big{[}h_{K}^{-2}\|\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0,K}^{2}+\|\nabla(\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi)\|_{0,K}^{2}\big{]}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,\mathcal{C}_{\mathtt{tr}}\,h\,\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\|\Phi\|_{2}\,.\end{array}$
(5.20)

For term $\mathrm{I}_{3}$, we use the defintion of $\mathcal{S}(\cdot,\cdot)$
given by (3.10), the continuity and orthogonality properties of operator
${\boldsymbol{\mathcal{P}}}_{b}^{K}$ and trace inequality, to get

$\begin{array}[]{c}\big{|}\mathrm{I}_{3}\big{|}\,=\,\big{|}\mathcal{S}(\mathcal{P}\mathcal{P}_{h}{\boldsymbol{\sigma}},~{}\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{|}\\\\[8.61108pt]
\,=\,\bigg{|}\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\left\langle{\boldsymbol{\mathcal{P}}}_{b}^{K}((\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}})\mathbf{n}-{\boldsymbol{\sigma}}\mathbf{n}),~{}{\boldsymbol{\mathcal{P}}}_{b}^{K}((\mathcal{P}\mathcal{P}_{0}^{K}\boldsymbol{\vartheta})\mathbf{n}-\boldsymbol{\vartheta}\mathbf{n})\right\rangle_{0,\partial
K}\bigg{|}\\\\[8.61108pt]
\,\leq\,\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\|{\boldsymbol{\mathcal{P}}}_{b}^{K}((\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}})\mathbf{n}-{\boldsymbol{\sigma}}\mathbf{n})\|_{0,\partial
K}^{2}\right)^{1/2}\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\|{\boldsymbol{\mathcal{P}}}_{b}^{K}((\mathcal{P}\mathcal{P}_{0}^{K}\boldsymbol{\vartheta})\mathbf{n}-\boldsymbol{\vartheta}\mathbf{n})\|_{0,\partial
K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\|(\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}})\mathbf{n}-{\boldsymbol{\sigma}}\mathbf{n}\|_{0,\partial
K}^{2}\right)^{1/2}\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}h_{K}\|(\mathcal{P}\mathcal{P}_{0}^{K}\boldsymbol{\vartheta})\mathbf{n}-\boldsymbol{\vartheta}\mathbf{n}\|_{0,\partial
K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}\|\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}\|_{0,K}^{2}+h_{K}^{2}\|\boldsymbol{\nabla}(\mathcal{P}\mathcal{P}_{0}^{K}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}})\|_{0,K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\times\,\left(\displaystyle\sum_{K\in\mathcal{K}_{h}}\|\mathcal{P}\mathcal{P}_{0}^{K}\boldsymbol{\vartheta}-\boldsymbol{\vartheta}\|_{0,K}^{2}+h_{K}^{2}\|\boldsymbol{\nabla}(\mathcal{P}\mathcal{P}_{0}^{K}\boldsymbol{\vartheta}-\boldsymbol{\vartheta})\|_{0,K}^{2}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,Ch^{k+2}\|{\boldsymbol{\sigma}}\|_{k+1}\|\boldsymbol{\vartheta}\|_{1}\,,\end{array}$
(5.21)

where the approximation property of the projector
$\mathcal{P}\mathcal{P}_{0}^{K}$ (cf. Lemma 3.3) was used in the last step.
Next, in order to estimate $\mathrm{I}_{4}$, we add zero in the form
$\pm\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})$, to
find that

$\mathrm{I}_{4}\,=\,\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\mathcal{P}\mathcal{P}_{h}^{K}\boldsymbol{\vartheta}_{h}-\boldsymbol{\vartheta})\,+\,\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})\,.$
(5.22)

The first term on the right-hand side of the above equation can be estimated
using the continuity of $\mathcal{A}_{h}$ given in Lemma 4.2 and the
approximation property of the projector $\mathcal{P}\mathcal{P}_{0}^{K}$ as

$\big{|}\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\mathcal{P}\mathcal{P}_{h}^{K}\boldsymbol{\vartheta}_{h}-\boldsymbol{\vartheta})\big{|}\,\leq\,c_{\mathcal{A}}\,\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\,\|\mathcal{P}\mathcal{P}_{h}^{K}\boldsymbol{\vartheta}_{h}-\boldsymbol{\vartheta}\|_{\mathbb{H},h}\,\leq\,Ch\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\|\boldsymbol{\vartheta}\|_{1}\,.$
(5.23)

Also, an application of the consistency property of $\mathcal{A}_{h}$ (cf.
Lemma 3.3) and testing the first row of (5.13) with
${\boldsymbol{\theta}}_{h}$, yields

$\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})\,=\,\mathcal{A}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})\,=\,-\displaystyle\int_{\Omega}\boldsymbol{\nabla}\Phi:{\boldsymbol{\theta}}_{0h}-\mathcal{C}(\textbf{u};\Phi,{\boldsymbol{\theta}}_{h})\,=\,-\mathrm{I}_{1}-\mathrm{I}_{2}-\mathcal{C}(\textbf{u};\Phi,{\boldsymbol{\theta}}_{h})\,,$
(5.24)

where the the equivalent expression of
$\int_{\Omega}\boldsymbol{\nabla}\Phi:{\boldsymbol{\theta}}_{0h}$ given by
(5.18) was used in the last step. Then, thanks to (5.19) and (5.20) it
suffices to estimate the last term on the right-hand side. For this purpose,
we rewrite the associated term by adding and subtracting some suitable term as

$\mathcal{C}(\textbf{u};\Phi,{\boldsymbol{\theta}}_{h})\,=\,\mathcal{C}(\textbf{u};\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi,{\boldsymbol{\theta}}_{h})\,+\,\mathcal{C}(\textbf{u};\boldsymbol{\mathcal{P}}_{h}\Phi,{\boldsymbol{\theta}}_{h})\,.$
(5.25)

To determine upper bounds for the right-hand side terms, we use Hölder
inequality, the approximation and the continuity properties of the projector
$\boldsymbol{\mathcal{P}}_{h}$, Gagliardo–Nirenberg inequality and inverse
inequality. This gives

$\big{|}\mathcal{C}(\textbf{u};\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi,{\boldsymbol{\theta}}_{h})\big{|}\,\leq\,\dfrac{1}{\nu}\|\textbf{u}\|_{\mathbf{Y}}\|\Phi-\boldsymbol{\mathcal{P}}_{h}\Phi\|_{\mathbf{Y}}\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\,\leq\,\dfrac{1}{\nu}\|\textbf{u}\|_{\mathbf{Y}}\,h\,\|\Phi\|_{2}\,\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\,,$

and

$\displaystyle\big{|}\mathcal{C}(\textbf{u};\boldsymbol{\mathcal{P}}_{h}\Phi,{\boldsymbol{\theta}}_{h})\big{|}$
$\displaystyle\,\leq\,\dfrac{1}{\nu}\|\textbf{u}\|_{\mathbf{Y}}\|\boldsymbol{\mathcal{P}}_{h}\Phi\|_{\mathbf{Y}}\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}$
$\displaystyle\,\leq\,\dfrac{1}{\nu}\|\textbf{u}\|_{\mathbf{Y}}\left(C_{\mathtt{GN}}\|\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0}^{1/2}\|\nabla(\boldsymbol{\mathcal{P}}_{h}\Phi)\|_{0}^{1/2}\right)\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}$
$\displaystyle\,\leq\,\dfrac{1}{\nu}\|\textbf{u}\|_{\mathbf{Y}}\mathcal{C}_{\mathtt{inv}}h^{1/2}\left(C_{\mathtt{GN}}\|\boldsymbol{\mathcal{P}}_{h}\Phi\|_{0,4}^{1/2}\|\nabla(\boldsymbol{\mathcal{P}}_{h}\Phi)\|_{0,4}^{1/2}\right)\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}$
$\displaystyle\,\leq\,Ch\|\textbf{u}\|_{\mathbf{Y}}\|\Phi\|_{2}\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\,.$

Combining the above bounds with (5.25) and (5.24), along with (5.19), (5.20)
yields

$\big{|}\mathcal{A}_{h}({\boldsymbol{\theta}}_{h},\boldsymbol{\vartheta})\big{|}\,\leq\,C\left(h^{k+2}\|{\boldsymbol{\sigma}}\|_{k+1}+h\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\right)\|\Phi\|_{2}\,,$

which, together with estimate (5.23) gives

$\big{|}\mathrm{I}_{4}\big{|}\,\leq\,C\left(h^{k+2}\|{\boldsymbol{\sigma}}\|_{k+1}+h\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\right)\|\Phi\|_{2}+Ch\|{\boldsymbol{\theta}}_{h}\|_{\mathbb{H},h}\|\boldsymbol{\vartheta}\|_{1}\,.$
(5.26)

In what follow, we focus on the estimate of $\mathrm{I}_{5}$. To this end,
using similar arguments as in the proof of Theorem 5.5 (cf. equation (5.11)),
we obtain

$\begin{array}[]{c}\mathrm{I}_{5}\,=\,\mathcal{C}(\mathbf{u};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})-\mathcal{C}(\mathbf{u}_{h};\mathbf{u}_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\\\\[8.61108pt]
\,=\,\mathcal{C}(\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})+\mathcal{C}(\mathbf{u}_{h};\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})+\mathcal{C}(\theta_{h};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})+\mathcal{C}(\mathbf{u}_{h};\theta_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\,,\end{array}$
(5.27)

and therefore, each of the terms above are estimated as follows:

$\displaystyle\big{|}\mathcal{C}(\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{|}$
$\displaystyle\,\leq\,\|\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\|_{0}\|\mathbf{u}\|_{0,4}\|\mathcal{P}\mathcal{P}_{0}\boldsymbol{\vartheta}\|_{0,4}\leq
Ch^{k+2}\|\mathbf{u}\|_{k+2}\|\mathbf{u}\|_{\mathbf{Y}}\|\boldsymbol{\vartheta}\|_{1}$
$\displaystyle\big{|}\mathcal{C}(\mathbf{u}_{h};\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{|}$
$\displaystyle\,\leq\,Ch^{k+2}\|\mathbf{u}\|_{k+2}\|\mathbf{u}_{h}\|_{\mathbf{Y}}\|\boldsymbol{\vartheta}\|_{1}$
$\displaystyle\big{|}\mathcal{C}(\theta_{h};\mathbf{u},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{|}$
$\displaystyle\,\leq\,\|\theta_{h}\|_{0,4}\|\mathbf{u}\|_{0,4}\|\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta}\|_{0}$
$\displaystyle\,\leq\,\|\theta_{h}\|_{0,4}\|\mathbf{u}\|_{\mathbf{Y}}\mathcal{C}_{\mathtt{inv}}h^{1/2}\|\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta}\|_{0,4}$
$\displaystyle\,\leq\,(C_{\mathtt{GN}}\|\theta_{h}\|_{0}^{1/2}\|\nabla\theta_{h}\|_{0}^{1/2})\|\mathbf{u}\|_{\mathbf{Y}}\mathcal{C}_{\mathtt{inv}}h^{1/2}\left(C_{\mathtt{GN}}\|\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta}\|_{0}^{1/2}\|\nabla(\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\|_{0}^{1/2}\right)$
$\displaystyle\,\leq\,h^{\frac{1}{4}}(C_{\mathtt{GN}}\|\theta_{h}\|_{0,4}^{1/2}\|\nabla\theta_{h}\|_{0}^{1/2})\|\mathbf{u}\|_{\mathbf{Y}}\mathcal{C}_{\mathtt{inv}}^{3}h^{\frac{3}{4}}\left(C_{\mathtt{GN}}\|\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta}\|_{0,4}^{1/2}\|\nabla(\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\|_{0}^{1/2}\right)$
$\displaystyle\,\leq\,Ch\|\nabla\theta_{h}\|_{0}\|\mathbf{u}\|_{\mathbf{Y}}\|\nabla(\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\|_{0}\leq
Ch\|\theta_{h}\|_{h}\|\mathbf{u}\|_{\mathbf{Y}}\|\boldsymbol{\vartheta}\|_{1}$
$\displaystyle\big{|}\mathcal{C}(\mathbf{u}_{h};\theta_{h},\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta})\big{|}$
$\displaystyle\,\leq\,Ch\|\theta_{h}\|_{h}\|\mathbf{u}_{h}\|_{\mathbf{Y}}\|\boldsymbol{\vartheta}\|_{1}$

which leads to

$\big{|}\mathrm{I}_{5}\big{|}\,\leq\,C\left(h^{k+2}\|\mathbf{u}\|_{k+2}+h\|\theta_{h}\|_{h}\right)(\|\mathbf{u}\|_{\mathbf{Y}}+\|\mathbf{u}_{h}\|_{\mathbf{Y}})\|\boldsymbol{\vartheta}\|_{1}\,.$
(5.28)

On the other hand, by arguments similar to those used in the estimation term
$\mathrm{I}_{1}$, we derive

$\begin{array}[]{c}\big{|}\mathrm{I}_{6}\big{|}\,=\,\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\bigg{|}\langle(\mathcal{P}\mathcal{P}_{0}\boldsymbol{\vartheta}-\mathcal{P}\mathcal{P}_{b}\boldsymbol{\vartheta})\mathbf{n},~{}\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\rangle_{\partial
K}\bigg{|}\\\\[4.30554pt]
\,=\,\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\bigg{|}\langle(\mathcal{P}\mathcal{P}_{0}\boldsymbol{\vartheta}-\boldsymbol{\vartheta})\mathbf{n}-(\mathcal{P}\mathcal{P}_{b}\boldsymbol{\vartheta}-\boldsymbol{\vartheta})\mathbf{n},~{}\mathbf{u}-\boldsymbol{\mathcal{P}}_{h}\mathbf{u}\rangle_{\partial
K}\bigg{|}\\\\[4.30554pt]
\,\leq\,\|\mathcal{P}\mathcal{P}_{h}\boldsymbol{\vartheta}-\boldsymbol{\vartheta}\|_{\mathbb{H},h}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\big{[}h_{K}^{-2}\|\textbf{u}-\boldsymbol{\mathcal{P}}_{h}\textbf{u}\|_{0,K}^{2}+\|\nabla(\textbf{u}-\boldsymbol{\mathcal{P}}_{h}\textbf{u})\|_{0,K}^{2}\big{]}\right)^{1/2}\\\\[4.30554pt]
\,\leq\,Ch^{k+2}\|\textbf{u}\|_{k+2}\|\boldsymbol{\vartheta}\|_{1}\,.\end{array}$

Similarly, using (5.4) we have

$\begin{array}[]{c}\big{|}\mathrm{I}_{7}\big{|}\,=\,\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\bigg{|}\langle\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n}),\theta_{h}\rangle_{\partial
K}\bigg{|}=\bigg{|}\langle\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n}),\theta_{h}-\bar{\theta}_{h}\rangle_{\partial
K}\bigg{|}\\\\[4.30554pt]
\,\leq\,\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}h_{K}\|\boldsymbol{\vartheta}\mathbf{n}-{\boldsymbol{\mathcal{P}}}_{b}^{K}(\boldsymbol{\vartheta}\mathbf{n})\|_{0,\partial
K}^{2}\right)^{1/2}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}h_{K}^{-1}\|\theta_{h}-\bar{\theta}_{h}\|_{0,\partial
K}^{2}\right)^{1/2}\\\\[4.30554pt]
\,\leq\,Ch\|\boldsymbol{\vartheta}\|_{1}\left(\displaystyle\sum\limits_{K\in\mathcal{K}_{h}}\big{[}h_{K}^{-2}\|\theta_{h}-\bar{\theta}_{h}\|_{0,K}^{2}+\|\nabla\theta_{h}\|_{0,K}^{2}\big{]}\right)^{1/2}\\\\[8.61108pt]
\,\leq\,Ch\|\theta_{h}\|_{1,h}\|\boldsymbol{\vartheta}\|_{1}\,.\end{array}$

Step 3: error estimate. We now insert the bounds on
$\mathrm{I}_{1}$-$\mathrm{I}_{7}$ in (5.17b), yielding

$\|\theta_{h}\|_{0}\,\leq\,\Big{(}Ch^{k+2}+h(\|\theta_{h}\|_{1,h}+\|{\boldsymbol{\theta}}\|_{\mathbb{H},h})\Big{)}(\|\boldsymbol{\vartheta}\|_{1}+\|\Phi\|_{2})\,.$
(5.29)

The sought result follows from Theorem 5.5 and employing regularity (5.14). ∎

We end this section by establishing the error estimate for the pressure. By
proceeding as in [17, Theorem 5.5, eqs. (5.38) and (5.39)] (see also [18, eq.
5.14]), we deduce the existence of a positive constant $C$, independent of
$h$, such that

$\|p-p_{h}\|_{0,\Omega}\,\leq\,C\,\Big{\\{}\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|_{0,\Omega}\,+\,\|{\mathbf{u}}-{\mathbf{u}}_{h}\|_{0,4;\Omega}\Big{\\}}\,.$
(5.30)

Therefore, thanks to Theorems 5.5 and 5.6 we get

$\|p-p_{h}\|_{0}\leq
C\,h^{k+1}\left(\|{\boldsymbol{\sigma}}\|_{k+1}+\|\mathbf{u}\|_{k+2}\right).$
(5.31)

##  6 Numerical Results

The purpose of this section is to illustrate the efficiency of the WG
pseudostress-based mixed-FEM for solving the Brinkman problem. Like in Ref.
[18], we use a real Lagrange multiplier to impose the zero integral mean
condition for ${\boldsymbol{\sigma}}_{h}$ of the discrete scheme. As a result,
Poblem 3 is rewritten as follow: find
$(({\boldsymbol{\sigma}}_{h},\mathbf{u}_{h}),\lambda)\in\mathbb{X}_{0,h}\times\mathbf{Y}_{h}\times\mathrm{R}$
such that

$\left\\{\begin{array}[]{rcll}\mathcal{A}_{h}({\boldsymbol{\sigma}}_{h},{\boldsymbol{\tau}}_{h})+\mathcal{C}(\mathbf{u}_{h},\mathbf{u}_{h};{\boldsymbol{\tau}}_{h})+\mathcal{B}_{h}({\boldsymbol{\tau}}_{h},\mathbf{u}_{h})+\lambda\int_{\Omega}\operatorname{tr}({\boldsymbol{\tau}}_{h})&=&\mathcal{G}({\boldsymbol{\tau}}_{h}),&\quad\forall\,{\boldsymbol{\tau}}_{h}\in\mathbb{X}_{h}\,,\\\\[8.53581pt]
\mathcal{B}_{h}({\boldsymbol{\sigma}}_{h},{\mathbf{v}}_{h})&=&\mathcal{F}({\mathbf{v}}_{h}),&\quad\forall{\mathbf{v}}_{h}\in\mathbf{Y}_{h}\,,\\\\[8.53581pt]
\xi\displaystyle\int_{\Omega}\operatorname{tr}({\boldsymbol{\sigma}}_{h})&=&0,&\quad\forall\,\xi\in\mathrm{R}\,.\end{array}\right.$
(6.1)

In addition, the individual errors associated to the main unknowns and the
postprocessed pressure are denoted and defined, as usual, by

$\begin{array}[]{c}\displaystyle{\tt
e}(\boldsymbol{\sigma})\,:=\,\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|_{{\mathbf{div}}_{4/3};\Omega}\,,\quad{\tt
e}(\mathbf{u})\,:=\,\|\mathbf{u}-\mathbf{u}_{h}\|_{0,4;\Omega}\,,{\quad\hbox{and}\quad}{\tt
e}(p)\,:=\,\|p-\widehat{p}_{h}\|_{0,\Omega}\,.\end{array}$

In turn, for all $\star\in\\{{\boldsymbol{\sigma}},{\mathbf{u}},p\\}$, we let
${\tt r}(\star)~{}:=~{}\dfrac{\text{log}({\tt e}(\star)/{\tt
e}^{\prime}(\star))}{\text{log}(h/h^{\prime})}$ be the experimental rates of
convergence, where $h$ and $h^{{}^{\prime}}$ denote two consecutive mesh sizes
with errors ${\tt e}(\star)$ and ${\tt e}^{{}^{\prime}}(\star)$, respectively.

In addition, a fixed point strategy with a fixed tolerance $\mathtt{Tol}=1$e-6
is utilized for the solving of the nonlinear equation (6.1). To that end, we
begin with vector all zeros as an initial guess and stop iterations when the
velocity’s error between two concessive iterations is adequately small,

$\|{\boldsymbol{\sigma}}^{m}-{\boldsymbol{\sigma}}_{h}^{m}\|_{0}+\|\mathbf{u}^{m}-\mathbf{u}_{h}^{m}\|_{\mathbf{Y}}\leq\mathtt{Tol}\,.$

In the first example, we mainly verify the accuracy of the constructed
numerical scheme. Examples 2 and 3 is utilized to evaluate the effectiveness
of the discrete scheme by simulation of practical problems for which no
analytical solutions. We performed our computations using the MATLAB 2020 b
software on an Intel Core i7 machine with 32 GB of memory.

Figure 1: Example 1, samples of the kind of meshes utilized.

###  6.1 Example 1: Accuracy assessment

We turn first to the numerical verification of the rates of convergence
anticipated by Theorems 5.5 and 5.6. To this end, we consider parameter
$\nu=0.1$ and design the exact solution as follows:

$\begin{array}[]{l}\mathbf{u}(x_{1},x_{2})\,=\,\begin{pmatrix}x_{1}^{2}\exp(-x_{1})(1+x_{2})\left(2\sin(1+x_{2})+(1+x_{2})\cos(1+x_{2})\right)\\\\[5.69054pt]
x_{1}(x_{1}-2)\exp(-x_{1})(1+x_{2})^{2}\sin(1+x_{2})\end{pmatrix}\,,\end{array}$
$p(x_{1},x_{2})\,=\,\sin(2\pi x_{1})\sin(2\pi x_{2})\,,$

for all $(x_{1},x_{2})^{\mathtt{t}}\in\Omega=(0,1)^{2}$. The model problem is
then complemented with the appropriate Dirichlet boundary condition. Using the
weak Galerkin spaces given in Sec. 3 with polynomial degree k$=0$, $1$, we
solve Problem 3 and obtain the approximated stress on a sequence of four
succesively refined polygonal meshes mades of hexagons, non-convex and
triangular elements (see Fig. 1). In addition, the discrete velocity field and
the discrete pressure are computed using post-processing approach stated in
Sec. 5. At each refinement level we compute errors between approximate and
smooth exact solutions in L2-norm. The results of this convergence study are
collected in Tables 1-3. One can see that the rate of convergence of
individual stress and pressure variables is $O(h^{k})$, whereas it is
$O(h^{k+1})$ for the velocity, which both are in agreement with the
theoretical analysis stated in Theorems 5.5 and 5.6. On the other hand, in
order to illustrate the accurateness of the discrete scheme, in Fig. 2 we
display components of the approximate velocity, stress, pressure on polygonal
mesh with $h=3.030$e-2 and k$=0$.

---
Figure 2: Example 1, snapshots of the numerical stress components (first row,
left to right), the velocity components and pressure (second row, left to
right), computed with $k=0$ in the mesh made of hexagons with $h=3.030$e-2.
Table 1: Example 1, history of convergence using hexagon

$k$ $h$ $\mathtt{e}(\boldsymbol{\sigma})$ ${\tt r}(\boldsymbol{\sigma})$
$\mathtt{e}(\mathbf{u})$ ${\tt r}(\mathbf{u})$ $\mathtt{e}(p)$ ${\tt r}(p)$
2.000e-01 1.27585e+00 - 1.63990e-01 - 3.92748e-01 - 1.111e-01 6.88337e-01
1.04985e+00 5.33247e-02 1.91124e+00 2.07409e-01 1.08623e+00 0 5.882e-02
3.14151e-01 1.23336e+00 1.58302e-02 1.90960e+00 9.47219e-02 1.23233e+00
3.030e-02 1.17182e-01 1.48674e+00 4.37012e-03 1.94051e+00 3.57120e-02
1.47063e+00 2.000e-01 3.34158e-01 - 3.63525e-02 - 9.11384e-02 - 1.111e-01
7.65257e-02 2.50769e+00 5.47756e-03 3.21988e+00 2.16022e-02 2.44916e+00 1
5.882e-02 1.64299e-02 2.41911e+00 7.67223e-04 3.09068e+00 5.13068e-03
2.26035e+00 3.030e-02 3.95366e-03 2.14755e+00 1.19847e-04 2.79901e+00
1.30515e-03 2.06382e+00

Table 2: Example 1, history of convergence using non-convex

$k$ $h$ $\mathtt{e}(\boldsymbol{\sigma})$ ${\tt r}(\boldsymbol{\sigma})$
$\mathtt{e}(\mathbf{u})$ ${\tt r}(\mathbf{u})$ $\mathtt{e}(p)$ ${\tt r}(p)$
2.500e-01 1.92411e+00 - 2.16399e-01 - 5.24963e-01 - 1.250e-01 9.61276e-01
1.00117e+00 6.67542e-02 1.69676e+00 2.69741e-01 9.60642e-01 0 6.250e-02
3.68609e-01 1.38286e+00 1.91857e-02 1.79883e+00 1.04967e-01 1.36163e+00
3.125e-02 1.21234e-01 1.60429e+00 5.08183e-03 1.91661e+00 3.54124e-02
1.56761e+00 2.500e-01 6.42714e-01 - 4.33653e-02 - 1.82378e-01 - 1.250e-01
2.05453e-01 1.64537e+00 6.90082e-03 2.65170e+00 6.23466e-02 1.54855e+00 1
6.250e-02 5.59423e-02 1.87679e+00 1.35203e-03 2.35164e+00 1.76339e-02
1.82196e+00 3.125e-02 1.48719e-02 1.91135e+00 3.12928e-04 2.11122e+00
4.77578e-03 1.88455e+00

Table 3: Example 1, history of convergence using triangular

$k$ $h$ $\mathtt{e}(\boldsymbol{\sigma})$ ${\tt r}(\boldsymbol{\sigma})$
$\mathtt{e}(\mathbf{u})$ ${\tt r}(\mathbf{u})$ $\mathtt{e}(p)$ ${\tt r}(p)$
1.768e-01 1.33737e+00 - 1.21833e-01 - 2.64186e-01 - 8.839e-02 7.00927e-01
9.32063e-01 3.13360e-02 1.95902e+00 1.22818e-01 1.10503e+00 0 4.419e-02
3.39792e-01 1.04461e+00 7.98954e-03 1.97164e+00 5.00518e-02 1.29503e+00
2.210e-02 1.63696e-01 1.05363e+00 2.01389e-03 1.98813e+00 1.99353e-02
1.32810e+00 1.768e-01 7.89707e-01 - 2.86149e-02 - 1.69212e-01 - 8.839e-02
2.50022e-01 1.65926e+00 5.94594e-03 2.26679e+00 5.68132e-02 1.57454e+00 1
4.419e-02 7.26955e-02 1.78212e+00 1.38897e-03 2.09789e+00 1.72263e-02
1.72162e+00 2.210e-02 2.02924e-02 1.84093e+00 3.40546e-04 2.02810e+00
4.95787e-03 1.79682e+00

###  6.2 Example 2 [24]: Lid-Driven Cavity Problem

The next example is chosen to illustrate the performance of the proposed
method for modeling the lid-driven cavity flow in the square domain
$\Omega=(0,1)^{2}$ with different values of $\nu$. For boundary conditions, we
set the inflow $\textbf{g}=(1,0)^{\mathtt{t}}$ at the top end of $\Omega$ and
no-slip condition everywhere on the boundary. In addition, the body force term
is $\textbf{f}=\textbf{0}$. In Fig. 3 we display the computed velocity
components and pressure on hexagon mesh with $h=3.03$e-2 and $k=0$, which
confirm the obtained results in [24].

---
Figure 3: Example 2. The numerical velocity and pressure for $\nu=1$ (first
row) and $\nu=0.01$ (second row) .

###  6.3 Example 4: Fluid flows with heterogeneous porous inclusions

Mathematical modeling and simulation of fluid flows in the presence of single
or multiple obstacles have been topics of interest for several decades due to
their wide applicability in various practical circumstances across
disciplines. Flow past solid bodies such as cylinders and airfoils has been
investigated broadly for a long time (see for instance, [15, 26, 30]) by using
Navier-Stokes equations. We study the time-dependent Navier-Stokes equation,
which is dependent on the porosity parameter, within the exterior flow domain
($\Omega_{f}$) characterized by large porosity and in the porous subdomains
($\Omega_{p}$) characterized by small porosity. A typical computational domain
in the present problem is illustrated in Fig. 4-(a). Here, denoting domain by
$\Omega=(-1,1)^{2}$ and the final time by $t_{F}$, we consider the time-
dependent Navier-Stokes equation with the porosity $\phi$ using the following
non-dimensionalization

Figure 4: Example 4. (a). Schematic of of the computational domain
$\Omega=\Omega_{p}\cup\Omega_{f}$ , where $\Omega_{f}$ is the large porosity
subdomain and $\Omega_{p}$ is subdomain with the small porosity. (b). The
polygonal mesh of domain.

$\begin{array}[]{rcll}\dfrac{1}{\phi}\left(\mathbf{u}_{t}+\mathbf{u}\cdot\boldsymbol{\nabla}\mathbf{u}-\nu\boldsymbol{\Delta}\mathbf{u}\right)+\nabla
p&=&\textbf{f}&\quad{\quad\hbox{in}\quad}\Omega\times(0,t_{F}]\,,\\\\[8.61108pt]
{\mathrm{div}}(\boldsymbol{u})&=&0&\quad{\quad\hbox{in}\quad}\Omega\times(0,t_{F}]\,,\\\\[8.61108pt]
{\mathbf{u}}(\cdot,0)&=&\mathbf{0}&\quad{\quad\hbox{in}\quad}\Omega\,.\end{array}$
(6.2)

On the left boundary, we prescribe a inflow velocity
$\mathbf{u}_{\mathtt{in}}=(0,1)^{\mathtt{t}}$. On the right boundary, we
impose a zero normal Cauchy stress, which means that we need to set

$\big{(}{\boldsymbol{\sigma}}+{\mathbf{u}}\otimes{\mathbf{u}}\big{)}{\mathbf{n}}\,=\,\mathbf{0}\qquad{\quad\hbox{on}\quad}\Gamma_{out}\,,$

and on the remainder of the boundary we set no-slip velocity
${\mathbf{u}}=\mathbf{0}$. In addition, we performed discretization of time by
Backward Euler method and taken 24 circular inclusions with different radii at
relatively random locations. The performance of the proposed method has been
tested with the following data:

$\begin{array}[]{llll}\phi\,=\,\left\\{\begin{array}[]{ll}0.2&{\quad\hbox{on}\quad}\Omega_{p}\,,\\\\[2.84526pt]
1&{\quad\hbox{on}\quad}\Omega_{f}\,,\end{array}\right.&\quad\nu\,=\,0.01\,,&t_{F}\,=\,0.1\,,&\Delta
t\,=\,2e-3\,.\end{array}$

Figure 5: Example 4. The numerical velocity and the first two components of
stress from left to right.

In Fig. 4-(b), we depicted the polygonal mesh of the computational domain. The
results of the desired case simulation are presented in Fig. 5. It can be seen
that the fluid is flowing faster in the area that has a large porosity; for
the area that has a small porosity, the fluid is flowing slowly. In the area
that has small porosity, we can see the gradation motion of the fluid clearly;
this fact emphasizes that the proposed scheme can deal with the irregular
pattern of porosity.

##  7 Declaration

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] L. Beirão da Veiga, C. Lovadina and G. Vacca, Virtual elements for the Navier–Stokes problem on polygonal meshes, SIAM J. Numer. Anal. 56 (2018) 1210–1242.
  * [2] L. Beirão da Veiga, D. Mora and G. Vacca, The Stokes complex for virtual elements with application to Navier–Stokes flows, J. Sci. Comput. 81 (2019) 990–1018.
  * [3] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, in: Springer Series in Computational Mathematics, vol. 15, Springer-Verlag, New York, 1991.
  * [4] D. Boffi, F. Brezzi and M. Fortin, Mixed Finite Element Methods and Applications, in: Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013.
  * [5] M. A. Behr, L.P. Franca and T.E. Tezduyar, Stabilized finite element methods for the velocity-pressure-stress formulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 104 (1) (1993) 31–48.
  * [6] J. Camaño, R. Oyarzúa, and G. Tierra,Analysis of an augmented mixed-FEM for the Navier-stokes problem, Math. Comput. 86(304) (2017), 589–615.
  * [7] J. Camaño, G. N. Gatica, R. Oyarzúa, and G. Tierra,An augmented mixed finite element method for the Navier-Stokes equations with variable viscosity, SIAM J. Numer. Anal. 54(2) (2016), 1069–1092.
  * [8] J. Camaño, C. García and R. Oyarzúa,Analysis of a momentum conservative mixed-FEM for the stationary Navier-Stokes problem, Numer. Methods Partial Differ. Equ. 37 (2021) 2895–2923.
  * [9] J. Camaño, G. N. Gatica, R. Oyarzúa and R. Ruiz-Baier, An augmented stress-based mixed finite element method for the steady state Navier–Stokes equations with nonlinear viscosity, Numer. Methods Partial Differ. Equ. 33 (2017) 1692–1725.
  * [10] Z. Cai and Y. Wang, Pseudostress-velocity formulation for incompressible Navier-stokes equations, Int. J. Numer. Methods Fluids 63(3) (2010), 341–356.
  * [11] G. Chen , M. Feng and X. Xie, Robust globally divergence-free weak Galerkin methods for Stokes equations, J. Comput. Math. 34 (5) (2016) 549–572.
  * [12] M. Dehghan and Z. Gharibi, Numerical analysis of fully discrete energy stable weak Galerkin finite element scheme for a coupled Cahn-Hilliard-Navier-Stokes phase-field model, Appl. Math. Comput. 410 (2021) 126487.
  * [13] M. Dehghan and Z. Gharibi, An analysis of weak Galerkin finite element method for a steady state Boussinesq problem, J. Comput. Appl. Math., 406 (2022), 114029.
  * [14] A. Ern and J.-L. Guermond: Theory and Practice of Finite Elements. Applied Mathematical Sciences, 159. Springer, New York (2004).
  * [15] B. Fornberg. A numerical study of steady viscous flow past a circular cylinder. Journal of Fluid Mechanics, 98 (4) (1980), 819–855.
  * [16] G. N. Gatica, A. Márquez and M.A. Sánchez, Analysis of a velocity-pressure-pseudostress formulation for the stationary Stokes equations, Comput. Methods Appl. Mech. Engrg. 199 (17-20) (2010) 1064–1079.
  * [17] G. N. Gatica, M. Munar and F. A. Sequeira, A mixed virtual element method for the Navier-Stokes equations, Math. Models Methods Appl. Sci. 28 (2018) 2719–2762.
  * [18] G. N. Gatica and F. A. Sequeira, An Lp spaces-based mixed virtual element method for the two-dimensional Navier–Stokes equations, Math. Models Methods Appl. Sci. 31 (14) (2021) 2937–2977.
  * [19] G. N. Gatica, A Simple Introduction to the Mixed Finite Element Method : Theory and Applications, Springer Briefs in Mathematics (Springer, 2014).
  * [20] V. Girault and P. Raviart, Finite Element Methods for Navier–Stokes Equations. Theory and Algorithms, Springer Series in Computational Mathematics, Vol. 5. Berlin: Springer, 1986.
  * [21] Z. Gharibi, M. Dehghan and M. Abbaszadeh, Numerical analysis of locally conservative weak Galerkin dual-mixed finite element method for the time-dependent Poisson-Nernst-Planck system, Comput. Math. Appl. 92 (2021) 88–108.
  * [22] Z. Gharibi, A weak Galerkin pseudostress-based mixed finite element method on polygonal meshes: application to the Brinkman problem appearing in porous media, Numer. Algorithms. (2024) 1–26.
  * [23] J. S. Howell and N. Walkington, Dual mixed finite element methods for the Navier–stokes equations, ESAIM: Math. Modell. Numer. Anal. 47 (2013), 789–805.
  * [24] X. Hu, L. Mu and X. Ye, A weak Galerkin finite element method for the Navier-Stokes equations, J. Comput. Appl. Math. 362 (2019) 614–625.
  * [25] X. Liu , J. Li and Z. Chen, A weak Galerkin finite element method for the Navier Stokes equations, J. Comput. Appl. Math. 333 (2019) 442–457.
  * [26] G. Li and J. AC Humphrey. Numerical modelling of confined flow past a cylinder of square cross-section at various orientations. International journal for numerical methods in fluids, 20 (11) (1995), 1215–1236.
  * [27] L. Mu , X. Wang and X. Ye, A modified weak Galerkin finite element method for the Stokes equations, J. Comput. Appl. Math. 275 (2015) 79–90.
  * [28] 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.
  * [29] L. Mu, A uniformly robust $H(\operatorname{div})$ weak Galerkin finite element method for Brinkman problems, SIAM J. Numer. Anal. 58 (3) (2020) 1422–1439.
  * [30] A. Sohankar, C. Norberg, and L. Davidson. Low-Reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet boundary condition. International journal for numerical methods in fluids, 26 (1) (1998), 39–56.
  * [31] J. Wang and X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115.
  * [32] J. Wang and X. Ye,A weak Galerkin mixed finite element method for second order elliptic problems, Math. Comput., 83 (2014), 2101–2126.
  * [33] T. Zhang and T. Lin, An analysis of a weak Galerkin finite element method for stationary Navier-Stokes problems, J. Comput. Appl. Math. 362 (2019) 484–497.
  * [34] B. Zhang, Y. Yang and M. Feng, A ${C^{0}}$-weak Galerkin finite element method for the two-dimensional Navier-Stokes equations in stream-function formulation, J. Comput. Math. 38 (2) (2020) 310–336.
  * [35] J. Zhang, K. Zhang, J. Li and X. Wang, A weak Galerkin finite element method for the Navier-Stokes equations, Commun. Comput. Phys. 23 (3) (2018) 706–746.
  * [36] Q. Zhai, R. Zhang and L. Mu,A new weak Galerkin finite element scheme for the Brinkman model, Commun. Comput. Phys. 19 (5) (2016) 1409–1434.