# Pressure-robust staggered DG methods for the Navier-Stokes equations on
general meshes††thanks: Submitted to the editors Month day, year.

Dohyun Kim School of Mathematics and Computing (Computational Science and
Engineering), Yonsei University, Seoul 03722, Korea (kim92n@yonsei.ac.kr).
Lina Zhao Department of Mathematics, City University of Hong Kong, Kowloon
Tong, Hong Kong SAR, China (linazha@cityu.edu.hk).  Eric Chung Department of
Mathematics, The Chinese University of Hong Kong, Hong Kong SAR, China
<EMAIL_ADDRESS> Eun-Jae Park School of Mathematics and Computing
(Computational Science and Engineering), Yonsei University, Seoul 03722, Korea
(ejpark@yonsei.ac.kr).

###### Abstract

In this paper, we design and analyze staggered discontinuous Galerkin methods
of arbitrary polynomial orders for the stationary Navier-Stokes equations on
polygonal meshes. The exact divergence-free condition for the velocity is
satisfied without any postprocessing. The resulting method is pressure-robust
so that the pressure approximation does not influence the velocity
approximation. A new nonlinear convective term that earning non-negativity is
proposed. The optimal convergence estimates for all the variables in $L^{2}$
norm are proved. Also, assuming that the rotational part of the forcing term
is small enough, we are able to prove that the velocity error is independent
of the Reynolds number and of the pressure. Furthermore, superconvergence can
be achieved for velocity under a suitable projection. Numerical experiments
are provided to validate the theoretical findings and demonstrate the
performances of the proposed method.

Keywords: Staggered grid, Discontinuous Galerkin method, Divergence free,
Pressure robustness, Navier-Stokes equations, General meshes, Superconvergence

##  1 Introduction

The Navier-Stokes equations play an important role in fluid dynamics. A number
of finite element methods were proposed to solve the incompressible Navier-
Stokes equations and accomplished many advances in the past half-century [31,
43, 26, 19, 21, 41, 9, 40]. Traditional inf-sup stable finite element methods
do not give robust velocity approximation when large irrotational force is
considered in general. Such methods yield velocity approximations which depend
on a pressure-dependent error contribution $\nu^{-1}\inf_{q\in
Q_{h}}\|p-q_{h}\|_{L^{2}(\Omega)}$ where $\nu$ is the viscosity and $Q_{h}$ is
the discrete pressure space. The influence of the pressure-dependent error
contribution is most pronouncing in the no-flow example which was first
considered in [25]. One of the main reason for the pressure-dependent error
contribution is that the velocity approximation does not satisfy the
incompressibility condition. To satisfy the incompressibility condition,
divergence-free finite element methods were proposed, see, for example [20,
45, 32, 33, 29, 44]. It is by no means trivial to construct finite element
spaces that satisfy inf-sup condition and at the same time yield divergence
free velocity. One of the approach exploited is to enrich the velocity space
without violating inf-sup condition. However, the construction of a suitable
space to enrich velocity space is tricky. As an alternative, one can use
divergence free velocity reconstruction operator to modify the right hand side
[4, 38, 27, 8]. When the nonlinear Navier-Stokes equations are considered, the
reconstruction operator should be applied to the velocity for each nonlinear
iteration and this leads to additional computational cost. For the case of
evolutionary incompressible Navier-Stokes equations, several pressure-robust
numerical approaches were considered such as $H(\mathrm{div}\,)$-conforming
discontinuous Galerkin (DG) [42, 34], $H^{1}$-conforming mixed finite element
method with grad-div stabilization [28, 1], and continuous interior penalty
methods [5]. While some of numerical simulations show optimal convergence
behavior, error estimates provided in the aforementioned works are suboptimal.

In recent years, a large effort has been devoted to the design and analysis of
discretization schemes that apply to general polygonal meshes. Among all the
methods, we mention several polygonal methods that have been designed for
arbitrary polygonal orders such as polygonal DG [6] methods, virtual element
methods (VEM) [2], hybrid high-order (HHO) [24] methods, and weak Galerkin
(WG) methods [46]. In the present work, we will devise a new approach within
the framework of staggered DG methods. Staggered DG methods as a new
generation of numerical schemes were firstly introduced by Chung and Engquist
for wave propagation on triangular meshes [14]. Since then it has been
successfully applied for various problems [15, 17, 16, 37, 12, 11, 13, 18].
Recently, Zhao and Park [49] extended this method to general polygonal meshes
for the Poisson equation. Then, a high order staggered DG method for general
second-order elliptic problems is developed in [51], and it is applied to
various physical problems arising from practical applications [55, 51, 36, 52,
47, 50, 54, 48]. For a pressure-robust method, Zhao et al. [53] proposed a
lowest-order pressure-robust staggered DG method for the Stokes equations.
They used a reconstruction operator based on an $H$(div) conforming function
space on polygonal meshes. However, its extension to high-order methods is not
trivial.

The objective of this paper is to design and analyze a pressure robust
staggered DG method of arbitrary polynomial orders on polygonal meshes for the
Navier-Stokes problem. The discrete formulation involves velocity gradient,
velocity and pressure, and the continuities for all the variables are
staggered on the inter-element boundaries in line with [47]. The staggered
continuity naturally gives an inter-element flux term which is free of
stabilization parameters and meanwhile ensures the inf-sup stability of the
resulting bilinear forms. The proposed method is pressure-robust without
enriching the velocity space nor introducing the divergence-preserving
reconstruction operator. To the best of the authors’ knowledge, this is the
first pressure-robust method on polygonal meshes without enriching the
discrete spaces and without establishing reconstruction operators. As such,
the construction of the method is relatively simple by using standard
polynomial spaces and the extension to arbitrary polynomial orders is
straightforward. Another novel contribution lies in the design of a new
nonlinear convective term that guarantees non-negativity. We prove that the
resulting method yields a divergence-free velocity approximation without any
postprocessing, which is particularly important for Navier-Stokes equations.
In addition, the unique solvability of the discrete solution is proved under a
smallness condition involving only the solenoidal part of the body force.
Importantly, the convergence estimate for velocity is proved to be independent
of the pressure and of the viscosity. The staggered DG method designed herein
offers some unique features which makes it advantageous over other polygonal
methods: It is locally conservative over each dual-element; superconvergence
can be obtained; it is stable without numerical flux nor penalty term; exact
divergence free condition is satisfied. It should be also emphasized that
developing a divergence free element on 3D is much more difficult and our
approach can be extended to 3D straightforwardly. The current paper focuses on
2D to present the core of the method and simplify the presentation.

The organization of the paper is as follows. In Section 2, we introduce the
pressure-robust staggered DG method on polygonal meshes and give some useful
lemmas including the divergence-free property of the velocity approximation.
In Section 3, we prove the existence and uniqueness of the nonlinear discrete
problem using fixed-point argument. Then a priori error estimates for all
variables are proved in Section 4. In Section 5, several numerical experiments
are conducted to verify the pressure robustness and the optimality of the
method. We end in Section 6 with some concluding remarks.

##  2 Staggered discontinuous Galerkin method

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded simply connected polygonal
domain with Lipschitz boundary $\partial\Omega$. For given data
$\bm{f}\in[L^{2}(\Omega)]^{2}$ and $\bm{g}\in[H^{1/2}(\partial\Omega)]^{2}$,
the incompressible Navier-Stokes problem in the conservative form seeks the
unknown velocity $\bm{u}$ and the pressure $p$ satisfying

$\displaystyle-\nu\Delta\bm{u}+\mathrm{div}\,(\bm{u}\otimes\bm{u})+\nabla p$
$\displaystyle=\bm{f}$ $\displaystyle\text{in }\Omega,$ (1)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0$ $\displaystyle\text{in
}\Omega,$ $\displaystyle\bm{u}$ $\displaystyle=\bm{g}$ $\displaystyle\text{on
}\partial\Omega,$ $\displaystyle\left\langle p\right\rangle_{\Omega}$
$\displaystyle=0.$

Here, $\left\langle
p\right\rangle_{\Omega}=\int_{\Omega}p\,\textup{d}\mathbf{x}$ and $\nu>0$ is a
real number representing the kinematic viscosity of the fluid. We introduce an
auxiliary variable $\bm{G}=\nu\nabla\bm{u}$, then (1) can be recast into the
following first order system

$\displaystyle\nu^{-1}\bm{G}-\nabla\bm{u}$ $\displaystyle=0,$
$\displaystyle\text{in }\Omega,$ (2)
$\displaystyle-\mathrm{div}\,\bm{G}+\mathrm{div}\,(\bm{u}\otimes\bm{u})+\nabla
p$ $\displaystyle=\bm{f},$ $\displaystyle\text{in }\Omega,$
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$ $\displaystyle\text{in
}\Omega,$ $\displaystyle\bm{u}$ $\displaystyle=\bm{g},$ $\displaystyle\text{on
}\partial\Omega,$ $\displaystyle\left\langle p\right\rangle_{\Omega}$
$\displaystyle=0.$

In the remainder of this paper, we assume $\bm{g}=\bm{0}$ for simplicity.

We denote the Sobolev space $W^{k,p}(\Omega)$ and we write $H^{k}(\Omega)$
when $p=2$. Those spaces are equipped with norm
$\left\|\cdot\right\|_{W^{k,p}(\Omega)}$ and
$\left\|\cdot\right\|_{H^{k}(\Omega)}$. $W^{k,p}_{0}(\Omega)$ and
$H^{k}_{0}(\Omega)$ denote the closure of $C^{\infty}_{c}(\Omega)$ with
respective norms. Here, $C^{\infty}_{c}(\Omega)$ is the infinitely
differentiable function spaces with compact support. We also define

$H(\mathrm{div}\,;\Omega)=\\{\bm{v}\in[L^{2}(\Omega)]^{2}:\nabla\cdot\bm{v}\in
L^{2}(\Omega)\\},\quad L^{2}_{0}(\Omega)=\\{q\in L^{2}(\Omega):\left\langle
q\right\rangle_{\Omega}=0\\}.$

In the sequel we use $C$ to denote a generic positive constant which may have
different values at different occurrences. The weak solution to (2) is defined
by

$\displaystyle(\nu^{-1}\bm{G}-\nabla\bm{u},\bm{H})$ $\displaystyle=0$
$\displaystyle\forall\bm{H}\in[L^{2}(\Omega)]^{2\times 2},$ (3)
$\displaystyle(\bm{G}-\bm{u}\otimes\bm{u}-p\bm{I},\nabla\bm{v})$
$\displaystyle=(\bm{f},\bm{v})$
$\displaystyle\forall\bm{v}\in[H^{1}_{0}(\Omega)]^{2},$
$\displaystyle(\nabla\cdot\bm{u},q)$ $\displaystyle=0$ $\displaystyle\forall
q\in L^{2}_{0}(\Omega).$

We can prove the following stability estimate proceeding similarly to [8] and
the proof is omitted for simplicity.

######  Lemma 2.1.

For $\bm{f}\in[L^{2}(\Omega)]^{2}$, let $\bm{f}=\bm{f}_{0}+\nabla\chi$ be the
Helmholtz-Hodge decomposition with $(\bm{f}_{0},\chi)\in
H(\mathrm{div}\,;\Omega)\times(H^{1}(\Omega)\cap L^{2}_{0}(\Omega))$. Further,
$\bm{f}_{0}$ satisfies $\nabla\cdot\bm{f}_{0}=0$. If $(\bm{G},\bm{u},p)$ is
the solution to (2), then, there holds

$\left\|\bm{u}\right\|_{H^{1}(\Omega)}\leq
C\nu^{-1}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}.$

In other words, $\bm{u}$ is independent of the irrotational force
$\nabla\chi$.

We now introduce some basic notations regarding staggered meshes that will be
exploited in the construction of the staggered DG method. Let $\mathcal{S}$ be
a star-shaped polygonal partition of $\Omega$, and the set of all the edges is
called primal-edges and is denoted by $\mathcal{F}_{pr}$. In addition, we use
$\mathcal{F}_{pr}^{i}$ and $\mathcal{F}_{pr}^{b}$ to stand for the interior
edges and boundary edges, respectively. For each primal-element
$T\in\mathcal{S}$, we select an interior point $\bm{x}$ and connect it to all
the vertices of $T$, thereby sub-triangular grids are generated during this
process. The resulting triangulation is denoted as $\mathcal{T}_{h}$ and the
edges generated during the subdivision process are called dual-edges and are
denoted by $\mathcal{F}_{dl}$. We rename the union of sub-triangles sharing
the common vertex $\bm{x}$ as $S(\bm{x})$. Here, we choose $\bm{x}$ to be an
interior point of the kernel of $S(\bm{x})\in\mathcal{S}$ and we use
$\mathcal{N}$ to represent the union of all interior points $\bm{x}$. For each
interior edge $e\in\mathcal{F}_{pr}^{i}$, we use $\mathcal{D}(e)$ to stand for
the dual-mesh, which is the union of two triangles in $\mathcal{T}_{h}$
sharing the common edge $e$. For each boundary edge
$e\in\mathcal{F}_{pr}^{b}$, we use $\mathcal{D}(e)$ to denote a triangle in
$\mathcal{T}_{h}$ having the edge $e$, see Figure 1 for an illustration. For
each edge $e\in\mathcal{F}_{h}=\mathcal{F}_{pr}\cup\mathcal{F}_{dl}$, we
define a unique normal vector $\bm{n}$ and tangential vector $\bm{t}$ by the
outward normal vector and counter-clockwise tangential vector of one of its
neighboring element, respectively.

Figure 1:  Schematic of primal- and dual-meshes. Solid lines are primal-edges
$\mathcal{F}_{pr}$ and dashed lines are dual-edges $\mathcal{F}_{dl}$. A
polygon surrounded by primal-edges is called primal-element $S(\bm{x})$ and
quadrilaterals surrounded by dual-edges are called dual-element
$\mathcal{D}(e)$.

For later analysis, we employ the general mesh regularity assumption (cf. [7,
49]): For every element $S(\bm{x})\in\mathcal{S}$ and every edge $e\in\partial
S(\bm{x})$, it satisfies $h_{e}\geq\rho_{E}h_{S(\bm{x})}$ for a positive
constant $\rho_{E}$, where $h_{e}$ is the length of edge $e$ and
$h_{S(\bm{x})}$ is the diameter of $S(\bm{x})$. Second, each element
$S(\bm{x})$ in $\mathcal{S}$ is star-shaped with respect to a ball of radius
$\geq\rho_{B}h_{S(\bm{x})}$, where $\rho_{B}$ is a positive constant. We
remark that the above assumptions ensure that the triangulation
$\mathcal{T}_{h}$ is shape regular. We employ the general mesh regularity
assumption just for the sake of simplicity. In fact, our numerical results
indicate that the proposed method allows elements with arbitrarily small edges
and a rigorous analysis for a mesh with small edges will be present in our
future work.

The jump $\left\llbracket\cdot\right\rrbracket$ and the average
$\left\\{\\!\\!\left\\{\cdot\right\\}\\!\\!\right\\}$ is defined by
$\left\llbracket v\right\rrbracket_{e}=v|_{\tau_{+}}-v|_{\tau_{-}}$ and
$\left\\{\\!\\!\left\\{v\right\\}\\!\\!\right\\}_{e}=(v|_{\tau_{+}}+v|_{\tau_{-}})/2$
where $\tau_{+}$ and $\tau_{-}$ are two elements in $\mathcal{T}_{h}$ sharing
the common edge $e\in\mathcal{F}_{pr}^{i}\cup\mathcal{F}_{dl}$. For
$e\in\mathcal{F}_{pr}^{b}$, we simply take $\left\llbracket
v\right\rrbracket_{e}=\left\\{\\!\\!\left\\{v\right\\}\\!\\!\right\\}_{e}=v|_{\tau_{+}}$.
The subscript $e$ will be omitted when there is no ambiguity. We denote the
$L^{2}$-inner product by $(f,g)_{\tau}=\int_{\tau}fg\,\textup{d}\mathbf{x}$
for 2D and $\left\langle f,g\right\rangle_{e}=\int_{e}fg\,\textup{d}s$ for 1D.
When $\bm{f}$ and $\bm{g}$ are vectors (or tensors), then
$(\cdot,\cdot)_{\tau}$ and $\left\langle\cdot,\cdot\right\rangle_{e}$ are
defined by their component-wise sum. That is, when
$\bm{f},\bm{g}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$,

$(\bm{f},\bm{g})_{\tau}=\int_{\tau}\bm{f}\cdot\bm{g}\,\textup{d}\mathbf{x},\quad\left\langle\bm{f},\bm{g}\right\rangle_{e}=\int_{e}\bm{f}\cdot\bm{g}\,\textup{d}s,$

and when $\bm{f},\bm{g}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2\times 2}$,

$(\bm{f},\bm{g})_{\tau}=\int_{\tau}\bm{f}:\bm{g}\,\textup{d}\mathbf{x},\quad\left\langle\bm{f},\bm{g}\right\rangle_{e}=\int_{e}\bm{f}:\bm{g}\,\textup{d}s.$

Here, $A:B=\sum_{ij}A_{ij}B_{ij}$ is the Frobenius inner product. The discrete
spaces are defined based on the definition described above. Let
$(H_{h},V_{h},Q_{h})$ be the discrete spaces defined by

$\displaystyle H_{h}$
$\displaystyle=\\{{\bm{H}}\in[\mathbb{P}_{k}(\mathcal{T}_{h})]^{2\times
2}:\left\llbracket\bm{H}\mathbf{n}\right\rrbracket_{e}=\bm{0}\;\forall
e\in\mathcal{F}_{pr}^{0},\;\left\llbracket\mathbf{t}\cdot\bm{H}\mathbf{n}\right\rrbracket_{e}=0\;\forall
e\in\mathcal{F}_{dl}\\},$ $\displaystyle V_{h}$
$\displaystyle=\\{\bm{v}\in[\mathbb{P}_{k}(\mathcal{T}_{h})]^{2}:\left\llbracket\bm{v}\cdot\mathbf{n}\right\rrbracket_{e}=0\;\forall
e\in\mathcal{F}_{dl}\\},$ $\displaystyle Q_{h}$
$\displaystyle=\\{q\in\mathbb{P}_{k}(\mathcal{T}_{h}):\left\llbracket
q\right\rrbracket_{e}=0\;\forall e\in\mathcal{F}_{pr}^{0}\\},$

where $\mathbb{P}_{k}(\mathcal{T}_{h})$ is a complete polynomial space of
degree less than or equal to $k$ on each triangle $\tau\in\mathcal{T}_{h}$.
The discrete space $V_{h}$ is equipped with the norm

$\left\|\bm{v}\right\|_{h}^{2}=\left\|\nabla\bm{v}\right\|_{L^{2}(\mathcal{T}_{h})}^{2}+\sum_{e\in\mathcal{F}_{pr}}h_{e}^{-1}\left\|\left\llbracket\bm{v}\right\rrbracket\right\|_{L^{2}(e)}^{2}+\sum_{e\in\mathcal{F}_{dl}}h_{e}^{-1}\left\|\left\llbracket\bm{v}\cdot\mathbf{t}\right\rrbracket\right\|_{L^{2}(e)}^{2}.$

Here, $\left\|\cdot\right\|_{L^{2}(\mathcal{T}_{h})}$ is the discrete
$L^{2}$-norm on the triangulation $\mathcal{T}_{h}$. We also define the
discrete $L^{4}$-norm for $V_{h}$ by

$\left\|\bm{v}\right\|_{0,4,h}^{4}=\left\|\bm{v}\right\|_{L^{4}(\mathcal{T}_{h})}^{4}+\sum_{e\in\mathcal{F}_{h}}h_{e}^{-1}\left\|\left\\{\\!\\!\left\\{\bm{v}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}^{4}$
(4)

and discrete $H^{1}$-seminorm for $Q_{h}$ by

$\left\|q\right\|_{1,h}^{2}=\left\|\nabla
q\right\|_{L^{2}(\mathcal{T}_{h})}^{2}+\sum_{e\in\mathcal{F}_{dl}}h_{e}^{-1}\left\|\left\llbracket
q\right\rrbracket\right\|_{L^{2}(e)}^{2}.$

To impose the mean zero condition for the pressure variable, we introduce

$Q_{h}^{0}=\\{q\in Q_{h}:\left\langle q\right\rangle_{\Omega}=0\\}.$

Note that $\left\|\cdot\right\|_{1,h}$ is a norm on $Q_{h}^{0}$.

In the following, we introduce the (discrete) trace inequality and the
discrete Sobolev embedding theorem.

######  Lemma 2.2.

There exists $C$ independent of $h$ such that for all $e\in\mathcal{F}_{h}$
and $e\subset\partial\tau\in\mathcal{T}_{h}$ (cf. [35])

$\displaystyle\left\|\bm{v}\right\|_{L^{4}(e)}\leq
C(h_{\tau}^{-1/4}\left\|\bm{v}\right\|_{L^{4}(\tau)}+h_{\tau}^{1/4}\left\|\nabla\bm{v}\right\|_{L^{2}(\tau)})\quad\forall\bm{v}\in
H^{1}(\tau).$ (5)

Also, we have (cf. [23])

$\left\|\bm{v}\right\|_{L^{2}(e)}\leq
Ch_{\tau}^{-1/2}\left\|\bm{v}\right\|_{L^{2}(\tau)}\quad\forall\bm{v}\in
V_{h}.$

######  Lemma 2.3.

[22] There exists $C$ independent of $h$ such that

$\left\|\bm{v}\right\|_{0,4,h}\leq
C\left\|\bm{v}\right\|_{h}\quad\forall\bm{v}\in V_{h}$

and

$\left\|\bm{v}\right\|_{L^{q}(\Omega)}\leq C\left\|\bm{v}\right\|_{h}\quad
1\leq q\leq 6,\;\forall\bm{v}\in V_{h}.$ (6)

We can define the following degrees of freedom in the spirit of [15, 47].

######  Lemma 2.4 (Degrees of freedom).

Any function $\bm{H}\in H_{h}$ is uniquely determined by the following degrees
of freedom:

  * (XD1)

For $e\in\mathcal{F}_{pr}$, we have

$\Phi^{pr}_{e}(\bm{H}):=\left\langle\bm{H}\mathbf{n},\bm{p}_{k}\right\rangle_{e}\quad\forall\bm{p}_{k}\in[\mathbb{P}_{k}(e)]^{2}.$

  * (XD2)

For $e\in\mathcal{F}_{dl}$, we have

$\Phi^{dl}_{e}(\bm{H}):=\left\langle\mathbf{t}\cdot\bm{H}\mathbf{n},p_{k}\right\rangle_{e}\quad\forall
p_{k}\in\mathbb{P}_{k}(e).$

  * (XD3)

For each $\tau\in\mathcal{T}_{h}$, we can obtain

$\Phi_{\tau}(\bm{H}):=(\bm{H},\bm{p}_{k-1})_{\tau}\quad\forall\bm{p}_{k-1}\in[\mathbb{P}_{k-1}(\tau)]^{2\times
2}.$

Similarly, any function $\bm{v}\in V_{h}$ is uniquely determined by the
following degrees of freedom:

  * (VD1)

For $e\in\mathcal{F}_{dl}$, we have

$\bm{\phi}_{e}(\bm{v}):=\left\langle\bm{v}\cdot\mathbf{n},p_{k}\right\rangle_{e}\quad\forall
p_{k}\in\mathbb{P}_{k}(e).$

  * (VD2)

For each $\tau\in\mathcal{T}_{h}$, we can obtain

$\bm{\phi}_{\tau}(\bm{v}):=(\bm{v},\bm{p}_{k-1})_{\tau}\quad\forall\bm{p}_{k-1}\in[\mathbb{P}_{k-1}(\tau)]^{2}.$

Finally, any function $q\in Q_{h}$ is uniquely determined by

  * (SD1)

For $e\in\mathcal{F}_{pr}$, we have

$\phi_{e}(q):=\left\langle q,p_{k}\right\rangle_{e}\quad\forall
p_{k}\in\mathbb{P}_{k}(e).$

  * (SD2)

For each $\tau\in\mathcal{T}_{h}$, we can obtain

$\phi_{\tau}(q):=(q,p_{k-1})_{\tau}\quad\forall
p_{k-1}\in\mathbb{P}_{k-1}(\tau).$

Using the degrees of freedom defined in Lemma 2.4, we define the interpolation
operators so that

$\bm{\Phi}(\bm{H}-\Pi_{h}\bm{H})=0,\quad\bm{\phi}(\bm{v}-J_{h}\bm{v})=0,\quad\phi(q-I_{h}q)=0.$

By the polynomial preserving property of the interpolation operators, we
obtain the following lemma with minor modification of [15].

######  Lemma 2.5.

There exists $C$ independent of $h$ such that

$\displaystyle\left\|\bm{H}-\Pi_{h}\bm{H}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch^{k+1}\left\|\bm{H}\right\|_{H^{k+1}(\Omega)}$
$\displaystyle\forall\bm{H}\in[H^{k+1}(\Omega)]^{2\times 2},$
$\displaystyle\left\|\bm{v}-J_{h}\bm{v}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch^{k+1}\left\|\bm{v}\right\|_{H^{k+1}(\Omega)}$
$\displaystyle\forall\bm{v}\in[H^{k+1}(\Omega)]^{2},$
$\displaystyle\left\|\bm{v}-J_{h}\bm{v}\right\|_{0,4,h}$ $\displaystyle\leq
Ch^{k+1}\left\|\bm{v}\right\|_{W^{k+1,4}(\Omega)}$
$\displaystyle\forall\bm{v}\in[W^{k+1,4}(\Omega)]^{2},$
$\displaystyle\left\|q-I_{h}q\right\|_{L^{2}(\Omega)}$ $\displaystyle\leq
Ch^{k+1}\left\|q\right\|_{H^{k+1}(\Omega)}$ $\displaystyle\forall q\in
H^{k+1}(\Omega).$

Then the staggered discontinuous Galerkin method for (1) is defined as
follows: Find $(\bm{G}_{h},\bm{u}_{h},p_{h})\in H_{h}\times V_{h}\times
Q_{h}^{0}$ such that

$\displaystyle\nu^{-1}(\bm{G}_{h},\bm{H})+B_{h}^{*}(\bm{u}_{h},\bm{H})$
$\displaystyle=\sum_{e\in\mathcal{F}_{pr}^{b}}\left\langle\bm{g},\bm{H}\mathbf{n}\right\rangle_{e}$
$\displaystyle\forall\bm{H}\in H_{h},$ (7a) $\displaystyle
B_{h}(\bm{G}_{h},\bm{v})+N_{h}(\bm{u}_{h};\bm{u}_{h},\bm{v})+b_{h}^{*}(p_{h},\bm{v})$
$\displaystyle=(\bm{f},\bm{v})+\sum_{e\in\mathcal{F}_{h}^{b}}\left\langle|\bm{g}\cdot\mathbf{n}|-\bm{g}\cdot\mathbf{n},\bm{g}\cdot\bm{v}\right\rangle_{e}$
$\displaystyle\forall\bm{v}\in V_{h},$ (7b) $\displaystyle
b_{h}(\bm{u}_{h},q)$
$\displaystyle=\sum_{e\in\mathcal{F}_{pr}^{b}}\left\langle\bm{g}\cdot\mathbf{n},q\right\rangle_{e}$
$\displaystyle\forall q\in Q_{h}^{0}.$ (7c)

Here,

$\displaystyle B_{h}(\bm{H},\bm{v})$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}(\bm{H},\nabla\bm{v})_{\tau}-\sum_{e\in\mathcal{F}_{pr}}\left\langle\bm{H}\mathbf{n},\left\llbracket\bm{v}\right\rrbracket\right\rangle-\sum_{e\in\mathcal{F}_{dl}}\left\langle\mathbf{t}\cdot\bm{H}\mathbf{n},\left\llbracket\bm{v}\cdot\mathbf{t}\right\rrbracket\right\rangle_{e},$
(8) $\displaystyle B_{h}^{*}(\bm{v},\bm{H})$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}(\bm{v},\mathrm{div}\,\bm{H})_{\tau}-\sum_{e\in\mathcal{F}_{dl}}\left\langle\bm{v}\cdot\mathbf{n},\left\llbracket\mathbf{n}\cdot\bm{H}\mathbf{n}\right\rrbracket\right\rangle_{e},$
(9) $\displaystyle b_{h}(\bm{v},q)$
$\displaystyle=-\sum_{\tau\in\mathcal{T}_{h}}(\bm{v},\nabla
q)_{\tau}+\sum_{e\in\mathcal{F}_{dl}}\left\langle\bm{v}\cdot\mathbf{n},\left\llbracket
q\right\rrbracket\right\rangle_{e},$ (10) $\displaystyle b_{h}^{*}(q,\bm{v})$
$\displaystyle=-\sum_{\tau\in\mathcal{T}_{h}}(q,\nabla\cdot\bm{v})_{\tau}+\sum_{e\in\mathcal{F}_{pr}}\left\langle
q,\left\llbracket\bm{v}\cdot\mathbf{n}\right\rrbracket\right\rangle_{e},$ (11)
$\displaystyle N_{h}(\bm{w};\bm{\psi},\bm{v})$
$\displaystyle=-\sum_{\tau\in\mathcal{T}_{h}}(\bm{\psi}\otimes\bm{w},\nabla\bm{v})_{\tau}+\sum_{e\in\mathcal{F}_{h}^{0}}\left\langle\left\\{\\!\\!\left\\{\bm{w}\cdot\mathbf{n}\right\\}\\!\\!\right\\},\left\\{\\!\\!\left\\{\bm{\psi}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}\right\rrbracket\right\rangle_{e}$
$\displaystyle\quad+\sum_{e\in\mathcal{F}_{h}}\left\langle|\left\\{\\!\\!\left\\{\bm{w}\cdot\mathbf{n}\right\\}\\!\\!\right\\}|,\left\llbracket\bm{\psi}\right\rrbracket\cdot\left\llbracket\bm{v}\right\rrbracket\right\rangle_{e},$
(12)

for any $(H,\bm{v},q)\in H_{h}\times V_{h}\times Q_{h}^{0}$ and
$(\bm{\psi},\bm{w})\in V_{h}\times V_{h}$.

In the remainder of this paper, we take $\bm{g}=\bm{0}$ for simplicity. Note
that the formulation is based on the conservative form of the Navier-Stokes
equation and the upwind term is added to ensure the non-negativity, see Lemma
2.9. Also, integration by parts and the definitions of the discrete spaces
lead to the following discrete adjoint properties

$\displaystyle B_{h}(\bm{H},\bm{v})$ $\displaystyle=-B_{h}^{*}(\bm{v},\bm{H})$
$\displaystyle\forall(\bm{H},\bm{v})\in H_{h}\times V_{h},$ (13)
$\displaystyle b_{h}(\bm{v},q)$ $\displaystyle=-b_{h}^{*}(q,\bm{v})$
$\displaystyle\forall(\bm{v},q)\in V_{h}\times Q_{h}.$

By the definition of the discrete bilinear forms and the interpolation
operators, we have the following lemma.

######  Lemma 2.6.

Assume that
$(\bm{G},\bm{u},p)\in([H(\mathrm{div}\,;\Omega)]^{2}\cap[H^{1/2+\epsilon}(\Omega)]^{2\times
2})\times[H^{1}(\Omega)]^{2}\times(L^{2}(\Omega)\cap
H^{1/2+\epsilon}(\Omega))$ with $\epsilon>0$. Then, the following equalities
hold:

$\displaystyle B_{h}(\Pi_{h}\bm{G}-\bm{G},\bm{v})$ $\displaystyle=0$
$\displaystyle\forall\bm{v}\in V_{h},$ $\displaystyle
b_{h}(J_{h}\bm{u}-\bm{u},q)$ $\displaystyle=0$ $\displaystyle\forall q\in
Q_{h},$ $\displaystyle B_{h}^{*}(J_{h}\bm{u}-\bm{u},\bm{H})$ $\displaystyle=0$
$\displaystyle\forall\bm{H}\in H_{h},$ $\displaystyle
b_{h}^{*}(I_{h}p-p,\bm{v})$ $\displaystyle=0$ $\displaystyle\forall\bm{v}\in
V_{h}.$

We have the discrete inf-sup condition for $B_{h}(\cdot,\cdot)$ and
$b_{h}(\cdot,\cdot)$ (cf. [15, 47]).

######  Lemma 2.7.

There exists $C_{B}$ and $C_{b}$ independent of $h$ such that

$0<C_{B}\leq\adjustlimits{\inf}_{\bm{v}\in V_{h}}{\sup}_{\bm{H}\in
H_{h}}\frac{B_{h}(\bm{H},\bm{v})}{\left\|\bm{H}\right\|_{L^{2}(\Omega)}\left\|\bm{v}\right\|_{h}},\quad
0<C_{b}\leq\adjustlimits{\inf}_{q\in Q_{h}}{\sup}_{\bm{v}\in
V_{h}}\frac{b_{h}^{*}(q,\bm{v})}{\left\|\bm{v}\right\|_{L^{2}(\Omega)}\left\|q\right\|_{1,h}}.$

Integration by parts reveals the consistency of the nonlinear trilinear form
$N_{h}(\cdot;\cdot,\cdot)$.

######  Lemma 2.8 (Consistency).

Let $\bm{u}\in[H^{1}_{0}(\Omega)]^{2}$ be the solution to (1). Then the
following holds:

$N_{h}(\bm{u};\bm{u},\bm{v}_{h})=(\mathrm{div}\,(\bm{u}\otimes\bm{u}),\bm{v}_{h})\quad\forall\bm{v}_{h}\in
V_{h}.$

In the next lemma, we state the non-negative property for
$N_{h}(\cdot;\cdot,\cdot)$, which is crucial for the subsequent analysis.

######  Lemma 2.9 (Non-negativity).

Let $\bm{w}\in H(\mathrm{div}\,;\Omega)$ with $\nabla\cdot\bm{w}=0$. Then we
have

$N_{h}(\bm{w};\bm{v},\bm{v})\geq 0\quad\forall\bm{v}\in
V_{h}+[H^{1}_{0}(\Omega)]^{2}.$

###### Proof.

Let $\tau\in\mathcal{T}_{h}$ be given. From integration by parts, we have

$(\bm{v}\otimes\bm{w},\nabla\bm{v})_{\tau}=-(\mathrm{div}\,(\bm{v}\otimes\bm{w}),\bm{v})_{\tau}+\left\langle(\bm{v}\otimes\bm{w})\mathbf{n},\bm{v}\right\rangle_{\partial\tau}.$
(14)

Since $\nabla\cdot\bm{w}=0$, integration by parts and straightforward
computation yield

$(\mathrm{div}\,(\bm{v}\otimes\bm{w}),\bm{v})_{\tau}=\frac{1}{2}\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{\partial\tau}$

and

$\left\langle(\bm{v}\otimes\bm{w})\mathbf{n},\bm{v}\right\rangle_{\partial\tau}=\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{\partial\tau}.$

Then (14) can be rewritten as

$(\bm{v}\otimes\bm{w},\nabla\bm{v})_{\tau}=\frac{1}{2}\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{\partial\tau}.$

Summing over all the elements $\tau\in\mathcal{T}_{h}$, we have

$\displaystyle\sum_{\tau\in\mathcal{T}_{h}}(\bm{v}\otimes\bm{w},\nabla\bm{v})_{\tau}$
$\displaystyle=\frac{1}{2}\sum_{\tau\in\mathcal{T}_{h}}\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{\partial\tau}$
$\displaystyle=\sum_{e\in\mathcal{F}_{h}^{0}}\left\langle\bm{w}\cdot\mathbf{n},\left\\{\\!\\!\left\\{\bm{v}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}\right\rrbracket\right\rangle_{e}+\frac{1}{2}\sum_{e\in\mathcal{F}_{h}^{b}}\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{e}.$

This and the definition of $N_{h}(\cdot;\cdot,\cdot)$ yield

$\displaystyle N_{h}(\bm{w};\bm{v},\bm{v})$
$\displaystyle=-\frac{1}{2}\sum_{e\in\mathcal{F}_{h}^{b}}\left\langle\bm{w}\cdot\mathbf{n},\bm{v}\cdot\bm{v}\right\rangle_{e}+\sum_{e\in\mathcal{F}_{h}}\left\langle|\bm{w}\cdot\mathbf{n}|,\left\llbracket\bm{v}\right\rrbracket\cdot\left\llbracket\bm{v}\right\rrbracket\right\rangle_{e}$
$\displaystyle\geq\frac{1}{2}\sum_{e\in\mathcal{F}_{h}}\left\langle|\bm{w}\cdot\mathbf{n}|,\left\llbracket\bm{v}\right\rrbracket\cdot\left\llbracket\bm{v}\right\rrbracket\right\rangle_{e}$
$\displaystyle\geq 0.$

∎

By using Hölder’s inequality and the discrete Sobolev embedding theorem, we
obtain the following lemma.

######  Lemma 2.10 (Boundedness).

For any $\bm{z}_{h},\bm{v}_{h},\bm{w}_{h}\in V_{h}$, it holds

$N_{h}(\bm{z}_{h};\bm{v}_{h},\bm{w}_{h})\leq
C_{N}\left\|\bm{z}_{h}\right\|_{h}\left\|\bm{v}_{h}\right\|_{h}\left\|\bm{w}_{h}\right\|_{h}.$

To ease later analysis, we define the divergence-free subspace of $V_{h}$ by

$V_{h}^{\mathrm{div}\,}=\\{\bm{v}\in V_{h}:\bm{v}\in
H(\mathrm{div}\,;\Omega),\;\nabla\cdot\bm{v}=0\\}.$

This subspace plays an important role in subsequent sections. We close this
section by observing the connection between $b_{h}(\cdot,\cdot)$ and
$V_{h}^{\mathrm{div}}\,$.

######  Lemma 2.11.

If $\bm{v}_{h}\in V_{h}$ satisfies

$b_{h}(\bm{v}_{h},q)=0\quad\forall q\in Q_{h}^{0},$ (15)

then $\bm{v}_{h}\in H(\mathrm{div}\,;\Omega)$ and it is divergence-free.

###### Proof.

Note that

$b_{h}(\bm{v}_{h},c)=0$

for any constant $c$ by the definition. Therefore, the condition (15) is
equivalent to

$b_{h}(\bm{v}_{h},q)=0\quad\forall q\in Q_{h}.$

By (13), we obtain

$0=b_{h}(\bm{v}_{h},q)=b_{h}^{*}(q,\bm{v}_{h})\quad\forall q\in Q_{h}.$

Take $q$ so that (cf. (SD1)-(SD2))

$\displaystyle(q+\nabla\cdot\bm{v}_{h},p_{k-1})_{\tau}$ $\displaystyle=0$
$\displaystyle\forall
p_{k-1}\in\mathbb{P}_{k-1}(\tau),\;\tau\in\mathcal{T}_{h},$
$\displaystyle\left\langle
q-\left\llbracket\bm{v}_{h}\cdot\mathbf{n}\right\rrbracket,p_{k}\right\rangle_{e}$
$\displaystyle=0$ $\displaystyle\forall
p_{k}\in\mathbb{P}_{k}(e),\;e\in\mathcal{F}_{pr}.$

Then we obtain

$0=b_{h}^{*}(q,\bm{v}_{h})=\sum_{\tau\in\mathcal{T}_{h}}\left\|\nabla\cdot\bm{v}_{h}\right\|_{L^{2}(\tau)}^{2}+\sum_{e\in\mathcal{F}_{pr}}\left\|\left\llbracket\bm{v}_{h}\cdot\mathbf{n}\right\rrbracket\right\|_{L^{2}(e)}^{2}.$

Therefore, $\left\llbracket\bm{v}_{h}\cdot\mathbf{n}\right\rrbracket=0$ for
all $e\in\mathcal{F}_{pr}$. Since
$\left\llbracket\bm{v}_{h}\cdot\mathbf{n}\right\rrbracket_{e}=0$ for all
$e\in\mathcal{F}_{dl}$ by definition of $V_{h}$, we have $\bm{v}_{h}\in
H(\mathrm{div}\,;\Omega)$. Furthermore, since
$\nabla\cdot(\bm{v}_{h}|_{\tau})=0$ for all $\tau\in\mathcal{T}_{h}$, we have
$\nabla\cdot\bm{v}_{h}=0$. ∎

######  Remark 2.12.

(divergence free velocity). Our proposed scheme can yield a divergence-free
velocity by following Lemma 2.11, which is a desirable feature. Thanks to the
divergence-free property and the specially designed term for the nonlinear
convective term, we are able to prove that the convergence estimates are
independent of the pressure variable and the coefficient $\nu$ under a
suitable assumption on the source term $\bm{f}$. Unlike the existing works on
polygonal meshes [10, 27, 39, 53], we do not require velocity reconstruction,
which can greatly reduce the computational complexity and ease the
construction of the method.

##  3 Existence and uniqueness

In this section, we discuss the existence and uniqueness of the solution to
(7). A solution operator $T_{h}:V_{h}^{\mathrm{div}}\,\rightarrow V_{h}$ is
defined as follows: For given $\bm{z}_{h}\in V_{h}$, find
$\bm{w}_{h}=T_{h}(\bm{v}_{h})\in V_{h}$ so that

$A_{h}((\bm{S}_{h},\bm{w}_{h},r_{h}),(\bm{H},\bm{v},q))+N_{h}(\bm{z}_{h};\bm{w}_{h},\bm{v})=(\bm{f},\bm{v})\quad\forall(\bm{H},\bm{v},q)\in
H_{h}\times V_{h}\times Q_{h}$ (16)

for some $\bm{G}_{h}\in H_{h}$ and $p_{h}\in Q_{h}^{0}$. Here,

$A_{h}((\bm{S},\bm{w},r),(\bm{H},\bm{v},q))=\nu^{-1}(\bm{S},\bm{H})+B_{h}^{*}(\bm{w},\bm{H})+B_{h}(\bm{S},\bm{v})+b_{h}^{*}(r_{h},\bm{v})+b_{h}(\bm{w},q).$

Observe that finding the solution to (7) is equivalent to finding a fixed-
point $\bm{u}_{h}$ of $T_{h}$ so that

$T_{h}(\bm{u}_{h})=\bm{u}_{h}$

with its corresponding $\bm{G}_{h}$ and $p_{h}$.

######  Lemma 3.1.

$T_{h}$ is well-defined on $V_{h}^{\mathrm{div}}\,$ and $T_{h}(\bm{z}_{h})\in
V_{h}^{\mathrm{div}}\,$ for all $\bm{z}_{h}\in V_{h}^{\mathrm{div}}\,$.

###### Proof.

Let $\bm{z}_{h}\in V_{h}^{\mathrm{div}}\,$ be given. From Lemma 2.11, it is
clear that $T_{h}(\bm{z}_{h})\in V_{h}^{\mathrm{div}}\,$. By Lemma 2.9 and the
definition of $A_{h}$, we obtain

$A_{h}(\varphi,\varphi)+N_{h}(\bm{z}_{h};\bm{v},\bm{v})\geq\nu^{-1}\left\|\bm{H}\right\|_{L^{2}(\Omega)}^{2}\quad\forall\varphi=(\bm{H},\bm{v},q)\in
H_{h}\times V_{h}\times Q_{h}^{0}.$

Let $\varsigma_{i}=(\bm{G}_{i},\bm{w}_{i},p_{i})$ be the solution to (16) with
$\bm{z}_{h}$ for $i=1,2$. Then

$\nu^{-1}\left\|\bm{G}_{1}-\bm{G}_{2}\right\|_{L^{2}(\Omega)}^{2}\leq
A_{h}(\varsigma_{1}-\varsigma_{2},\varsigma_{1}-\varsigma_{2})+N_{h}(\bm{z}_{h};\bm{w}_{1}-\bm{w}_{2},\bm{w}_{1}-\bm{w}_{2})=0.$

Therefore, we have $\bm{G}_{1}=\bm{G}_{2}$. Furthermore, Lemma 2.7 and taking
$(\bm{v},q)=(\bm{0},0)$ in (16) lead to

$\left\|\bm{w}_{1}-\bm{w}_{2}\right\|_{h}\leq C\sup_{\bm{H}\in
H_{h}}\frac{B_{h}(\bm{H},\bm{w}_{1}-\bm{w}_{2})}{\left\|\bm{H}\right\|_{L^{2}(\Omega)}}=C\sup_{\bm{H}\in
H_{h}}\frac{B_{h}^{*}(\bm{w}_{1}-\bm{w}_{2},\bm{H})}{\left\|\bm{H}\right\|_{L^{2}(\Omega)}}=0.$

Finally, the discrete adjoint property (13) and (16) with
$(\bm{H},q)=(\bm{0},0)$ read

$0=b_{h}^{*}(p_{1}-p_{2},\bm{v})=b_{h}(\bm{v},p_{1}-p_{2})\quad\forall\bm{v}\in
V_{h}.$

By the inf-sup condition (2.7), we can obtain

$\displaystyle\left\|p_{1}-p_{2}\right\|_{1,h}\leq C\sup_{\bm{v}\in
V_{h}}\frac{b_{h}(\bm{v},p_{1}-p_{2})}{\left\|\bm{v}\right\|_{L^{2}(\Omega)}}=0.$

Therefore, we have $\varsigma_{1}=\varsigma_{2}$. Since $V_{h}$ is finite
dimensional space and $T_{h}$ is linear, $T_{h}$ is well-defined.

∎

Then, we have the following stability estimate.

######  Lemma 3.2.

For any $\bm{z}_{h}\in V_{h}^{\mathrm{div}}\,$, we have

$\left\|T_{h}(\bm{z}_{h})\right\|_{h}=\left\|\bm{w}_{h}\right\|_{h}\leq
C_{B}^{-2}\nu^{-1}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)},$

where $\bm{f}_{0}$ is from the Helmholtz-Hodge decomposition
$\bm{f}=\bm{f}_{0}+\nabla\chi$.

###### Proof.

Let $\bm{z}_{h}\in V_{h}^{\mathrm{div}}\,$ and $\bm{w}_{h}=T_{h}(\bm{z}_{h})$.
Then there exists $\varsigma=(\bm{G}_{h},\bm{w}_{h},p_{h})$ such that

$A_{h}(\varsigma,\varphi)+N_{h}(\bm{z}_{h};\bm{w}_{h},\bm{v})=(\bm{f},\bm{v})\quad\forall\varphi=(\bm{H},\bm{v},q)\in
H_{h}\times V_{h}\times Q_{h}^{0}.$

Taking $\varphi=\varsigma$, we obtain

$\nu^{-1}\left\|\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}\leq
A_{h}(\varsigma,\varsigma)+N_{h}(\bm{z}_{h};\bm{w}_{h},\bm{w}_{h})=(\bm{f},\bm{w}_{h}).$

From Lemma 2.7 and (16), we have

$C_{B}^{2}\left\|\bm{w}_{h}\right\|_{h}^{2}\leq\nu^{-2}\left\|\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}\leq\nu^{-1}(\bm{f},\bm{w}_{h}).$

Finally, recall that $\bm{w}_{h}=T_{h}(\bm{z}_{h})\in V_{h}^{\mathrm{div}}\,$
by Lemma 3.1. Then the Helmholtz-Hodge decomposition
$\bm{f}=\bm{f}_{0}+\nabla\chi$ yields

$C_{B}^{2}\left\|\bm{w}_{h}\right\|_{h}\leq\nu^{-1}\frac{(\bm{f}_{0}+\nabla\chi,\bm{w}_{h})}{\left\|\bm{w}_{h}\right\|_{h}}=\nu^{-1}\frac{(\bm{f}_{0},\bm{w}_{h})}{\left\|\bm{w}_{h}\right\|_{h}}\leq\nu^{-1}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}.$

Therefore, the proof is completed. ∎

######  Remark 3.3.

Lemma 3.2 implies that $T_{h}(\bm{z}_{h})$ is independent of the irrotational
part of $\bm{f}$.

By Lemma 3.2 and the Brouwer fixed point theorem, the existence of the fixed-
point $\bm{u}_{h}=T_{h}(\bm{u}_{h})$ is guaranteed. To show that the fixed-
point is unique, it suffices to show that $T_{h}$ is a contraction mapping.

######  Theorem 3.4.

Assume that $\bm{f}=\bm{f}_{0}+\nabla\chi$ satisfies

$\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}<C_{B}^{4}C_{N}^{-1}\nu^{2}.$

Then $T_{h}$ has a unique fixed-point in

$W_{h}^{\rho}=\\{\bm{v}\in
V_{h}^{\mathrm{div}}\,:\left\|\bm{v}\right\|_{h}\leq\rho<C_{B}^{2}C_{N}^{-1}\nu\\}.$

###### Proof.

Let $\bm{\varsigma}_{i}=(\bm{G}_{i},\bm{w}_{i},p_{i})$ be the solution to (16)
with $\bm{z}_{i}\in V_{h}^{\mathrm{div}}\,$ for $i=1,2$. Then we have

$A_{h}(\bm{\varsigma}_{i},\bm{\varphi})+N_{h}(\bm{z}_{i};\bm{w}_{i},\bm{v})=(\bm{f},\bm{v})\quad\forall\varphi=(\bm{H},\bm{v},q)\in
H_{h}\times V_{h}\times Q_{h}^{0}.$

Taking $\bm{\varphi}=\bm{\varsigma}_{1}-\bm{\varsigma}_{2}$, we obtain

$\displaystyle\nu^{-1}\left\|\bm{G}_{1}-\bm{G}_{2}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq
A_{h}(\bm{\varsigma}_{1}-\bm{\varsigma}_{2},\bm{\varsigma}_{1}-\bm{\varsigma}_{2})+N_{h}(\bm{z}_{1};\bm{w}_{1}-\bm{w}_{2},\bm{w}_{1}-\bm{w}_{2})$
$\displaystyle=N_{h}(\bm{z}_{2}-\bm{z}_{1};\bm{w}_{2},\bm{w}_{1}-\bm{w}_{2})$
$\displaystyle\leq
C_{N}\left\|\bm{z}_{1}-\bm{z}_{2}\right\|_{h}\left\|\bm{w}_{2}\right\|_{h}\left\|\bm{w}_{1}-\bm{w}_{2}\right\|_{h}.$

From Lemma 2.7 and Lemma 3.2, we obtain

$\displaystyle\left\|\bm{w}_{1}-\bm{w}_{2}\right\|_{h}^{2}$ $\displaystyle\leq
C_{B}^{-2}\nu^{-2}\left\|\bm{G}_{1}-\bm{G}_{2}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq
C_{B}^{-2}C_{N}\nu^{-1}\left\|\bm{z}_{2}-\bm{z}_{1}\right\|_{h}\left\|\bm{w}_{2}\right\|\left\|\bm{w}_{1}-\bm{w}_{2}\right\|_{h}$
$\displaystyle\leq
C_{B}^{-4}C_{N}\nu^{-2}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|\bm{z}_{2}-\bm{z}_{1}\right\|_{h}\left\|\bm{w}_{1}-\bm{w}_{2}\right\|_{h}.$

Recalling that $\bm{w}_{i}=T_{h}(\bm{z}_{i})$, there holds

$\displaystyle\left\|T_{h}(\bm{z}_{1})-T_{h}(\bm{z}_{2})\right\|_{h}\leq
C_{B}^{-4}\nu^{-2}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|\bm{z}_{2}-\bm{z}_{1}\right\|_{h}.$

By assuming
$\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}<C_{B}^{4}C_{N}^{-1}\nu^{2}$,
$T_{h}$ is a contraction mapping. Therefore, an application of the Brouwer
fixed-point theorem implies that there exists a unique fixed-point in
$W_{h}^{\rho}$. ∎

##  4 A priori error estimates

In this section, we derive a priori error estimates of the discrete solution.
In particular, the $L^{2}$ error estimates for all the variables are derived,
and a superconvergent result for velocity can be designed under
$\|\cdot\|_{h}$ norm.

Firstly, we have the following a priori bounds for
$\left\|\bm{u}_{h}-J_{h}\bm{u}\right\|_{h}$.

######  Lemma 4.1.

Let $(\bm{u},\bm{G})$ be the solution of (1) and $(\bm{u}_{h},\bm{G}_{h})\in
V_{h}\times H_{h}$ be the discrete solution of (7). Then the following
estimate holds:

$\left\|J_{h}\bm{u}-\bm{u}_{h}\right\|_{h}\leq
C_{B}\nu^{-1}\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}.$

###### Proof.

The Galerkin orthogonality reads

$\nu^{-1}(\bm{G}-\bm{G}_{h},\bm{H})+B_{h}^{*}(\bm{u}-\bm{u}_{h},\bm{H})=0\quad\forall\bm{H}\in
H_{h}.$

By Lemma 2.6, we have

$B_{h}(\bm{H},J_{h}\bm{u}-\bm{u}_{h})=-B_{h}^{*}(\bm{u}-\bm{u}_{h},\bm{H})=\nu^{-1}(\bm{G}-\bm{G}_{h},\bm{H}).$

This and Lemma 2.7 imply that

$\displaystyle\left\|J_{h}\bm{u}-\bm{u}_{h}\right\|_{h}$ $\displaystyle\leq
C_{B}\sup_{\bm{H}\in
H_{h}}\frac{B_{h}(\bm{H},J_{h}\bm{u}-\bm{u}_{h})}{\left\|\bm{H}\right\|_{L^{2}(\Omega)}}$
$\displaystyle=C_{B}\sup_{\bm{H}\in
H_{h}}\frac{\nu^{-1}(\bm{G}-\bm{G}_{h},\bm{H})}{\left\|\bm{H}\right\|_{L^{2}(\Omega)}}$
$\displaystyle\leq
C_{B}\nu^{-1}\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}.$

∎

We decompose the velocity error $\bm{u}-\bm{u}_{h}$ into two parts,

$\bm{u}-\bm{u}_{h}=(\bm{u}-J_{h}\bm{u})+(J_{h}\bm{u}-\bm{u}_{h})=:\bm{e}_{*}+\bm{e}_{h}.$

######  Lemma 4.2.

Let $\bm{f}=\bm{f}_{0}+\nabla\chi$ be the Helmholtz-Hodge decomposition of
$\bm{f}$. Then we have

$|N_{h}(\bm{u};\bm{u},\bm{v}_{h})-N_{h}(\bm{u}_{h};\bm{v}_{h},\bm{v}_{h})|\leq\nu^{-1}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}(C_{-1}\left\|\bm{e}_{*}\right\|_{0,4,h}+C_{*}\left\|\bm{e}_{h}\right\|_{h})\left\|\bm{v}_{h}\right\|_{h}$

for all $\bm{v}_{h}\in V_{h}$.

For the proof, see Appendix A. With the aid of the aforementioned lemmas, we
can prove the convergence of $\bm{G}_{h}$.

######  Theorem 4.3.

Assume that
$\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\leq\nu^{2}\min\\{\frac{1}{2}C_{*}^{-1}C_{B}^{2},C_{B}^{4}C_{N}^{-1}\\}$.
Then we have

$\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}\leq
C(\nu^{2}\left\|\bm{u}-J_{h}\bm{u}\right\|_{0,4,h}+\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}).$

###### Proof.

Let $(\bm{G}_{h},\bm{u}_{h},p_{h})\in H_{h}\times V_{h}\times Q_{h}^{0}$ be
the solution to (7). Then we have

$\displaystyle\nu^{-1}\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\nu^{-1}(\bm{G}_{h}-\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})$
$\displaystyle=B_{h}^{*}(\bm{u}_{h}-\bm{u},\bm{G}_{h}-\Pi_{h}\bm{G})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G}).$

By Lemma 2.6 and the adjoint properties of $B_{h}(\cdot,\cdot)$ and
$B_{h}^{*}(\cdot,\cdot)$, we obtain

$\displaystyle\nu^{-1}\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=B_{h}^{*}(\bm{u}_{h}-J_{h}\bm{u},\bm{G}_{h}-\Pi_{h}\bm{G})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})$
$\displaystyle=B_{h}(\bm{G}_{h}-\Pi_{h}\bm{G},\bm{u}_{h}-J_{h}\bm{u})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})$
$\displaystyle=B_{h}(\bm{G}_{h}-\bm{G},\bm{u}_{h}-J_{h}\bm{u})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G}).$

The consistency of discrete operators and (7b) yield

$\displaystyle\nu^{-1}\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=N_{h}(\bm{u};\bm{u},\bm{u}_{h}-J_{h}\bm{u})-N_{h}(\bm{u}_{h};\bm{u}_{h},\bm{u}_{h}-J_{h}\bm{u})$
$\displaystyle\quad+b_{h}^{*}(p-p_{h},\bm{u}_{h}-J_{h}\bm{u})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G}).$

Lemmas 2.6, 4.2, (13) and (7c) lead to

$\displaystyle\nu^{-1}\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=N_{h}(\bm{u};\bm{u},\bm{u}_{h}-J_{h}\bm{u})-N_{h}(\bm{u}_{h};\bm{u}_{h},\bm{u}_{h}-J_{h}\bm{u})$
$\displaystyle\quad+b_{h}(\bm{u}_{h}-J_{h}\bm{u},I_{h}p-p_{h})+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})$
$\displaystyle=N_{h}(\bm{u};\bm{u},\bm{u}_{h}-J_{h}\bm{u})-N_{h}(\bm{u}_{h};\bm{u}_{h},\bm{u}_{h}-J_{h}\bm{u})$
$\displaystyle\quad+\nu^{-1}(\bm{G}-\Pi_{h}\bm{G},\bm{G}_{h}-\Pi_{h}\bm{G})$
$\displaystyle\leq
C\nu^{-1}(\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|\bm{u}-J_{h}\bm{u}\right\|_{0,4,h}+\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)})\left\|\bm{G}_{h}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}$
$\displaystyle\quad+C_{*}\nu^{-1}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|J_{h}\bm{u}-\bm{u}_{h}\right\|_{h}^{2}.$

Using Lemma 4.1, we obtain

$\nu^{-1}(1-C_{*}C_{B}^{-2}\nu^{-2}\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)})\left\|\bm{G}_{h}-\Pi_{h}\bm{G}_{h}\right\|_{L^{2}(\Omega)}\leq
C\nu^{-1}(\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|\bm{u}-J_{h}\bm{u}\right\|_{0,4,h}+\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}).$

This coupled with the assumption
$\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\leq\frac{1}{2}C_{*}^{-1}C_{B}^{2}\nu^{2}$
yields

$\left\|\bm{G}_{h}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}\leq
C(\left\|\bm{f}_{0}\right\|_{L^{2}(\Omega)}\left\|\bm{u}-J_{h}\bm{u}\right\|_{0,4,h}+\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}).$

Therefore, the proof is completed. ∎

Now, we are ready to prove the main theorem of this paper.

######  Theorem 4.4.

Assume that $\bm{f}_{0}$ satisfies the assumption in Theorem 4.3 and the weak
solution $(\bm{u},p)$ to (1) belongs to $[H^{k+2}(\Omega)]^{2}\times
H^{k+1}(\Omega)$. Then it holds

$\displaystyle\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq C\nu h^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)},$
$\displaystyle\left\|\bm{u}-\bm{u}_{h}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq Ch^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)},$
$\displaystyle\left\|p-p_{h}\right\|_{L^{2}(\Omega)}$ $\displaystyle\leq
Ch^{k+1}(\nu\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}+\left\|p\right\|_{H^{k+1}(\Omega)}).$

In addition, the following superconvergence holds

$\displaystyle\left\|J_{h}\bm{u}-\bm{u}_{h}\right\|_{h}$ $\displaystyle\leq
Ch^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}.$

###### Proof.

From Theorem 4.4 and Lemma 2.5, we have

$\displaystyle\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}+\left\|\bm{G}_{h}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C(\nu^{2}\left\|\bm{u}-J_{h}\bm{u}\right\|_{0,4,h}+\left\|\bm{G}-\Pi_{h}\bm{G}\right\|_{L^{2}(\Omega)})$
$\displaystyle\leq
Ch^{k+1}(\nu^{2}\left\|\bm{u}\right\|_{W^{k+1,4}(\Omega)}+\left\|\bm{G}\right\|_{H^{k+1}(\Omega)}).$

Note that $\bm{G}=\nu\nabla\bm{u}$, the Sobolev embedding theorem yields

$\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}\leq
Ch^{k+1}(\nu^{2}\left\|\bm{u}\right\|_{W^{k+1,4}(\Omega)}+\nu\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}).$

Then the assumption on $\bm{f}_{0}$ implies

$\left\|\bm{G}-\bm{G}_{h}\right\|_{L^{2}(\Omega)}\leq C\nu
h^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}.$

This and Lemma 4.1 yield the superconvergence

$\left\|J_{h}\bm{u}-\bm{u}_{h}\right\|_{h}\leq
C\nu^{-1}\left\|\bm{G}-\bm{G}_{h}\right\|_{h}\leq
Ch^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}.$

The second assertion can be derived by using the triangle inequality and Lemma
2.3.

$\displaystyle\left\|\bm{u}-\bm{u}_{h}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq\left\|\bm{u}-J_{h}\bm{u}\right\|_{L^{2}(\Omega)}+\left\|\bm{u}_{h}-J_{h}\bm{u}\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq\left\|\bm{u}-J_{h}\bm{u}\right\|_{L^{2}(\Omega)}+C\left\|\bm{u}_{h}-J_{h}\bm{u}\right\|_{h}$
$\displaystyle\leq Ch^{k+1}\left\|\bm{u}\right\|_{H^{k+2}(\Omega)}.$

Finally, we prove the error estimates for the pressure variable. By the inf-
sup condition [3] and integration by parts, the following holds

$\displaystyle C\left\|p_{h}-p\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{(\nabla\cdot\bm{v},p_{h}-p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}$
$\displaystyle\leq\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{(\nabla\cdot\bm{v},p_{h}-I_{h}p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}+\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{(\nabla\cdot\bm{v},I_{h}p-p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}$
$\displaystyle\leq\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{-b_{h}(\bm{v},p_{h}-I_{h}p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}+Ch^{k+1}\left\|p\right\|_{H^{k+1}(\Omega)}.$

Therefore, it is enough to bound the first term on the right-hand side. Using
the definition and the stability estimate of $J_{h}$, we obtain

$\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{-b_{h}(\bm{v},p_{h}-I_{h}p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}\leq
C\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{-b_{h}(J_{h}\bm{v},p_{h}-I_{h}p)}{\left\|J_{h}\bm{v}\right\|_{h}}.$

The error equation reads

$\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{-b_{h}(\bm{v},p_{h}-I_{h}p)}{\left\|\bm{v}\right\|_{H^{1}(\Omega)}}\leq
C\sup_{\bm{v}\in
H^{1}_{0}(\Omega)}\frac{B_{h}(\bm{G}_{h}-\bm{G},J_{h}\bm{v})+N_{h}(\bm{u}_{h};\bm{u}_{h},J_{h}\bm{v})-N_{h}(\bm{u};\bm{u},J_{h}\bm{v})}{\left\|J_{h}\bm{v}\right\|_{h}}.$

Combining Lemma 4.2 and the above estimates completes the proof. ∎

######  Remark 4.5.

The usual finite element methods for the (Navier-)Stokes problem yields
velocity error estimates that depend on $\nu^{-1}\left\|p-I_{h}p\right\|_{0}$.
Thanks to the special finite element pairs designed, our method yields
divergence free velocity, thereby the velocity error estimate is independent
of $\nu$. The error estimates for velocity gradient and pressure also depend
on the positive power of $\nu$, which is desirable when $\nu\ll 1$.

##  5 Numerical Experiments

In this section, we verify our theoretical findings and demonstrate the
performance of the proposed method. In particular, the pressure-robustness and
optimality of the proposed method will be illustrated. To solve the nonlinear
equation, we used the fixed-point iteration: For $n=1,2,\cdots$, find
$(\bm{G}_{h}^{n},\bm{u}_{h}^{n},p_{h}^{n})\in H_{h}\times V_{h}\times Q_{h}$
such that

$\displaystyle\nu^{-1}(\bm{G}_{h}^{n},\bm{H})+B_{h}^{*}(\bm{u}_{h}^{n},\bm{H})$
$\displaystyle=\sum_{e\in\mathcal{F}_{pr}^{b}}\left\langle\bm{g},\bm{H}\mathbf{n}\right\rangle_{e}$
$\displaystyle\forall\bm{H}\in H_{h},$ $\displaystyle
B_{h}(\bm{G}_{h}^{n},\bm{v})+N_{h}(\bm{u}_{h}^{n-1};\bm{u}_{h}^{n},\bm{v})+b_{h}^{*}(p^{n},\bm{v})$
$\displaystyle=(\bm{f},\bm{v})+\sum_{e\in\mathcal{F}_{h}^{b}}\left\langle|\bm{g}\cdot\mathbf{n}|,\bm{g}\cdot\bm{v}\right\rangle_{e}$
$\displaystyle\forall\bm{v}\in V_{h},$ $\displaystyle b_{h}(\bm{u}_{h}^{n},q)$
$\displaystyle=\sum_{e\in\mathcal{F}_{pr}^{b}}\left\langle\bm{g}\cdot\mathbf{n},q\right\rangle_{e}$
$\displaystyle\forall q\in Q_{h}$

with known $\bm{u}_{h}^{n-1}$ and the initial guess is chosen to be zero. The
stopping criteria is set to a successive maximum error at nodal points with
the tolerance $10^{-7}$. In this section, $\|\cdot\|$ in the legend of each
figure represents $\|\cdot\|_{L^{2}(\Omega)}$.

###  5.1 Taylor vortex

Figure 2: Convergence history with respect to degrees of freedom with $k=1$
(left) and $k=2$ (right) for the Taylor vortex example.

Figure 3: $L^{2}$-error with various $\nu$ for the Taylor vortex example of
the Navier-Stokes equations. Here, $h=1/16$ is used with $k=1$ (left) and
$k=2$ (right).

Figure 4: $L^{2}$-error with various $\nu$ for the Taylor vortex example of
the Stokes equations. Here, $h=1/16$ is used with $k=1$ (left) and $k=2$
(right).

In this example, the optimal convergence and the pressure-robustness of the
solution are studied. Consider a solution $(\bm{u},p)$ defined by

$\bm{u}=\begin{bmatrix}-\cos(\pi x)\sin(\pi y)\\\ \phantom{-}\sin(\pi
x)\cos(\pi y)\end{bmatrix},\quad p=-\frac{\cos(2\pi x)+\cos(2\pi y)}{4}.$

The velocity gradient $\bm{G}$ and the force $\bm{f}$ can be computed from the
solution.

In the first experiment, we set $\nu=10^{-1}$. To observe the convergence
behavior of the discrete solution, uniform rectangular grids with
$h=2^{-2},\;2^{-3},\;\cdots,\;2^{-5}$ are considered. In Figure 2, the
convergence history against the number of degrees of freedom with polynomial
orders $k=1$ and $k=2$ is depicted. Optimal convergence rates derived in
Theorem 4.4 for all variables are observed for both $k=1$ and $k=2$.

The second experiment is performed to observe the pressure-robustness of the
proposed method. In Theorem 4.4, both velocity and gradient errors are
independent of pressure error. Also, we expect that the error is bounded by
interpolation error independent of $\nu$ when $\bm{f}$ satisfies the small
data assumption. In Figure 3, $L^{2}$-errors with varying $\nu$ are observed.
We can observe that the pressure error depends on $\nu$ when $\nu$ is not
sufficiently small which can be attributed to the dominant
$\nu\left\|\bm{u}-I_{h}\bm{u}\right\|_{L^{2}(\Omega)}$. However, the error
shows $\nu$-independent behavior for sufficiently small $\nu$, i.e., the
interpolation error of the pressure dominates the pressure error. For the
velocity error, we can observe that the velocity error is independent of $\nu$
when $\nu$ is sufficiently small. While the convergence analysis without small
data assumption for $\bm{f}$ implies that the error may depend on $\nu^{-1}$,
the error shows independent behavior for sufficiently small $\nu$.

The dependence of gradient error on $\nu^{-1}$ comes from the nonlinear
convective term $N_{h}(\cdot;\cdot,\cdot)$. Therefore, we can expect that our
formulation leads to $\nu$-independent error for velocity and velocity
gradient when the Stokes problem is considered. In Figure 4, we indeed have
$\nu$-independent error behavior for both $k=1$ and $k=2$.

###  5.2 No-flow example

Figure 5: Convergence history of $L^{2}$-norm for the no-flow example with
$\lambda=10^{7}$ and $\nu=1$. The formulation in [55] (left) and the current
formulation (right) are used, respectively.

Figure 6: No-flow example. Velocity magnitudes obtained with the formulation
given in [55] (left) and current (right). Both approximations are obtained
with a quasi-uniform polygonal mesh of $h\approx 1/8$.

In this example, we investigate the pressure-robustness more thoroughly.
Consider the following no-flow example

$\bm{u}=\bm{0},\quad p=\lambda(y^{3}-y^{2}/2+y-7/12),$

where $\lambda=10^{7}$. Since $\bm{u}\in V_{h}$, we have
$\left\|\bm{u}-J_{h}\bm{u}\right\|_{L^{2}(\Omega)}=0$. Therefore, one can
expect that a pressure-robust method produces a discrete velocity that
coincides with the exact velocity. In this example, we also implemented the
formulation derived in [55] for comparison. Since [55] considered the Stokes
equations with $k=0$, we extended the formulation proposed therein to
arbitrary order polynomial and handled the nonlinear term as in the current
formulation. In Figure 5, the convergence history of $L^{2}$-errors obtained
with quasi-uniform polygonal meshes and quadratic polynomial spaces is
depicted. While the discrete velocity obtained from [55] converges, its
magnitude is almost the same as the pressure error. This shows that the
velocity error depends on not only $(\bm{u}-J_{h}\bm{u})$ but also
$(p-I_{h}p)$. On the other hand, the discrete velocity produced by the
proposed method is around $10^{-10}$ regardless of the mesh size. Here,
$10^{-10}$ is compatible with the machine precision multiplied by the
condition number of the resulting matrix. We also included the velocity
magnitude profiles, in Figure 6, obtained from the two formulations when
$h\approx 1/8$. From those two observations, we can conclude that the proposed
method indeed leads to a pressure-robust velocity approximation.

###  5.3 Lid-driven cavity

Figure 7: Streamline of lid-driven cavity with $h=1/32$ and $k=2$.

In the last example, we consider the famous lid-driven cavity problem to
demonstrate the performance of the proposed method. Consider a rectangular
domain $(0,1)\times(0,1)$. We choose $\bm{f}=\bm{0}$ and the Dirichlet
boundary condition is given by

$\bm{u}=\begin{cases}[1,0]^{T}&\text{on }(0,1)\times\\{y=1\\},\\\
\bm{0}&\text{otherwise.}\end{cases}$

Piecewise quadratic polynomials on a uniform rectangular grid with $h=1/32$ is
used. In Figure 7, contour plots of streamfunctions for
$\nu^{-1}=400,\;1000,\;5000$ are displayed. Here, the contour levels are
chosen as in [30]. We can observe that the streamlines are qualitatively
similar to that of [30].

##  6 Conclusion

In this study, we developed a high-order polygonal staggered DG method for the
Navier-Stokes equations. A new finite element pair is introduced in the spirit
of [47], which enables us to achieve pressure robustness. It is worth
mentioning that the proposed scheme yields a divergence free velocity
approximation without any postprocessing which is a desirable feature in the
applications of (Navier-)Stokes equations. Another novel contribution lies in
the design a new nonlinear convective term that earns non-negativity. A full
convergence analysis is carried out, showing that the error estimate for
velocity is independent of pressure and of the viscosity. Finally, several
numerical experiments are carried out, and we can observe that our scheme is
indeed pressure robust. The numerical experiments indicate that our method is
a good candidate for practical applications in particular for problems with
high Reynolds number.

##  Appendix A Proof of Lemma 4.2

In this appendix, we prove Lemma 4.2.

###### Proof.

By adding and subtracting $N_{h}(J_{h}\bm{u};J_{h}\bm{u},\bm{v}_{h})$, we
obtain

$\displaystyle
N_{h}(\bm{u};\bm{u},\bm{v}_{h})-N_{h}(\bm{u}_{h};\bm{u}_{h},\bm{v}_{h})$
$\displaystyle=N_{h}(\bm{e}_{*};\bm{u},\bm{v}_{h})+N_{h}(J_{h}\bm{u};\bm{e}_{*},\bm{v}_{h})$
(17) $\displaystyle\quad
N_{h}(\bm{e}_{h};\bm{u}_{h},\bm{v}_{h})+N_{h}(J_{h}\bm{u};\bm{e}_{h},\bm{v}_{h}).$

We will bound each term separately. First we have

$\displaystyle N_{h}(\bm{e}_{*};\bm{u},\bm{v}_{h})$
$\displaystyle=-(\bm{u}\otimes\bm{e}_{*},\nabla\bm{v}_{h})+\left\langle\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\},\left\\{\\!\\!\left\\{\bm{u}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}_{h}\right\rrbracket\right\rangle_{\mathcal{F}_{h}}.$

The first term on the right-hand side can be bounded using Hölder’s inequality
and the Sobolev inequality

$-(\bm{u}\otimes\bm{e}_{*},\nabla\bm{v}_{h})\leq
C\left\|\bm{u}\right\|_{H^{1}(\Omega)}\left\|\bm{e}_{*}\right\|_{L^{4}(\mathcal{T}_{h})}\left\|\nabla\bm{v}_{h}\right\|_{L^{2}(\mathcal{T}_{h})}.$

Let us focus on the second term. For $e\in\mathcal{F}_{h}$, the trace
inequality (5) implies that

$\displaystyle\left\langle\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\},\left\\{\\!\\!\left\\{\bm{u}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}_{h}\right\rrbracket\right\rangle_{e}$
$\displaystyle\leq\left\|\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}\left\|\left\\{\\!\\!\left\\{\bm{u}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}\left\|\left\llbracket\bm{v}_{h}\right\rrbracket\right\|_{L^{2}(e)}$
$\displaystyle\leq
C\left(\left\|\bm{u}\right\|_{L^{4}(T)}+h_{T}^{1/2}\left\|\bm{u}\right\|_{H^{1}(T)}\right)h_{T}^{1/4}\left\|\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}h_{T}^{-1/2}\left\|\left\llbracket\bm{v}_{h}\right\rrbracket\right\|_{L^{2}(e)},$

where $T$ is a neighboring triangle of $e$. By using shape regularity
$h_{T}\leq\rho h_{e}$, summing over edges yields

$\displaystyle\left\langle\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\},\left\\{\\!\\!\left\\{\bm{u}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}_{h}\right\rrbracket\right\rangle_{\mathcal{F}_{h}}\leq
C\left\|\bm{u}\right\|_{L^{4}(\Omega)}\left(\sum_{e\in\mathcal{F}_{h}}h_{e}\left\|\left\\{\\!\\!\left\\{\bm{e}_{*}\cdot\mathbf{n}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}^{4}\right)^{1/4}\left(\sum_{e\in\mathcal{F}_{h}}h_{e}^{-1}\left\|\left\llbracket\bm{v}_{h}\right\rrbracket\right\|_{L^{2}(e)}^{2}\right)^{1/2}.$

By combining these, we obtain

$N_{h}(\bm{e}_{*};\bm{u},\bm{v}_{h})\leq
C\left\|\bm{u}\right\|_{H^{1}(\Omega)}\left\|\bm{e}_{*}\right\|_{0,4,h}\left\|\bm{v}_{h}\right\|_{h},$

where $\left\|\cdot\right\|_{0,4,h}$ is defined in (4). Similarly, we obtain

$N_{h}(J_{h}\bm{u};\bm{e}_{*},\bm{v}_{h})\leq
C\left\|J_{h}\bm{u}\right\|_{h}\left\|\bm{e}_{*}\right\|_{0,4,h}\left\|\bm{v}_{h}\right\|_{h}\leq
C\left\|\bm{u}\right\|_{H^{1}(\Omega)}\left\|\bm{e}_{*}\right\|_{0,4,h}\left\|\bm{v}_{h}\right\|.$

It remains to estimate the last two terms on the right hand side of (17) .
Observe that

$\displaystyle N_{h}(\bm{e}_{h};\bm{u}_{h},\bm{v}_{h})$
$\displaystyle=-(\bm{u}_{h}\otimes\bm{e}_{h},\nabla\bm{v}_{h})+\left\langle\left\\{\\!\\!\left\\{\bm{e}_{h}\cdot\mathbf{n}\right\\}\\!\\!\right\\},\left\\{\\!\\!\left\\{\bm{u}_{h}\right\\}\\!\\!\right\\}\cdot\left\llbracket\bm{v}_{h}\right\rrbracket\right\rangle_{\mathcal{F}_{h}}+\left\langle|\left\\{\\!\\!\left\\{\bm{e}_{h}\cdot\mathbf{n}\right\\}\\!\\!\right\\}|,\left\llbracket\bm{u}_{h}\right\rrbracket\cdot\left\llbracket\bm{v}_{h}\right\rrbracket\right\rangle_{\mathcal{F}_{h}}$
$\displaystyle\leq\left\|\bm{u}_{h}\right\|_{L^{4}(\Omega)}\left\|\bm{e}_{h}\right\|_{L^{4}(\Omega)}\left\|\nabla\bm{v}_{h}\right\|_{L^{2}(\mathcal{T}_{h})}$
$\displaystyle\quad+\sum_{e\in\mathcal{F}_{h}}\left(\left\|\left\\{\\!\\!\left\\{\bm{u}_{h}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}+\left\|\left\llbracket\bm{u}_{h}\right\rrbracket\right\|_{L^{4}(e)}\right)\left\|\left\\{\\!\\!\left\\{\bm{e}_{h}\cdot\mathbf{n}\right\\}\\!\\!\right\\}\right\|_{L^{4}(e)}\left\|\left\llbracket\bm{v}_{h}\right\rrbracket\right\|_{L^{2}(e)}.$

The discrete trace inequality (2.2) yields

$\displaystyle N_{h}(\bm{e}_{h};\bm{u}_{h},\bm{v}_{h})$
$\displaystyle\leq\left\|\bm{u}_{h}\right\|_{L^{4}(\Omega)}\left\|\bm{e}_{h}\right\|_{L^{4}(\Omega)}\left\|\nabla\bm{v}_{h}\right\|_{L^{2}(\mathcal{T}_{h})}$
$\displaystyle\quad+C\left\|\bm{e}_{h}\right\|_{L^{4}(\Omega)}\left\|\bm{u}_{h}\right\|_{L^{4}(\Omega)}\left(\sum_{e\in\mathcal{F}_{h}}h_{e}^{-1}\left\|\left\llbracket\bm{v}_{h}\right\rrbracket\right\|_{L^{2}(e)}^{2}\right)^{1/2}.$

By (6), we obtain

$N_{h}(\bm{e}_{h};\bm{u}_{h},\bm{v}_{h})\leq
C\left\|\bm{u}_{h}\right\|_{h}\left\|\bm{e}_{h}\right\|_{h}\left\|\bm{v}_{h}\right\|_{h}.$

Similarly, we obtain

$N_{h}(J_{h}\bm{u};\bm{e}_{h},\bm{v}_{h})\leq
C\left\|J_{h}\bm{u}\right\|_{h}\left\|\bm{e}_{h}\right\|_{h}\left\|\bm{v}_{h}\right\|_{h}\leq
C\left\|\bm{u}\right\|_{H^{1}(\Omega)}\left\|\bm{e}_{h}\right\|_{h}\left\|\bm{v}_{h}\right\|_{h}.$

Combining the preceding estimates with Lemma 2.1 and Theorem 3.4 concludes the
proof. ∎

## Acknowledgments

The research of Eric Chung is partially supported by the Hong Kong RGC General
Research Fund (Project numbers 14304719 and 14302018) and the CUHK Faculty of
Science Direct Grant 2020-21. The research of Eun-Jae Park was supported by
the National Research Foundation of Korea (NRF) grant funded by the Ministry
of Science and ICT (NRF-2015R1A5A1009350 and NRF-2019R1A2C2090021).

## References

  * [1] D. Arndt, H. Dallmann, and G. Lube, Local projection FEM stabilization for the time‐dependent incompressible Navier–Stokes problem, Numer. Methods Partial Differential Equations, (2015), pp. 1224–1250.
  * [2] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci., 23 (2013), pp. 199–214.
  * [3] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Elasticity Theory, Cambridge University Press, 2007.
  * [4] C. Brennecke, A. Linke, C. Merdon, and J. Schöberl, Optimal and pressure-independent $L^{2}$ velocity error estimates for a modified Crouzeix-Raviart Stokes element with BDM reconstructions, J. Comput. Math., 33 (2015), pp. 191–208.
  * [5] E. Burman and M. A. Fernández, Continuous interior penalty finite element method for the time-dependent Navier–Stokes equations: space discretization and convergence, Numer. Math., (2007), pp. 39–77.
  * [6] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston, hp-version Discontinuous Galerkin Methods on Polygonal and Polyhedral Meshes, Springer Briefs in Mathematics, Springer, Cham, 2017.
  * [7] A. Cangiani, E. H. Georgoulis, T. Pryer, and O. J. Sutton, A posteriori error estimates for the virtual element method, Numer. Math., 137 (2017), pp. 857–893.
  * [8] D. Castanon Quiroz and D. A. Di Pietro, A hybrid high-order method for the incompressible Navier–Stokes problem robust for large irrotational body forces, Comput. Math. Appl., 79 (2020), pp. 2655–2677.
  * [9] A. Cesmelioglu, B. Cockburn, and W. Qiu, Analysis of a hybridizable discontinuous Galerkin method for the steady-state incompressible Navier-Stokes equations, Math. Comp., 86 (2017), pp. 1643–1670.
  * [10] L. Chen and F. Wang, A divergence free weak virtual element method for the Stokes problem on polytopal meshes, Journal of Scientific Computing, (2019), pp. 864––886.
  * [11] S. W. Cheung, E. Chung, H. H. Kim, and Y. Qian, Staggered discontinuous Galerkin methods for the incompressible Navier-Stokes equations, J. Comput. Phys., 302 (2015), pp. 251–266.
  * [12] E. T. Chung and P. Ciarlet Jr., A staggered discontinuous Galerkin method for wave propagation in media with dielectrics and meta-materials, J. Comput. Appl. Math., 239:189–207 (2013).
  * [13] E. T. Chung, J. Du, and C. Y. Lam, Discontinuous Galerkin methods with staggered hybridization for linear elastodynamics, Comput. Math. Appl., 74 (2017), pp. 1198–1214.
  * [14] E. T. Chung and B. Engquist, Optimal discontinuous Galerkin methods for wave propagation, SIAM J. Numer. Anal., 44 (2006), pp. 2131–2158.
  * [15] E. T. Chung and B. Engquist, Optimal discontinuous Galerkin methods for the acoustic wave equation in higher dimensions, SIAM J. Numer. Anal., 47 (2009), pp. 3820–3848.
  * [16] E. T. Chung, H. H. Kim, and O. B. Widlund, Two-level overlapping Schwarz algorithms for a staggered discontinuous Galerkin method, SIAM J. Numer. Anal., 51 (2013), pp. 47–67.
  * [17] E. T. Chung and C. S. Lee, A staggered discontinuous Galerkin method for the curl–curl operator, IMA J. Numer. Anal., 32 (2011), pp. 1241–1265.
  * [18] E. T. Chung and W. Qiu, Analysis of an SDG method for the incompressible Navier–Stokes equations, SIAM J. Numer. Anal., 55 (2017), pp. 543–569.
  * [19] B. Cockburn, G. Kanschat, and D. Schötzau, A locally conservative LDG method for the incompressible Navier-Stokes equations, Math. Comp., 74 (2005), pp. 1067–1095.
  * [20] B. Cockburn, G. Kanschat, and D. Schötzau, A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes equations, J. Sci. Comput., 31 (2007), pp. 61–73.
  * [21] B. Cockburn, G. Kanschat, and D. Schötzau, An equal-order DG method for the incompressible Navier-Stokes equations, J. Sci. Comput., 40 (2009), pp. 188–210.
  * [22] D. A. Di Pietro and A. Ern, Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier–Stokes equations, Math. Comp., 79 (2010), pp. 1303–1330.
  * [23] D. A. Di Pietro and A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Springer-Verlag Berlin Heidelberg, 2012.
  * [24] D. A. Di Pietro and A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Meth. Appl. Mech. Engrg., (2015), pp. 1–21.
  * [25] O. Dorok, W. Grambow, and L. Tobiska, Aspects of Finite Element Discretizations for Solving the Boussinesq Approximation of the Navier-Stokes Equations, Vieweg+Teubner Verlag, Wiesbaden, 1994, pp. 50–61.
  * [26] L. P. Franca and S. L. Frey, Stabilized finite element methods: II. The incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg., 99 (1992), pp. 209–233.
  * [27] D. Frerichs and C. Merdon, Divergence-preserving reconstructions on polygons and a really pressure-robust virtual element method for the Stokes problem, IMA J. Numer. Anal., (2020), p. draa073.
  * [28] J. d. Frutos, B. García-Archilla, V. John, and J. Novo, Analysis of the grad-div stabilization for the time-dependent Navier–Stokes equations with inf-sup stable finite elements, Adv. Comput. Math., (2018), pp. 195–225.
  * [29] G. Fu, Y. Jin, and W. Qiu, Parameter-free superconvergent H(div)-conforming HDG methods for the Brinkman equations, IMA J. Numer. Anal., 39 (2018), pp. 957–982.
  * [30] U. Ghia, K. N. Ghia, and C. T. Shin, High-resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys., 48 (1982), pp. 387–411.
  * [31] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag Berlin Heidelberg, 1986.
  * [32] J. Guzmán and M. Neilan, Conforming and divergence-free Stokes elements in three dimensions, IMA J. Numer. Anal., 34 (2013), pp. 1489–1508.
  * [33] J. Guzmán and M. Neilan, Conforming and divergence-free Stokes elements on general triangular meshes, Math. Comp., 83 (2014), pp. 15–36.
  * [34] Y. Han and Y. Hou, Semirobust analysis of an H(div)-conforming DG method with semi-implicit time-marching for the evolutionary incompressible Navier–Stokes equations, IMA J. Numer. Anal., (2021), pp. 1–30.
  * [35] O. A. Karakashian and W. N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM J. Numer. Anal., 35 (1998), pp. 93–120.
  * [36] D. Kim, L. Zhao, and E.-J. Park, Staggered DG methods for the pseudostress-velocity formulation of the Stokes equations on general meshes, SIAM J. Sci. Comput., 42 (2020), pp. A2537–A2560.
  * [37] H. H. Kim, E. T. Chung, and C. S. Lee, A staggered discontinuous Galerkin method for the Stokes system, SIAM J. Numer. Anal., 51 (2013), pp. 3327–3350.
  * [38] A. Linke, G. Matthies, and L. Tobiska, Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors, ESIAM: Math. Model. Numer. Anal., 50 (2016), pp. 289–309.
  * [39] L. Mu, Pressure robust weak Galerkin finite element methods for Stokes problems, SIAM J. Sci. Comput., (2020), pp. B608–B629.
  * [40] S. Rhebergen and G. N. Wells, A hybridizable discontinuous Galerkin method for the Navier–Stokes equations with pointwise divergence-free velocity field, J. Sci. Comput., 76 (2018), pp. 1484–1501.
  * [41] B. Riviere and S. Sardar, Penalty-free discontinuous Galerkin methods for incompressible Navier–Stokes equations, Math. Models Methods Appl., 24 (2014), pp. 1217–1236.
  * [42] P. W. Schroeder, C. Lehrenfeld, A. Linke, and G. Lube, Towards computable flows and robust estimates for inf-sup stable FEM applied to the time-dependent incompressible Navier-Stokes equations, SeMA J., (2018), pp. 629–653.
  * [43] C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Comput. Fluids, 1 (1973), pp. 73–100.
  * [44] J. Wang, Y. Wang, and X. Ye, A robust numerical method for Stokes equations based on divergence-free H(div) finite element methods, SIAM J. Sci. Comput., 31 (2009), pp. 2784–2802.
  * [45] J. Wang and X. Ye, New finite element methods in computational fluid dynamics by H(div) elements, SIAM J. Numer. Anal., 45 (2007), pp. 1269–1286.
  * [46] J. Wang and X. Ye, A weak Galerkin mixed finite element method for second order elliptic problems, Math. Comp., 241 (2013), pp. 103–115.
  * [47] L. Zhao, E. Chung, and M. F. Lam, A new staggered DG method for the Brinkman problem robust in the Darcy and Stokes limits, Comput. Methods Appl. Mech. Engrg., 364 (2020), p. 112986.
  * [48] L. Zhao, E. Chung, E.-J. Park, and G. Zhou, Staggered DG method for coupling of the Stokes and Darcy–Forchheimer problems, SIAM J. Numer. Anal., 59 (2021), pp. 1–31.
  * [49] L. Zhao and E.-J. Park, A staggered discontinuous Galerkin method of minimal dimension on quadrilateral and polygonal meshes, SIAM J. Sci. Comput., 40 (2018), pp. A2543–A2567.
  * [50] L. Zhao and E.-J. Park, A lowest-order staggered DG method for the coupled Stokes–Darcy problem, IMA J. Numer. Anal., 40 (2020), pp. 2871–2897.
  * [51] L. Zhao and E.-J. Park, A new hybrid staggered discontinuous Galerkin method on general meshes, J. Sci. Comput., 82 (2020), pp. 1–33.
  * [52] L. Zhao and E.-J. Park, A staggered cell-centered DG method for linear elasticity on polygonal meshes, SIAM J. on Sci. Comput., 42 (2020), pp. A2158–A2181.
  * [53] L. Zhao, E.-J. Park, and E. Chung, A pressure robust staggered discontinuous Galerkin method for the Stokes equations, arXiv preprint, (2020), https://doi.org/arXiv:2007.00298.
  * [54] L. Zhao, E.-J. Park, and E. Chung, Staggered discontinuous Galerkin methods for the Helmholtz equation with large wave number, Comput. Math. Appl., 80 (2020), pp. 2676–2690.
  * [55] L. Zhao, E.-J. Park, and D. Shin, A staggered DG method of minimal dimension for the Stokes equations on general meshes, Comput. Methods Appl. Mech. Engrg., 345 (2019), pp. 854–875.