A new approach, iteration projection method, is proposed to solve the saddle point problem obtained after the full discretization of the unsteady Navier-Stokes equations. It is motivated by the Uzawa method, an iterative method for saddle point problems, and the convectional projection method for the Navier-Stokes equations that projects the intermediate velocity to the divergence free space only once per time step. The proposed method iterates projections in each time step with a proper convection form. We prove that the projection iterations converge with a certain parameter regime.  Optimal iteration convergence can be achieved with the modulation of parameter values.

This new method has several significant improvements over the Uzawa method and the projection method both theoretically and practically. First, when the iterative projections are fully convergent in each time step, the numerical velocity is weakly divergence free (pointwise divergence free in divergence free finite element methods), and  the stability and error estimate are rigorously proven.  With proper parameters, this method converges much faster than the Uzawa method. Second, numerical simulations show that with rather relaxed stopping criteria which require only a few iterations each time step, the numerical solution preserves stability and accuracy for high Reynolds numbers, where the convectional projection method would fail. Furthermore, this method retains the efficiency of the traditional projection method by decoupling the velocity and pressure fields, which splits the saddle point system into small elliptic problems. Three dimensional simulations with Taylor-Hood P2/P1 finite elements are presented to demonstrate the performance and efficiency of this method.

More importantly, this method is a generic approach and thus there are many potential improvements and extensions of the iterative projection method, including utilization of various convection forms, association with stabilization techniques for high Reynolds numbers, applications to other saddle point problems.

Keywords.  Navier-Stokes equations, iterative projection method, saddle-point problems,  divergence free velocity,  skew-symmetric convection, fully implicit scheme

§ INTRODUCTION

§.§ Problem and motivation

The incompressible Navier-Stokes equations are predominant in viscous fluid dynamics, including problems with high Reynolds numbers such as turbulence flows. However, it is extremely challenging to design efficient and robust numerical schemes for these problems. A major difficulty is the treatment of the incompressibility constraint,
or the divergence free condition, which couples the velocity and pressure.
Consider a bounded three-dimensional domain $\Omega$ and the Navier-Stokes equations

\begin{align}
{u}_t - \nu \Delta u + \text{NL}(u,u)  + \nabla p &= f,
\label{ns001}
\\
\nabla\cdot{u} &= 0,
\label{ns002}
\end{align}

with the nonlinear convection $\text{NL}(u,u)=(u\cdot \nabla)u$,  the Dirichlet boundary condition
$u|_{\partial\Omega} = 0$,
where $u=(u_1, u_2, u_3)$ is the velocity,  $p$ is  the kinematic pressure,  $\nu$ is the kinematic viscosity, and $f$ is the force. There are various equivalent forms of the convection term (e.g., see [9]), and this work focuses on the skew-symmetric form

\begin{align}
\text{NL}(u,v) = (u\cdot\nabla) v + \frac12 (\nabla\cdot u) v.
\end{align}

It is well-known that this form can be used to prove the stability of finite element methods for Navier-Stokes equations (see e.g., <cit.>).
This work adopts the popular BDF2 (Backward Differentiation Formula 2) time integration for Navier-Stokes equations.
Let $k>0$ be the uniform time step size. Denote the numerical solutions of  velocity and pressure at time step $t_n=n k$ as ${u}^n$ and $p^n$, respectively. One obtains the following scheme,

\begin{align}
(1.5-k\nu \Delta ) u^{n+1}
+ k \text{NL}(\tilde{w}, u^{n+1})
+ k\nabla p^{n+1} &= F^{n+1},
\label{eqn01}\\
\nabla\cdot u^{n+1} &= 0,\label{eqn02}
\end{align}

where $F^{n+1}=2{ u}^n-0.5{ u}^{n-1} + k f^{n+1}$.
The choice of $\tilde{w}$ leads to two different problems.

When $\tilde{w}=w^{n+1}\triangleq 2{u}^{n}-{u}^{n-1}$, a second order extrapolation of $u^{n+1}$, the convection term is semi-implicit. Then, a spatial discretization of this scheme using, e.g.,  the mixed finite elements methods yields a generalized saddle point problem (see Remark <ref>).
Some extensive studies of the solvers of the saddle point problems can be found in  [5, 1].

When $\tilde{w}=u^{n+1}$, the convection term is fully implicit and this scheme is nonlinear. Typically, a nonlinear solver such as the Picard or Newton method is employed to linearize it to a series of saddle point problems, which are then computed by saddle point solvers as mentioned above.

One oldest but still popular method for the saddle-point problem is the  Uzawa method [4] (see also [13, 29]) published in 1958. For   the problem (<ref>),(<ref>), the Uzawa method can be written as, for $s=0,1,\cdots$,

\begin{align}
(1.5 - k\nu) u^{n+1,s} + k \text{NL}(\tilde{w}^s, u^{n+1,s})
&= - k\nabla p^{n+1,s} +  F^{n+1},
\label{uzawa01}\\
p^{n+1,s+1} &= p^{n+1,s} -\rho  \nabla\cdot u^{n+1,s},  \label{uzawa02}
\end{align}

with the initial guess $p^{n+1,0}=2p^n-p^{n-1}$and
the boundary condition $u^{n+1,s}|_{\partial\Omega}=0$.
When $\tilde{w}^s=w^{n+1}$, this scheme corresponds to the semi-implicit convection term.
When $\tilde{w}^s = u^{n+1,s-1}$ with $u^{n+1,-1}=w^{n+1}$, this scheme relates to the fully implicit convection.
Here, the parameter $\rho$ is chosen in the range of $0<\rho<2\nu$ in order to make the iterations convergent, as discussed in <cit.>.

The projection method is a classic and popular method to solve the Navier-Stokes equations first developed in [10, 28] in 1968, where a systematic review can be found in [17]. Although it is not particularly targeting the saddle point problem  (<ref>) and (<ref>), it can be regarded as a one-step approximation of it. The crucial operation of this method is projecting the intermediate velocity, which is not divergence free, to the divergence free space in a certain sense.
One of the most recommended projection scheme in [17],
rotational incremental pressure-correction method,
can be put as follows,

\begin{align}
(1.5-k\nu) u^{n+1,*} + k \text{NL}(w^{n+1}, u^{n+1,*})  &=-k\nabla \tilde{p}^{n+1}+ F^{n+1},
\label{pp001}\\
-\Delta \phi^{n+1} &= -\frac{1}{k} \nabla\cdot u^{n+1,*},
\label{pp002}\\
p^{n+1} &= \tilde{p}^{n+1}  +1.5\phi^{n+1} - \nu  \nabla\cdot { u}^{n+1,*},
\label{pp003}\\
u^{n+1} &= u^{n+1,*} - k \nabla \phi^{n+1},
\label{pp004}
\end{align}

where $\tilde{p}^{n+1}=2p^n-p^{n-1}$,
$\frac{\partial\phi^{n+1}}{\partial{ n}}\Big|_{\partial\Omega} = 0$.
The variable $\phi^{n+1}$ is from the Hodge decomposition $u^{n+1,*}= u^{n+1} + k \nabla \phi^{n+1}$, where $u^{n+1}$ is supposed to be divergence free.
This method was first proposed in [30] in 1996 and a similar one was  introduced in [7].
It was proven in [19] that the velocity in $H^1$ norm and pressure in $L^2$ norm have convergence order $3/2$ in time. For simplicity, we call this method “rotational projection method” in this paper.

Another popular projection scheme replaces the pressure update  (<ref>) with

\begin{equation}
p^{n+1} = \tilde{p}^{n+1}  +1.5\phi^{n+1}  \label{standard_proj}
\end{equation}

and is called the standard incremental pressure-correction method in [17]. We simply call it “standard projection method” in this work.  This one is also widely used in literature, such as [16] and [2].
As shown in [17], this scheme produces first order convergence in time of the velocity in $H^1$ norm and the pressure in $L^2$ norm.

The main strength of the projection method lies in its simplicity and efficiency. It decouples the velocity and pressure fields in the Navier-Stokes equations and splits the original system into two smaller standard problems: one advection-diffusion equation for $u^{n+1,*}$ and one Poisson equation for $\phi^{n+1}$, followed by two updates to get the end-of-step velocity and pressure. Therefore, its computational cost is far less than that of any iterative methods including the Uzawa and Newton methods.
However,  it has some serious defects. First, the velocity field obtained is not divergence free even in the weak sense in the mixed finite element implementations (see [18] and Appendix <ref>), and even when the finite element method admits divergence free velocity field (see Definition <ref> and  [26]). This would induce mass loss when the
velocity is used to transport a density function.
Second, the stability proof for the BDF2 schemes is not available, largely because of the complexity of this multi-step approach.
Third, this scheme is unable to treat high Reynolds numbers due to the violation of the divergence free condition  (see Section <ref> and <ref>).  In contrast, the Uzawa method can achieve weakly divergence free velocity solution in the mixed finite element implementations (see e.g. [13, 29] and Section <ref>).

§.§ New approach: iterative projection method

Motivated by the Uzawa iterations, we are interested in  whether the iterative projections could efficiently solve the Navier-Stokes equations and its impact on the incompressibility, stability, and error estimate of the Navier-Stokes solutions, especially when the Reynolds number is high.
As for the saddle point problem (<ref>) and (<ref>), we propose the following iterative projection method, with the iteration index $s=0,1,\cdots$,

\begin{align}
(1.5 -k\nu \Delta) u^{n+1,s}
+ k \text{NL}(\tilde{w}^{s}, u^{n+1,s})
& = - k\nabla p^{n+1,s} + F^{n+1},
\label{mol001}\\
-\Delta \phi^{n+1,s} &= -\frac{1}{k} \nabla\cdot u^{n+1,s},
\label{mol002}\\
p^{n+1,s+1} &= {p}^{n+1,s} + \alpha \phi^{n+1,s} - \rho  \nabla\cdot u^{n+1,s},
\label{mol003}
\end{align}

with boundary conditions $u^{n+1,s}|_{\partial\Omega}=0$ and
$\frac{\partial\phi^{n+1,s}}{\partial{ n}}\Big|_{\partial\Omega} = 0$.
When the iterations converge, the numerical solution of the Navier-Stokes equations at time step $n+1$ is defined as $(u^{n+1}, p^{n+1})=\lim_{s\to\infty}(u^{n+1,s}, p^{n+1,s})$.
The choice of $\tilde{w}^s$ results in different schemes.

When $\tilde{w}^n=w^{n+1}$,  this scheme is called Semi-Implicit Skew-Symmetric (SI-SKEW) Iterative Projection Method, because it  converges to the semi-implicit version of (<ref>) and (<ref>).

When $\tilde{w}^s = u^{n+1,s-1}$, it is named Fully Implicit Skew-Symmetric (FI-SKEW) Iterative Projection Method,  because it is derived from the fully implicit version of (<ref>) and (<ref>).

The control parameters $\alpha$ and $\rho$ are  used to optimize the convergence speed. All the schemes in Section <ref> are special cases indicated as below.

(1) When $\alpha=1.5$ and $\rho=\nu$, it is the iteration of the rotational projection scheme.

(2) When $\alpha=1.5$ and $\rho=0$, it is the iteration of the standard projection scheme.

(3) When $\alpha=0$, it is the Uzawa iteration.

It is important to point out that the proposed iterative projection method is a new solver of the saddle point problem (<ref>) and (<ref>), different than any existing approaches (e.g, [5, 1]).
Indeed, this strategy, the iterative projections to the divergence free space,  has been utilized to solve the Navier-Stokes equations in some previous work including  [31, 12, 3].
In [31], it was found that the repeated projections can reduce the spurious errors generated from the singular surface tension forces in the free boundary problems.
In [12],  the iterative projections provide accurate velocity fields in the stratification of temperature in the Boussinesq flows [12].
In [3], the projection  is iterated along with implicit-explicit (IMEX) time-integration solvers, where it is found that iterations could improve splitting errors of projection method.  All these three papers use the fully  explicit form of the convection. In [3], the authors proved the iterative convergence at each time step when the parameter $\alpha$ is fixed according to the time discretization schemes and $\rho$ is fixed as $\nu$. So far, there have been no study on the iterative projection schemes with semi-implicit or fully implicit convection terms, nor the theoretical results on the stability and error estimates of this strategy, which are  the targets of this work.

A closely related fact is that the conventional projection method has been used as preconditioners for Krylov subspace methods, such as the GMRES (generalized minimal residual) algorithms. For example, this approach is used in [15, 8] to solve velocity and pressure simultaneously in a system similar to (<ref>) and  (<ref>). In contrast, this work uses the iterative projection as a standalone solver.

In Section <ref>, we analyze the convergence of the proposed iterative projection methods for the saddle point problem (<ref>) and (<ref>), using normal modes and mixed finite element methods.
In Section,<ref>, we prove the stability and provide the error estimates of the proposed methods in the full temporal and spatial discretizations under the condition that the iterative projections are convergent in each time step.
In Section <ref>, we present the numerical simulations based on the Taylor-Hood P2/P1 finite elements to test the iterative projection method with high Reynolds numbers. In particular, Section <ref> introduces two incompressibilfity measures to evaluate the divergence free condition, Section <ref>  tests the iterative projection convergence at one single time step, Section <ref> discusses the optimal parameters values for iteration convergence, Section <ref> examines the performance of the proposed method of treating high Reynolds number in Problem 1 with the prescribed solution, and Section <ref> investigates the method on a classic lid driven cavity flow problem.  The conclusions and discussions are given in Section <ref>.

§ CONVERGENCE OF ITERATIVE PROJECTION METHOD FOR THE SADDLE POINT PROBLEM

This section is devoted to the iteration analysis of the proposed iterative projection method (<ref>), (<ref>), (<ref>) for the saddle point problem
(<ref>) and (<ref>), including both the semi-implicit and fully implicit convection forms. To simplify notations, we delete the time step superscript $(n+1)$  in this section.

§.§ Normal mode analysis without convection

In order to perform the normal mode analysis, we neglect the convection term  and assume the solution is smooth and periodic in the region $\Omega=[0,2\pi]^3$.

Note that the intermediate value $\phi^s$ in (<ref>) can be written as
\phi^s = -\frac{1}{k} (-\Delta)^{-1} (\nabla\cdot { u}^s),
where $(-\Delta)^{-1}$ refers to the inverse operator of the  Neumann problem (<ref>).
So the iterative scheme (<ref>), (<ref>), (<ref>) when the convection is removed can be simplified to, for $s=0, 1, \cdots$,

\begin{align}
1.5 { u}^{s}
- k\nu \Delta { u}^s + k\nabla p^{s} &= F,
\hspace{.3 in} { u}^s|_{\partial\Omega}=0,
\label{weiwei001}\\
p^{s+1} &= {p}^{s} -\big(\rho +  \frac{\alpha}{k} (-\Delta)^{-1} \big)  (\nabla\cdot { u}^{s}).
\label{weiwei002}
\end{align}

Consider the solution $({u},p)$ satisfying (<ref>) and (<ref>) without the convection, and the sequence $({u}^s,p^s)$ of (<ref>) and (<ref>).
Let $\bar{u}^s={u}^s-{u}$ and $\bar{p}^s=p^s-p$. Then

\begin{align}
1.5 \bar{{u}}^s - k\nu \Delta \bar{u}^s + k \nabla\bar{p}^s &=0,
\label{pup01}\\
\bar{p}^{s+1} &=\bar{p}^{s} - \big(\rho + \frac{\alpha}{k} (-\Delta)^{-1}\big) (\nabla\cdot \bar{u}^{s}). \label{pup02}
\end{align}

Denote $\bar{u}^s=(\bar{u}^s_1, \bar{u}^s_2, \bar{u}^s_3)$, the multi-wavenumber $\xi=(\xi_1, \xi_2, \xi_3)\in\mathbb{Z}^3$,
$\xi\cdot x=\xi_1 x_1 + \xi_2 x_2 + \xi_3 x_3$.
Denote the Fourier series of $\bar{u}^s_j$ and $\bar{p}^s$ as
$\bar{u}_j^s(x) = \sum_{\xi\in \mathbb{Z}^3} \widehat{\bar{u}_j^s}(\xi) e^{i\xi\cdot x}$, $j=1,2,3$,
$\bar{p}^s(x) = \sum_{\xi\in\mathbb{Z}^3} \widehat{\bar{p}^s}(\xi) e^{i\xi\cdot x}$, where $i=\sqrt{-1}$.
Then (<ref>) and (<ref>) become

\begin{align}
1.5 \widehat{\bar{u}_j^s}(\xi) + k\nu |\xi|^2 \widehat{\bar{u}_j^s}(\xi)
+ ki\xi_j \widehat{\bar{p}^s}(\xi)&=0,\hspace{.3 in} j=1,2,3;
\label{mop001}
\\
\widehat{\bar{p}^{s+1}}(\xi) &=
\widehat{\bar{p}^{s}}(\xi) -
\Big(
\rho + \frac{\alpha}{k|\xi|^2} \Big)
\sum_{j=1}^3 i\xi_j \widehat{\bar{u}^s_j}(\xi).
\label{mop002}
\end{align}

It can be calculated from (<ref>) that
\widehat{\bar{u}^s_j}(\xi) = \frac{-ki\xi_j \widehat{\bar{p}^s}}{1.5+k\nu |\xi|^2}
$, which is substituted into (<ref>) to achieve

\begin{equation}
\widehat{\bar{p}^{s+1}}(\xi) =
\frac{ (1.5-\alpha) + k(\nu-\rho)|\xi|^2}{1.5+k\nu|\xi|^2}
\cdot \widehat{\bar{p}^{s}}(\xi)
\triangleq C(\alpha,\rho)  \cdot \widehat{\bar{p}^{s}}(\xi).
\end{equation}

For the sequence $\widehat{\bar{p}^{s}}$ to converge to zero, it is equivalent to set $|C(\alpha,\rho)|<1$.
Note the optimal convergence occurs when the constant $C=0$, that is, $\alpha=1.5$ and $\rho=\nu$, which corresponds to the iterative rotational projection method.  Therefore, this method obtains the exact solution in just one iteration. However, this analysis is based on the smooth solutions and is without convection, which is not true in the full Navier-Stokes equations with the convection term and  spatial discretizations of  limited solution regularity.

If $\alpha=1.5$, then $|C(\alpha,\rho)|<1$ leads to

\begin{equation}
-\frac{1.5}{k |\xi|^2} < \rho < 2\nu + \frac{1.5}{k|\xi|^2}.
\end{equation}

For these inequalities to hold for all $\xi\in\mathbb{Z}^3$, it is equivalent to have $0\le \rho\le 2\nu$.

If $\alpha=1.5$ and $\rho=0$, that is, when the iterative standard projection method is used, then

\begin{equation}
C(1.5, 0)=\frac{k\nu|\xi|^2}{1.5+k\nu|\xi|^2}.
\end{equation}

In this case, the iteration always converges but the convergent constant $C$ is close to 0 when $k\nu|\xi|^2$ is small, and approaching 1 when $k\nu|\xi|^2$ is large. Thus, this iterative scheme is preferred only when $\nu$ is sufficiently small when the time step $k$ and largest wavenumber magnitude $|\xi|$ are fixed.

If $\alpha=0$, then this scheme reduces to the Uzawa scheme and the convergence constant is

\begin{equation}
\frac{1.5 + k(\nu-\rho) |\xi|^2}{1.5+k\nu|\xi|^2}.
\end{equation}

The convergence requires
$ 0 < \rho \le 2\nu$.
Let $z=k\nu |\xi|^2$, then

\begin{equation}
C_u(z;\rho) = \frac{1.5+ z (1-\frac{\rho}{\nu})}{1.5+z}.
\end{equation}

Its graph is shown in Figure <ref>.
When $z$ is large, the best choice to obtain fast convergence is $\rho=\nu$. When $z$ is small, the best choice of $\rho$ should be closer to $2\nu$. But when $z\ll 1.5$, $|Cu|$ is very close to 1 regardless the choice of $\rho$ in the range $(0,2\nu]$. Thus, the Uzawa iterations would be applicable only when $\nu$ is sufficiently large. This is consistent with the fact that the Uzawa method is mainly used for the Stokes equations whose viscosity is considerably large.

Uzawa method: the convergence constant $C_u(z;\rho)$ as a function of $z=k\nu|\xi|^2$ with parameter $\rho\in (0,2\nu]$.

The above analysis is summarized in Table <ref>. The rotational projection iteration method takes the advantages of the other two methods: the fast convergence of standard projection for small viscosity values and the Uzawa method for large viscosity values. Therefore, the rotational projection iteration is comparatively indifferent to viscosity variations.

Normal mode analysis when the convection is neglected. The  pressure update is $p^{s+1}={p}^{s} -(\rho +  \frac{\alpha}{k} (-\Delta)^{-1} )  (\nabla\cdot { u}^{s})$.

method     pressure iteration      Convergence constant     parameter range    notes

rotational projection     $\alpha=1.5$, $\rho=\nu$

standard projection     $\alpha=1.5$, $\rho=0$
Convergence is fast when $k\nu|\xi|^2$ is small

Uzawa     $\alpha=0$
$\frac{1.5 + k(\nu-\rho) |\xi|^2}{1.5+k\nu|\xi|^2}
$0<\rho\le 2\nu$     Convergence is fast when $\rho=\nu$ and $k\nu|\xi|^2$ is large.

§.§ Iteration analysis with mixed finite element method

If there is no convection or the convection term is treated utterly explicitly, the iteration convergence  for a special case when $\rho=\nu$ has been proven in [3], whose proof utilizes the eigenvalue analysis in [8]. Here, we present a self-contained  convergence study for the general case with parameters $\alpha$ and $\rho$ in Section <ref>.
More efforts are spent on the semi-implicit and fully implicit convection cases, as shown in Section <ref> and Section <ref>, respectively.

We denote $L^2(\Omega)$ as the Lebesgue space of square integrable functions with inner product $(u,v) = \int_{\Omega} u(x) v(x) \mathrm{d}x$,
and the $L^2$ norm as $\|u\|_{L^2(\Omega)} = \sqrt{(u,u)}$. Denote
$H^1_0(\Omega) = \{f: f, \frac{\partial f}{\partial x_i}\in L^2(\Omega), i=1,2, 3, f|_{\partial\Omega}=0 \}$. We introduce

\begin{align}
V&= \{u=(u_1, u_2, u_3): u_i \in H^1_0(\Omega), i=1,2,3 \} ,
\label{def_V}\\
Q&=L^2(\Omega)/\mathbb{R} =\{q\in L^2(\Omega): \frac{1}{|\Omega|} \int_\Omega q(x) \mathrm{d} x=0 \},
\label{def_Q}
\end{align}

with the norm $|u|_1 = (\int_{\Omega} |\nabla u(x)|^2 \mathrm{d}x)^{1/2}$,
where $|\nabla u|^2=\sum_{i,j=1}^3 |\frac{\partial u_i}{\partial x_j}|^2
$ and $\frac{1}{|\Omega|} \int_\Omega q(x) \mathrm{d}x $ is the average of $q$ over $\Omega$. Let $V_h\subset V$ and $Q_h\subset Q$ be a pair of conforming mixed finite element spaces  approximating  the velocity and pressure, respectively,
which satisfies the inf-sup condition (a.k.a. the Ladyzhenskaya-Babuka-Brezzi (LBB) condition, see e.g. [13]), i.e.,

\begin{equation}
\inf_{q\in Q_h} \sup_{v\in V_h}
\frac{(\nabla\cdot v, q)}{\|q\|_{L^2(\Omega)/\mathbb{R}} |v|_1} > 0.
\label{inf_sup}
\end{equation}

Denote the basis of $V_h$ as $w_i, i=1, \cdots, l$, and the basis of $Q_h$ as $q_j, j=1,\cdots, m$.

A velocity field $u_h\in V_h$ is weakly divergence free if $(\nabla\cdot u_h, q_h)=0$ for all $q_h\in Q_h$ and strongly divergence free if $\nabla\cdot u_h$ is almost everywhere zero in $\Omega$.
A pair of mixed finite element spaces $V_h$ and $Q_h$ is called  divergence free if $\nabla\cdot V_h\subset Q_h$, where $\nabla\cdot V_h$ is the range of the divergence of $V_h$. Furthermore, a finite element method with a pair of divergence free finite element spaces is called a div-free FEM, and otherwise a non-div-free FEM.

In this work, we use a weak incompressibility measure (<ref>) (a strong incompressiblity measure (<ref>))  to check whether a velocity is weakly (strongly) divergence free. A weakly divergence free velocity may be not pointwise divergence free, as seen in Section <ref>.
However, under a pair of divergence free FEM spaces, a weakly divergence free velocity is also strongly divergence free. Indeed, if $u_h$ satisfies $(\nabla\cdot u_h, q_h)=0$ for all $q_h\in Q_h$ and $\nabla\cdot V_h\subset Q_h$, then $\|\nabla\cdot u_h\|_{L^2(\Omega)}=0$ by choosing $q_h=\nabla\cdot u_h$, which suggests
$\nabla\cdot u_h$ is zero almost everywhere in $\Omega$.
One popular example of divergence free elements is Scott–Vogelius pair [27]. However, many $H^1$-conforming, inf-sup stable mixed finite elements pairs are not  divergence free, including the very popular Tayor-Hood elements.

§.§.§ Iteration analysis without convection

Without the convection term, the weak form of the problem (<ref>), (<ref>) and
(<ref>) in the mixed finite element method is finding
${ u}^s\in V_h$,
$\phi^s\in Q_h$,
$p^s\in L^2(\Omega)$,  such that

\begin{align}
1.5({u^s}, {w_i}) + k\nu (\nabla {u^s}, \nabla{w_i})
%+ k(({w}\cdot\nabla) {u^s}, {w_i})
- k(\nabla\cdot {w_i}, p^s) = (F,{w_i}),
&\qquad\forall\, i=1,\cdots, l,
\label{ip001}\\
(\nabla \phi^s, \nabla q_i) =
-\frac{1}{k}(\nabla \cdot u^s, q_i),  &\qquad\forall\, i=1,\cdots, m,
\label{ip002} \\
(p^{s+1}, q_i) = (p^s, q_i) + \alpha (\phi^s, q_i) - \rho (\nabla\cdot{ u}^s, q_i), &\qquad\forall\, i=1, \cdots, m.
\label{ip003}
\end{align}

The matrices $A_{l\times l}$, $B_{m\times l}$,
$G_{m\times m}$, and $M_{m\times m}$ are introduced as follows,

\begin{align}
A_{ij} &= 1.5 (w_i, w_j) + k\nu (\nabla w_i, \nabla w_j),
%+ k (w_i,(w\cdot \nabla)w_j),
\qquad \forall\, i,j=1,\cdots, l;
\label{matrixA} \\
B_{ij} &= - (q_i, \nabla\cdot w_j), \qquad \forall\, i=1,\cdots, m,\ \forall\, j=1,\cdots, l;
\label{matrixB}\\
{G}_{ij} &= (\nabla q_i, \nabla q_j), \qquad \forall\, i,j=1, \cdots, m;
\label{matrixG}\\
M_{ij} &= (q_j, q_j), \qquad \forall\, i,j=1, \cdots, m.
\label{matrixM}
\end{align}

The transpose of $B$ is written as $B^t$.
Denote the vector $\vec{f}$ as
$\vec{f}=((F,w_1)  \cdots  (F,w_l))^t$.
We express the following functions as the linear combinations of basis functions, i.e.,  $u^s = \sum_{i=1}^l u^s_i w_i$,
$\phi^s = \sum_{i=1}^m \phi^s_i q_i$,
$p^s = \sum_{i=1}^m p^s_i q_i$ and
denote the corresponding vector forms as

\begin{align}
\vec{u}^{s} = \left(
\begin{array}{c}
u^s_1 \\ \vdots \\ u^s_l
\end{array}
\right),\quad
\vec{\phi}^{s} = \left(
\begin{array}{c}
\phi^s_1 \\ \vdots \\ \phi^s_m
\end{array}
\right),\quad
\vec{p}^{\,s} = \left(
\begin{array}{c}
p^s_1 \\ \vdots \\ p^s_m
\end{array}
\right).
\end{align}

Thereafter, the transformation between the function form and the vector form of a quantity will be used in this way.
Thus, the numerical scheme,
(<ref>), (<ref>), and (<ref>),
is turned into the following matrix-vector form,

\begin{align}
A\vec u^s + k B^t \vec{p}^{\,s} &= \vec f,
\label{wei041}\\
G\vec{\phi}^s &=\frac{1}{k} B\vec{u^s},
\label{wei042}\\
M\vec{p}^{\,s+1} &= M\vec{p}^{\,s} + \alpha M\vec{\phi}^s + \rho B\vec{u}^s.
\label{wei043}
\end{align}

Without the convection,  the weak solution of (<ref>) and (<ref>)  is $(u,p)\in V_h\times Q_h$ satisfying

\begin{align}
1.5(u, {w_i}) + k\nu (\nabla u, \nabla{w_i})
%+ k(({w}\cdot\nabla) u, {w_i})
- k(\nabla\cdot {w_i}, p) = (F,{w_i}),
&\qquad\forall\, i=1,\cdots, l,
\label{ip011}\\
(\nabla\cdot u, q_i) = 0, &\qquad\forall\, i=1, \cdots, m.
\label{ip013}
\end{align}

The corresponding vector form is

\begin{align}
A\vec u + k B^t \vec p &= \vec f,
\label{wei044}\\
G\vec{\phi} &=\frac{1}{k} B\vec{u},
\label{wei045}\\
M\vec{p} &= M\vec{p} + \alpha M\vec{\phi} + \rho B\vec{u}.
\label{wei046}
\end{align}

The last two equations (<ref>) and (<ref>) can be reduced to $\left(G + \frac{\alpha}{k \rho} M \right)\vec{\phi} = 0$, whose solution is $\vec{\phi}=0$ since both $G$ and $M$ are symmetric and positive definite. Thus, (<ref>) and (<ref>) are equivalent to

\begin{align}
\label{Bu0}
\end{align}

But we keep the current form to compare with (<ref>) and (<ref>).

Denote the errors in the function form as

\begin{equation}
{\bar{u}}^{s} = {u}^s - {u}, \qquad
{\bar{\phi}}^{s} = {\phi}^s - {\phi}, \qquad
{\bar{p}}^{s} = {p}^s - {p},
\end{equation}

and the corresponding errors in the vector form as

\begin{equation}
\vec{\bar{u}}^{s} = \vec{u}^s - \vec{u}, \qquad
\vec{\bar{\phi}}^{s} = \vec{\phi}^s - \vec{\phi}, \qquad
\vec{\bar{p}}^{\,s} = \vec{p}^{\,s} - \vec{p}.
\end{equation}

The subtractions of these equations,
(<ref>)-(<ref>),  give rise to

\begin{align}
A\vec{\bar{u}}^s + k B^t \vec{\bar{p}}^{\,s} &= 0,
\label{kun001}\\
G \vec{\bar{\phi}}^s &= \frac{1}{k} B \vec{\bar{u}}^s,
\label{kun002}\\
\vec{\bar{p}}^{\,s+1} &= \vec{\bar{p}}^{\,s} + \alpha \vec{\bar{\phi}}^s + \rho M^{-1} B \vec{\bar{u}}^s.
\label{kun003}
\end{align}

From (<ref>) and (<ref>) we get

\begin{align}
\vec{\bar{u}}^s &= -k A^{-1}B^t \vec{\bar{p}}^{\,s}, \label{ubyp}\\
\vec{\bar{\phi}}^s &= - G^{-1}BA^{-1} B^t \vec{\bar{p}}^{\,s}.
\label{phibyp}
\end{align}

\begin{equation}
\end{equation}

which is the negative Schur complement of $A$ (e.g., see [5]),

\begin{equation}
\vec{\bar{p}}^{\,s+1} = ( I-(\alpha G^{-1} + \rho k M^{-1} )D )
\vec{\bar{p}}^{\,s}
\triangleq (I-K) \vec{\bar{p}}^{\,s},
\label{keyrelation01}
\end{equation}

\begin{equation}
K_{m\times m}=(\alpha G^{-1}  + \rho k M^{-1} )D.
\label{def_K}
\end{equation}

To study the convergence, we will analyze the eigenvalues of the matrix $K$.
First, we introduce the following lemmas.

Let $S$, $T$ be symmetric matrices. If $S$ is positive definite, then $S^{-1}T$ is diagonalizable.

Since $S$ is symmetric and positive definite (SPD), its inverse $S^{-1}$ is also SPD. Denote its square root as $S^{-1/2}$. Then $S^{1/2} S^{-1}TS^{-1/2}=S^{-1/2}TS^{-1/2}$,
where the latter is symmetric because $T$ and $S^{-1/2}$ are both symmetric. Therefore, $S^{-1}T$ is similar to a symmetric matrix and hence diagonalizable.

All the eigenvalues of $K$ are positive and $K$ has a set of linearly independent eigenvectors that span $\mathbb{R}^m$, which is called an eigenbasis of $K$.

First, note that without the convection term the matrix $A$ is SPD. The inf-sup condition guarantees $B$ has the full row rank and then the matrix $D$ is also SPD. In addition, the matrices $G$ and $M$ are also SPD.

Suppose $\lambda$ is an eigenvalue of $K$ and $\vec{\xi}\in \mathbb{R}^m$ is an associated eigenvector, i.e.,
$K\vec{\xi}=\lambda\vec{\xi}$.   Left-multiplying both sides by $\vec{\xi}^{\,t}D$ leads to

\begin{equation}
\alpha \vec{\xi}^{\,t} DG^{-1}D \vec{\xi}
+ \rho k \vec{\xi}^{\,t} DM^{-1}D \vec{\xi} = \lambda \vec{\xi}^{\,t} D \vec{\xi}.
\label{kun006}
\end{equation}

Since the quantities $\alpha$, $\rho$, $k$, $\vec{\xi}^{\,t} DG^{-1}D \vec{\xi}$, $\vec{\xi}^{\,t} DM^{-1}D \vec{\xi}$, $\vec{\xi}^{\,t} D \vec{\xi}$ are all positive, the value of $\lambda$ must be positive.

Furthermore, a matrix  is diagonalizable if and only if there is a set of eigenbasis of this matrix (e.g., [20] Theorem 1.3.7, p. 59).
In the matrix $K$, the left factor $(\alpha G^{-1} + \rho k M^{-1})$ is SPD and the right factor $D$ is symmetric.
According to Lemma <ref>,  $K$ has an eigenbasis.

Denote an eigenbasis of $K$ as
$\{ \vec{\xi}_1, \cdots, \vec{\xi}_m \}$ and the corresponding eigenvalues as $\lambda_j, j=1,\cdots, m$.
Let $\rho(I-K)$ be the spetral radius of $I-K$, that is,

\begin{equation}
\rho(I-K)=  \max_{j=1,\cdots, m} |1-\lambda_j|.
\label{kun000}
\end{equation}

Because all eigenvalues are positive according to Lemma <ref>,  it is easy to deduce from (<ref>) that

$\rho(I-K)<1$ if and only if $\lambda_{\max}<2$.

The next lemma gives an upper bound of the largest eigenvalue $K$ dependent on the parameters in the numerical scheme.

Let $\lambda_{\max}$ be the largest eigenvalue of $K$, then

\begin{equation}
\lambda_{\max}
\le  \max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right).
\label{cat104}
\end{equation}

Let $\lambda$ be an eigenvalue of $K$ and $\vec{\xi}$ be an associated eigenvector. According to (<ref>),

\begin{equation}
\lambda = \frac{\alpha \langle G^{-1}D\vec{\xi}, D\vec{\xi}\rangle
+ \rho k \langle M^{-1}D\vec{\xi}, D\vec{\xi}\rangle}{\langle D\vec{\xi}, \vec{\xi}\rangle},
\label{kun007}
\end{equation}

where the notation $\langle \vec{\xi}, \vec{q}\rangle\triangleq \vec{\xi}^{\,t} \vec{q}$ for any two vectors in the same Euclidean space $\mathbb{R}^m$ or $\mathbb{R}^l$.
Denote $\vec{v}=A^{-1}B^t \vec{\xi}$, i.e., $A\vec{v}=B^t\vec{\xi}$ or
Then (<ref>) becomes

\begin{equation}
\lambda
=\frac{\alpha \langle G^{-1}B\vec{v}, B\vec{v}\rangle
+ \rho k \langle M^{-1} B\vec{v}, B\vec{v}\rangle}{\langle \vec{v}, A\vec{v}\rangle}
=\frac{\alpha \langle G^{-1}B\vec{v}, B\vec{v}\rangle
+ \rho k \langle M^{-1} B\vec{v}, B\vec{v}\rangle}
{1.5\|v\|^2_{L^2(\Omega)} + \nu k\|\nabla v\|^2_{L^2(\Omega)}}
\triangleq R(\vec{v}).
\label{kun008}
\end{equation}

This equates the eigenvalues to the stationary values of the Rayleigh quotient $R(\vec{v})$ defined in (<ref>).

To estimate the two numerator terms in (<ref>), let $\vec{\phi}=G^{-1}B\vec{v}$,
then $G\vec{\phi}=B\vec{v}$. Using the connection between the vector and function forms:
$\vec{\phi}=(\phi_1 \cdots  \phi_m )^t$ and $\phi=\sum_{i=1}^m \phi_i q_i $,
$\vec{v}=(v_1 \cdots v_l )^t$ and $v=\sum_{i=1}^l v_i w_i$, we see that $G\vec{\phi}=B\vec{v}$  is equivalent to

\begin{equation}
\int_\Omega\nabla \phi \cdot \nabla \eta \mathrm{d}x = -\int_\Omega (\nabla\cdot v) \eta \mathrm{d}x, \quad
\forall \eta \in Q_h.
\label{kun009}
\end{equation}

\begin{equation}
\langle G^{-1}B\vec{v}, B\vec{v}\rangle = \langle \vec{\phi}, B\vec{v}\rangle
= -\int_\Omega \phi (\nabla\cdot v) \mathrm{d}x.
\end{equation}

Using $\eta=\phi$ in (<ref>), integration by parts, and Young's inequality, we get

\begin{equation}
\|\nabla \phi\|_{L^2(\Omega)}^2 =
-\int_\Omega \phi (\nabla\cdot v)\mathrm{d}x
=\int_\Omega v\cdot  \nabla \phi \mathrm{d}x
\le \frac{1}{2}\big( \|v\|^2_{L^2(\Omega)} + \|\nabla \phi\|^2_{L^2(\Omega)}
\big),
\end{equation}

which further gives $\|\nabla \phi\|_{L^2(\Omega)}^2\le \|v\|^2_{L^2(\Omega)}$.

\begin{equation}
\langle G^{-1}B\vec{v}, B\vec{v}\rangle  = \|\nabla \phi\|^2_{L^2(\Omega)} \le  \|v\|^2_{L^2(\Omega)}.
\label{kun009a}
\end{equation}

To study $\langle M^{-1} B\vec{v}, B\vec{v}\rangle $, we let $\vec{y}=M^{-1}B\vec{v}$, i.e., $M\vec{y}=B\vec{v}$.
Suppose the vector $\vec{y}=(y_1 \cdots y_m )^t$ and  the corresponding function  $y=\sum_{i=1}^m y_i q_i$. Then $M\vec{y}=B\vec{v}$ can be rewritten as

\begin{equation}
\int_\Omega y \eta \mathrm{d}x = -\int_\Omega (\nabla\cdot v) \eta \mathrm{d}x,
\quad \forall \eta \in Q_h.
\label{kun010}
\end{equation}

Letting $\eta=y$, then

\begin{equation}
\|y\|^2_{L^2(\Omega)} = -\int_\Omega (\nabla\cdot v)y \mathrm{d}x
\le \frac{1}{2}\big(
\|\nabla\cdot v\|^2_{L^2(\Omega)} + \|y\|^2_{L^2(\Omega)}
\big).
\end{equation}

This results in $\|y\|_{L^2(\Omega)} \le \|\nabla\cdot v\|_{L^2(\Omega)}$.
It is also known that for any $v\in [H_0^1(\Omega)]^3$, for example,  [29] p. 93,

\begin{equation}
\|\nabla\cdot v\|_{L^2(\Omega)} \le \|\nabla v\|_{L^2(\Omega)}.
\label{kun204}
\end{equation}

\begin{equation}
\langle M^{-1} B\vec{v}, B\vec{v}\rangle  = \langle \vec{y}, B\vec{v}\rangle
=  -\int_\Omega (\nabla\cdot v) y \mathrm{d}x
= \|y\|^2_{L^2(\Omega)} \le \|\nabla\cdot v\|_{L^2(\Omega)}^2
\le \|\nabla v\|_{L^2(\Omega)}^2.
\label{kun011}
\end{equation}

Collecting (<ref>), (<ref>), (<ref>), we obtain

\begin{equation}
\lambda_{\max}
\le \frac{\alpha \|v\|^2_{L^2(\Omega)}
+ \rho k \|\nabla v\|^2_{L^2(\Omega)} }
{1.5\|v\|^2_{L^2(\Omega)} + \nu k\|\nabla v\|^2_{L^2(\Omega)}}
\le \max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right).
\label{cat004}
\end{equation}

We introduce the maximum vector norm with respect to the eigenbasis  $\{\vec{\xi}_i\}_{i=1}^m$ of $K$, that is,

\begin{align}
\|\vec{q}\|_{K,\infty} = \max_{1\le j\le m} |q_{K,j}| \text{ when }
\vec{q} = \sum_{j=1}^m q_{K,j} \vec{\xi}_j.
\end{align}

The pressure solution error $\vec{\bar{p}}^{\,s}$ between the scheme (<ref>), (<ref>),  (<ref>) and the system (<ref>), (<ref>), (<ref>) satisfies

\begin{equation}
\| \vec{\bar{p}}^{\,s+1}  \|_{K,\infty}
\le  \rho(I-K)  \|\vec{\bar{p}}^{\,s} \|_{K,\infty}.
\label{kun033}
\end{equation}

Thus, if $\rho(I-K)<1$, then the iterative solution (${u}^s$, ${p}^{\,s}$) of  (<ref>), (<ref>), (<ref>) converges linearly to (${u}$, ${p}$) of (<ref>), (<ref>).
In particular, if $\max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right)<2$, then $\rho(I-K)<1$. Furthermore, the convergence of iteration is independent of the time step size $k$.

Express the error
$\vec{\bar{p}}^{\,s}= \sum_{j=1}^m \bar{p}^s_{K,j} \vec{\xi}_j$. Then
(<ref>) becomes

\begin{eqnarray}
\sum_{j=1}^m \bar{p}^{s+1}_{K,j} \vec{\xi}_j
= \sum_{j=1}^m (1-\lambda_j) \bar{p}^s_{K,j} \vec{\xi}_j.
\end{eqnarray}

Due to the linear independency of $\{\vec{\xi}_j\}_{j=1}^m$, we arrive at

\begin{equation}
\bar{p}^{s+1}_{K,j} = (1-\lambda_j) \bar{p}^s_{K,j},
\quad j=1,\cdots, m.
\end{equation}

Taking maximum values on both sides yields (<ref>).
By Lemma <ref> and Lemma <ref>, if $\max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right)<2$, then $\rho(I-K)<1$.

§.§.§ Iteration analysis with semi-implicit convection

When the semi-implicit convection form $\text{NL}(w,u)$ is used, the weak form of the saddle point problem (<ref>) and (<ref>) is composed by (<ref>) and

\begin{align}
1.5({u}, {w_i}) + k\nu (\nabla u, \nabla{w_i})
+ k( \text{NL}(w, u), w_i)
- k(\nabla\cdot {w_i}, p) = (F,{w_i}),
\quad \forall\, i=1,\cdots, l.
\label{cat005}
\end{align}

Accordingly, the weak form  of (<ref>) turns to

\begin{align}
1.5({u^s}, {w_i}) + k\nu (\nabla {u^s}, \nabla{w_i})
+ k( \text{NL}(w, u^s), w_i)
- k(\nabla\cdot {w_i}, p^s) = (F,{w_i}),
\quad \forall\, i=1,\cdots, l.
\label{cat001}
\end{align}

The entire semi-implicit iterative projection scheme includes (<ref>), (<ref>), (<ref>).
Define the matrices $N(w)$ and $A_{N}(w)$ as, respectively,

\begin{align}
N(w)_{ij} &= (\text{NL}(w,w_j), w_i),
\quad \forall\, i,j=1,\cdots, l,
\label{matrix_Ac}\\
A_{N}(w) &= A+kN(w).
\label{matrix_Atilde}
\end{align}

Then (<ref>) and (<ref>) can be written as

\begin{align}
A_{N}(w) \vec{u} + kB^t \vec{p} &= \vec{f},
\label{cat003} \\
A_{N}(w) \vec{u}^s + kB^t \vec{p}^{\,s} &= \vec{f}.
\label{cat002}
\end{align}

The proofs of iteration convergence with both SI and FI convection forms are by writing the system as a perturbation with the coefficient $k$ of the non-convection system, (<ref>). This is done for SI schemes in (<ref>), (<ref>), and for FI schemes in (<ref>). The perturbation terms with the coefficient $k$ depend on the convection matrix $N(w)$, which has an upper bound given in Lemma <ref>. Then, a sufficiently small $k$ keeps the system as a contraction.
To earn this upper bound, we adopt the spectral matrix norm $\vertiii{\cdot}_2$ for square matrices, which is induced by the Euclidean metric in $\mathbb{R}^l$. For more details of this norm, see, e.g., <cit.>. As for the spectral norm of $N(w)$, the following upper bound on $\|w\|_{L^2(\Omega)}$ when $w\in V_h$ is provided.

There exists a constant $C$ only dependent on $h$ and the basis of $V_h$ such that

\begin{align}
|(\text{NL}(u,v),w)| &\le  C \|u\|_{L^2(\Omega)}
\|v\|_{L^2(\Omega)} \|w\|_{L^2(\Omega)},
\quad \forall u,v,w\in V_h,
\label{kun300}\\
\vertiii{ N(w) }_2 &\le C  \|w\|_{L^2(\Omega)},
\quad \forall w\in V_h.
\label{kun301}
\end{align}

Using <cit.>,
$\vertiii{N(w)}_2 = \max_{\|\vec{\xi}\|_2=\|\vec{\eta}\|_2=1} {\vec{\xi}}^t N(w) \vec{\eta}$, where $\vec{\xi}, \vec{\eta}\in \mathbb{R}^l$.
Using the connection between the vector and function forms:
$\vec{\xi}=(\xi_1 \cdots  \xi_l )^t$,
$\xi=\sum_{i=1}^l \xi_i w_i $,
$\vec{\eta}=(\eta_1 \cdots  \eta_l )^t$,
$\eta=\sum_{i=1}^l \eta_i w_i$,
we obtain
$\vec{\xi}^t N(w) \vec{\eta} =
\sum_{i,j=1}^l \xi_i ((w\cdot \nabla)w_i, w_i) \eta_j
=((w\cdot\nabla)\xi, \eta)$.

\begin{equation}
\vertiii{N(w)}_2  =\max_{\|\vec{\xi}\|_2=\|\vec{\eta}\|_2=1}  ((w\cdot\nabla)\xi, \eta).
\label{kun207}
\end{equation}

Note that

\begin{align}
|((w\cdot\nabla)\xi, \eta)|  &=
\bigg| \int_\Omega ((w\cdot \nabla) \xi) \cdot \eta \mathrm{d}x \bigg|,  \label{kun200}\\
&\le  \|w\|_{L^4(\Omega)} \|\eta\|_{L^4(\Omega)} \|\nabla \xi\|_{L^2(\Omega)} \label{kun201}\\
&\le 2 \|w\|^{\frac14}_{L^2(\Omega)} \|\nabla w\|^{\frac34}_{L^2(\Omega)}
\|\eta\|^{\frac14}_{L^2(\Omega)}
\|\nabla \eta\|^{\frac34}_{L^2(\Omega)}
\|\nabla \xi\|_{L^2(\Omega)},
\label{kun202}
\end{align}

where (<ref>) is obtained by Hölder's inequality and
(<ref>) by Ladyzhenskaya's inequality in three-dimensional case,
$\|v\|_{L^4(\Omega)} \le \sqrt{2} \|v\|^{1/4}_{L^2(\Omega)}
\|\nabla v\|^{3/4}_{L^2(\Omega)}$, $\forall v\in V$ (see, e.g., <cit.>).
Next, by applying the inverse inequality $\|\nabla v\|_{L^2(\Omega)}
\le Ch^{-1} \|v\|_{L^2(\Omega)}$ for any $v\in V_h$ (e.g., <cit.>), where $C$ only depends on $\Omega$, the above estimates continue as

\begin{align}
|((w\cdot \nabla)\xi, \eta)| &\le 2 \|w\|^{\frac14}_{L^2(\Omega)}
(Ch^{-1})^{\frac34} \|w\|^{\frac34}_{L^2(\Omega)}
\|\eta\|^{\frac14}_{L^2(\Omega)}
\|\eta\|^{\frac34}_{L^2(\Omega)}
&\le Ch^{-\frac52} \|w\|_{L^2(\Omega)}
\|\xi\|_{L^2(\Omega)}
\|\eta\|_{L^2(\Omega)}\\
&= C \|w\|_{L^2(\Omega)}
\|\xi\|_{L^2(\Omega)}
\|\eta\|_{L^2(\Omega)},
\label{kun203}
\end{align}

where the last $C$ absorbs $h^{-2.5}$ in the previous step.
Note the inequality (<ref>)  holds for any $\xi,\eta \in V_h$, which leads to (<ref>).

But $\|\xi\|^2_{L^2(\Omega)} = (\sum_{i=1}^l \xi_i w_i,\sum_{j=1}^l \xi_j w_j) = \sum_{i,j=1}^l \xi_i \xi_j (w_i, w_j) = \vec{\xi}^t A_0 \vec{\xi}$, where $A_0$ is the mass matrix of $V_h$ with $(A_0)_{ij}=(w_i,w_j)$.  Clearly, $A_0$ is symmetric and positive definite. Since $\|\vec{\xi}\|_2=1$, $\vec{\xi}^t A_0 \vec{\xi}\le \rho(A_0)$, the spectral radius of $A_0$. Thus, $\|\xi\|_{L^2(\Omega)} \le \sqrt{\rho(A_0)}$ and similarly
$\|\eta\|_{L^2(\Omega)} \le \sqrt{\rho(A_0)}$. Therefore,

\begin{align}
|((w\cdot \nabla)\xi, \eta)| \le C \rho(A_0)
\|w\|_{L^2(\Omega)}
\le C \|w\|_{L^2(\Omega)},
\label{kun205}
\end{align}

where $C\rho(A_0)$ is denoted as a new $C$ in the final step,  which  only depends on $h$ and the basis of $V_h$. Finally,  (<ref>) is achieved from (<ref>) and (<ref>).

If the time step size $k$ is sufficiently small, then the symmetric part of $A_N(w)$, $H\triangleq \frac12 (A_N(w) + A_N^t(w))$, is positive definite.

This implies the system of (<ref>) and  (<ref>)
is a generalized saddle point problem according to  [5].

Let $\vec{v}=(v_1 \cdots v_l)^t\in\mathbb{R}^l$ and $v=\sum_{i=1}^l v_i w_i$, then

\begin{align}
\vec{v}^t H \vec{v} &= \vec{v}^t A\vec{v}
+  \frac{k}{2} \vec{v}^t (N(w) + N^t(w)) \vec{v}
\nonumber \\
&= \|v\|^2_{L^2(\Omega)} + k \|\nabla v\|^2_{L^2(\Omega)}
+ \frac{k}{2} ( (\text{NL}(w),v),v) + \text{NL}(v,w),w)) \nonumber \\
&\ge \|v\|^2_{L^2(\Omega)}  - kC \|w\|_{L^2(\Omega)} \|v\|^2_{L^2(\Omega)}
= (1-kC\|w\|_{L^2(\Omega)}) \|v\|^2_{L^2(\Omega)},
\nonumber
\end{align}

where $(\text{NL}(w),v),v), (\text{NL}(v),w),v)\ge - C  \|w\|_{L^2(\Omega)} \|v\|^2_{L^2(\Omega)}$ from Lemma <ref> is used in the last step. If $k<1/(C\|w\|_{L^2(\Omega)})$, then $H$ is positive definite.

Using the same process of obtaining (<ref>) in the case without convection, we obtain that the pressure solution error between the scheme (<ref>), (<ref>), (<ref>) and the system (<ref>),  (<ref>), (<ref>) satisfies

\begin{equation}
\vec{\bar{p}}^{\,s+1}  = (I-K_N(w)) \vec{\bar{p}}^{\,s},
\label{kun112}
\end{equation}

\begin{equation}
K_N(w)\triangleq (\alpha G^{-1} + \rho k M^{-1}) B A^{-1}_N(w) B^t.
\label{kun124}
\end{equation}

The equation (<ref>) can be regarded as a first order perturbation  in time step size $k$ of (<ref>).
The following lemma indicates that $I-K_N(w)$ is also a first order perturbation of $I-K$.

For each $w\in V_h$, there exists $k_w>0$, such that when  $0<k<\min(k_w,1)$, the matrix $A_N(w)$ is invertible,

\begin{equation}
\vertiii{A^{-1}_N(w)}_2 \le C(1+k \|w\|_{L^2(\Omega)}),
\label{kun500}
\end{equation}

\begin{equation}
I-K_N(w) = I-K - kA_{N2}(w),
\label{kun501}
\end{equation}

\begin{equation}
\vertiii{A_{N2}(w)}_2 \le  C \|w\|_{L^2(\Omega)}.
\label{kun502}
\end{equation}

Here, the constants $C$ are independent of $w$, but depend on $h$ and bases of $V_h$ and $Q_h$.

Note $A_N(w) = A+kN(w) = A (I+kA^{-1} N(w))$. When $w$ and the basis of $V_h$ are fixed, the matrices $A$ and $N(w)$ are also fixed. Note $A_N$ is invertible if and only if $(I+kA^{-1} N(w))$ is invertible.
According to <cit.>, $(I+kA^{-1} N(w))$ being invertible is equivalent to $k\vertiii{A^{-1} N(w)}_2 <1$.

\begin{equation}
k_w = \frac{0.5}{C \vertiii{A^{-1}}_2  \|w\|_{L^2(\Omega)}},
\label{def_kw}
\end{equation}

where the value of $C$ is the one in (<ref>).
If $0<k<k_w$,  then
$k \vertiii{A^{-1} N(w)}_2  \le
k \vertiii{A^{-1} }_2  \vertiii{N(w)}_2
\le k C \vertiii{A^{-1} }_2  \|w\|_{L^2(\Omega)} <0.5$.
Thus, the matrix $(I+kA^{-1} N(w))$ is invertible and  $(I+kA^{-1} N(w))^{-1}=I+kA_{N1}(w)$,  where

\begin{eqnarray}
A_{N1}(w) \triangleq  \sum_{j=1}^{\infty} (-1)^j k^{j-1} (A^{-1} N(w))^j
\label{def_A1}
\end{eqnarray}

\begin{equation}
\vertiii{A_{N1}(w)}_2\le \frac{\vertiii{A^{-1} N(w)}_2}{1-k\vertiii{A^{-1} N(w)}_2} \le
2 \vertiii{A^{-1}}_2 \vertiii{ N(w)}_2.
\label{kun111}
\end{equation}

Thus, when $0<k<k_w$, $A_N(w)$ is invertible and

\begin{align}
A^{-1}_N(w) = (I+k  A_{N1}(w)) A^{-1}
= A^{-1} + kA_{N1}(w) A^{-1}.
\label{formula_Atilde}
\end{align}

\begin{align}
\vertiii{A^{-1}_N(w)}_2
&= \vertiii{ (I+k  A_{N1}(w)) A^{-1}}_2
\le \vertiii{A^{-1}}_2 (1 + k \vertiii{A_{N1}(w)}_2) \label{kun499}\\
&\le \vertiii{A^{-1}}_2 (1 + 2k \vertiii{A^{-1}}_2
\vertiii{ N(w)}_2) \label{kun498}\\
&\le \vertiii{A^{-1}}_2 (1 + 2C k \vertiii{A^{-1}}_2
\|w\|_{L^2(\Omega)})
\label{kun497}\\
&\le  C' (1+k \|w\|_{L^2(\Omega)}),
\label{kun496}
\end{align}

where (<ref>) is used in (<ref>) to get (<ref>), and (<ref>) is applied in (<ref>) to get (<ref>). The $C'$ in (<ref>) is $C'=\max(\vertiii{A^{-1}}_2,2C\vertiii{A^{-1}}^2_2)$ from (<ref>). Changing $C'$ to a new $C$ results in (<ref>).

Then (<ref>) and (<ref>)  lead to

\begin{align}
I-(\alpha G^{-1} + \rho k M^{-1}) B A^{-1} B^t
- k (\alpha G^{-1} + \rho k M^{-1}) BA_{N1}(w) A^{-1} B^t  \\
&\triangleq I-K - k A_{N2}(w),
\label{kun122}
\end{align}

where $K$ is defined in (<ref>) and

\begin{eqnarray}
A_{N2}(w) \triangleq (\alpha G^{-1} + \rho k M^{-1}) BA_{N1}(w) A^{-1} B^t.
\label{def_A2}
\end{eqnarray}

\begin{align}
\vertiii{A_{N2}(w)}_2 &\le
( \alpha\vertiii{ G^{-1}}_2 +   \rho k \vertiii{M^{-1}}_2 )
\vertiii{B}_2
\vertiii{B^t}_2 \vertiii{A^{-1}}_2 \vertiii{A_{N1}(w)}_2 \\
&\le ( \alpha\vertiii{ G^{-1}}_2 +   \rho  \vertiii{M^{-1}}_2 )
\vertiii{B}_2
\vertiii{B^t}_2
\cdot 2 \vertiii{A^{-1}}^2_2  \vertiii{N(w)}_2
\label{kun120}\\
&=  C_1 \vertiii{N(w)}_2 \le C  \|w\|_{L^2(\Omega)},
\label{kun121}
\end{align}

where we have used $0<k<1$ and (<ref>) to get  (<ref>) and Lemma <ref> to obtain (<ref>).
The constant $C_1$ in (<ref>) denotes the product of all terms in front of $\vertiii{N(w)}_2$ in (<ref>), which only depends on the bases of $V_h$ and $Q_h$.
The identity (<ref>) and estimate (<ref>)  conclude this lemma.

The next result presents the recursive relation of the pressure error and the iteration convergence conditions.

If $k$ is sufficiently small, then the pressure error between (<ref>), (<ref>), (<ref>) and the system (<ref>), (<ref>), (<ref>), where $w$ is  explicit,  satisfies

\begin{equation}
\| \vec{\bar{p}}^{\,s+1}  \|_{K,\infty}
\le \big( \rho(I-K)
+ k C  \|w\|_{L^2(\Omega)}  \big)
\cdot \|\vec{\bar{p}}^{\,s} \|_{K,\infty},
\label{kun503}
\end{equation}

where $C$ is a constant that only depends on $h$  and the bases of $V_h$ and $Q_h$. If $\rho(I-K)<1$, which can be guaranteed by enforcing $\max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right)<2$,  and $k$ is such that $\rho(I-K) + k C  \|w\|_{L^2(\Omega)}<1$, then the iteration converges. In this case, the solution $(u^s,p^s)$ of (<ref>), (<ref>), (<ref>) converges linearly to $(u,p)$ of (<ref>), (<ref>).

According to (<ref>) and Lemma <ref>,

\begin{equation}
\vec{\bar{p}}^{\,s+1} = (I-K) \vec{\bar{p}}^{\,s}
- kA_{N2}(w) \vec{\bar{p}}^{\,s}.
\label{kun402}
\end{equation}

We express $A_{N2}(w) \vec{\bar{p}}^{\,s}$ as a linear combination of the eigenbasis $\{\vec{\xi}_j \}_{j=1}^m$, i.e.,

\begin{equation}
A_{N2}(w) \vec{\bar{p}}^{\,s} = \sum_{j=1}^m \eta_{K,j} \vec{\xi}_j
\label{kun401}
\end{equation}

for some real numbers $\eta_{K,j}\in \mathbb{R}$, $j=1, \cdots, m$. Then there exists a constant $C$ merely dependent on the  bases of $V_h$ and $Q_h$ such that

\begin{align}
\|\vec{\eta}\|_{K,\infty} &=
\max_{1\le j\le m} |\eta_{K,j}| \le
C \| A_{N2}(w) \vec{\bar{p}}^{\,s} \|_2
\le C \vertiii{A_{N2}(w)}_2 \|\vec{\bar{p}}^{\,s}\|_2\\
&\le C  \|w\|_{L^2(\Omega)}
\|\vec{\bar{p}}^{\,s} \|_{K,\infty},
\label{kun403}
\end{align}

where (<ref>) is used in the last step and the equivalence of norms in $\mathbb{R}^m$ is repeatedly used.

Using the expressions $\vec{\bar{p}}^{\,s+1}
=\sum_{j=1}^m \bar{p}^{s+1}_{K,j} \vec{\xi}_j$
and $\vec{\bar{p}}^{\,s}
=\sum_{j=1}^m \bar{p}^{s}_{K,j} \vec{\xi}_j$,
(<ref>) and (<ref>) yield
$\sum_{j=1}^m \bar{p}^{s+1}_{K,j} \vec{\xi}_j =
\sum_{j=1}^m ( (1-\lambda_j) \bar{p}^s_{K,j} - k\eta_{K,j})
\vec{\xi}_j$.
The linear independency of the eigenbasis $\{\vec{\xi} \}_{j=1}^m$ implies

\begin{equation}
\bar{p}^{s+1}_{K,j}=
(1-\lambda_j) \bar{p}^s_{K,j} - k\eta_{K,j}, \quad j=1, \cdots, m.
\end{equation}

This suggests

\begin{eqnarray}
\| \vec{\bar{p}}^{\,s+1}  \|_{K,\infty}
&\le&  \max_{1\le j \le m} |1-\lambda_j|  \cdot
\| \vec{\bar{p}}^{\,s}  \|_{K,\infty}
+ k \| \vec{\bar{\eta}}  \|_{K,\infty} \nonumber \\
&\le& \rho(I-K)
\| \vec{\bar{p}}^{\,s}  \|_{K,\infty}
+ k C  \|w\|_{L^2(\Omega)}
\|\vec{\bar{p}}^{\,s} \|_{K,\infty},
\end{eqnarray}

where (<ref>) is used in the last step.

When $\max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right)<2$, Lemma <ref> and Lemma <ref> imply $\rho(I-K)<1$.
In the recursive relation (<ref>),  if $k$ is sufficiently small, the iteration ratio satisfies $0<  \rho(I-K)
+ k C  \|w\|_{L^2(\Omega)}  <1$. In this case, the iteration converges linearly.

§.§.§ Iteration analysis with fully implicit convection

The fully implicit and  nonlinear system of (<ref>) and (<ref>) in the weak form is,  after deleting the time step index $(n+1)$,

\begin{align}
1.5 (u, v) + k\nu (\nabla u, \nabla v)
+ k (\text{NL}(u, u), v) - k(p, \nabla\cdot v) &= (F,v),
\quad \forall v\in V_h,
\label{aa001}\\
(\nabla\cdot u, q) &=0, \quad \forall q\in Q_h.
\label{aa002}
\end{align}

The corresponding FI-Iterative Projection Method is, for $s=0, 1, \cdots$,

\begin{align}
1.5 (u^s, v) + k\nu (\nabla u^s, \nabla v)
+ k (\text{NL}(u^{s-1}, u^s), v) - k(p^s, \nabla\cdot v) &= (F,v),
\quad \forall v\in V_h,
\label{aa003}
\end{align}

along with (<ref>) and (<ref>).
Here, $(u^{-1}, p^0)$ is the initial guess of the solution.
The matrix form of (<ref>) is

\begin{align}
A_N(u)  \vec{u} + k B^t \vec{p} = \vec{f},
\label{kunaa011}
\end{align}

and that of (<ref>) is

\begin{align}
A_N(u^{s-1}) \vec{u}^{s} + k B^t {\vec p}^{\,s} = \vec{f},
\label{kunaa010}
\end{align}

where the matrix $A_N$ is defined in (<ref>) and (<ref>).

The convergence of this scheme is given below.

If $\rho(I-K)<1$ and the time step size $k$ is sufficient small, then there exists a constant $0<\gamma<1$ such that for any $s\ge 0$,
the pressure and velocity errors between (<ref>),
(<ref>), (<ref>) and the system (<ref>), (<ref>), (<ref>) satisfy

\begin{align}
\|\vec{\bar{p}}^{\,s} \|_{K,\infty} \le
+ \|\vec{\bar{u}}^{-1}\|_2) \gamma^s,
\quad
\|\vec{\bar{u}}^{\,s} \|_{2} \le  \sqrt{k}
+ \|\vec{\bar{u}}^{-1}\|_2) \gamma^{s+1}.
\label{kun604}
\end{align}

Thus, the solution $(u^s, p^s)\in V_h\times Q_h$ of (<ref>), (<ref>) and (<ref>) converges linearly to the solution $(u,p)\in V_h\times Q_h$ of (<ref>) and (<ref>), as $s\to \infty$.
In particular, if $\max\left( \frac{\alpha}{1.5}, \frac{\rho}{\nu}\right)<2$, then $\rho(I-K)<1$.

Step 1. Derive the relations (<ref>) and (<ref>).
Subtracting (<ref>) from (<ref>) gives, for any $v\in V_h$,

\begin{eqnarray}
1.5(\bar{u}^s,v)  + k\nu (\nabla\bar{u}^s, \nabla v) +
k(\text{NL}(u^{s-1}, \bar{u}^s), v)
+ kb(\text{NL}(\bar{u}^{s-1}, u), v)
-k(\bar{p}^s, \nabla\cdot v)=0.
\end{eqnarray}

This results in, by using the notations of vectors
$\vec{u}^s$, $\vec{\phi}^s$, $\vec{p}^{\,s}$ and errors
$\vec{\bar{u}}^s$, $\vec{\bar{\phi}}^s$, $\vec{\bar{p}}^{\,s}$ in Section <ref>,

\begin{eqnarray}
A\vec{\bar{u}}^s + kN(u^{s-1}) \vec{\bar{u}}^s
+ kB^t\vec{\bar{p}}^{\,s} = -kA_d(u) \vec{\bar{u}}^{s-1},
\end{eqnarray}

where  $A+kN(u^{s-1}) = A_N(u^{s-1})$ and  the new matrix $A_d(u)$ is defined as

\begin{eqnarray}
(A_d(u))_{ij} &\triangleq &
(\text{NL}(w_j, u), w_i), \quad i,j=1,\cdots, l.
\end{eqnarray}

According to Lemma <ref>, $\vertiii{N(u^{s-1})}_2  \le C \|u^{s-1}\|_{L^2(\Omega)}$.
According to Lemma <ref>, for each $u^{s-1}$, there exists $k_{u^{s-1}}>0$ that depends on the norm $\|u^{s-1}\|_{L^2(\Omega)}$ (see (<ref>)) such that
when $0<k<\min(k_{u^{s-1}},1)$,
$A_N(u^{s-1})$ is invertible.
In Step 2 of this proof, a non-zero lower bound of $k_{u^{s}}$ with respective to $s=-1, 0, \cdots$ is found in (<ref>), $\frac{0.5}{CC_2 \vertiii{A^{-1}}_2}$, such that $A_N(u^s)$ is invertible for all $s=-1, 0, \cdots$ when (<ref>) is satisfied.
According to Lemma <ref> and (<ref>),
$\vertiii{A_d(u)}_2 \le C \|u\|_{L^2(\Omega)}$.
In this problem, $u$ is the solution of (<ref>) and (<ref>), therefore it is a fixed quantity and thus in the following we denote it as $A_d$ and

\begin{equation}
\vertiii{A_d}_2 \le C.
\label{kun601}
\end{equation}

\begin{eqnarray}
\vec{\bar{u}}^s= - k A^{-1}_N(u^{s-1})
(B^t \vec{\bar{p}}^{\,s} + A_d \vec{\bar{u}}^{s-1}),
\end{eqnarray}

\begin{align}
\|\vec{\bar{u}}^s\|_2 &\le k
\vertiii{ A^{-1}_N(u^{s-1})}_2
(\vertiii{B^t}_2 \|\vec{\bar{p}}^{\,s}\|_2
+ \vertiii{A_d}_2 \|\vec{\bar{u}}^{s-1}\|_2 )\nonumber\\
&\le k C (1+k\|u^{s-1}\|_{L^2(\Omega)} )
( \|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+  \|\vec{\bar{u}}^{s-1}\|_2),
\label{kun602}
\end{align}

where (<ref>) and (<ref>) are used in the last step.

Similar to the derivation of (<ref>) and (<ref>) in Section <ref>, we obtain

\begin{align}
\vec{\bar{p}}^{\,s+1}
&= \vec{\bar{p}}^{\,s} +
\bigg(\frac{\alpha}{k} G^{-1} +\rho M^{-1}\bigg) B\vec{\bar{u}}^{\,s}  \\
&= \vec{\bar{p}}^{\,s} -
\big(\alpha G^{-1} +\rho  k M^{-1}\big)
B  A^{-1}_N(u^{s-1})
(B^t \vec{\bar{p}}^{\,s} + A_d \vec{\bar{u}}^{s-1})
\label{vi006}
\\
&= \vec{\bar{p}}^{\,s} -
K_N(u^{s-1})  \vec{\bar{p}}^{\,s}
- (\alpha G^{-1} +\rho  k M^{-1} )
A^{-1}_N(u^{s-1}) A_d \vec{\bar{u}}^{s-1}
\label{vi007}
\\
&= ( I - K - k A_{N2}(u^{s-1})) \vec{\bar{p}}^{\,s}
- B_1(u^{s-1})  \vec{\bar{u}}^{s-1},
\label{vi008}
\end{align}

where the definition $K_N(u^{s-1})$ in (<ref>) is used in (<ref>) to obtain (<ref>),
the formula (<ref>) is employed to go from (<ref>) to (<ref>),
and $B_1(u^{s-1})$ is defined as

\begin{equation}
B_1(u^{s-1}) \triangleq (\alpha G^{-1} +\rho  k M^{-1} )
A^{-1}_N(u^{s-1}) A_d.
\end{equation}

\begin{align}
\vertiii{B_1(u^{s-1})}_2 &\le (\alpha \vertiii{G^{-1}}_2
+ \rho k \vertiii{M^{-1}}_2)
\vertiii{(A^{-1}_N(u^{s-1}) }_2
\vertiii{A_d}_2 \nonumber\\
&\le C (1+k\|u^{s-1}\|_{L^2(\Omega)}),
\end{align}

where (<ref>) and (<ref>) are used in the last step.

Expressing $A_{N2}(u^{s-1}) \vec{\bar{p}}^{\,s}$
and $B_1(u^{s-1})  \vec{\bar{u}}^{s-1}$ in (<ref>)
with the eigenbasis $\{\vec{\xi}_j\}_{j=1}^m$ of the matrix $K$:

\begin{align}
A_{N2}(u^{s-1}) \vec{\bar{p}}^{\,s} =
\sum_{j=1}^m \eta_{K,j} \vec{\xi}_{j} = \vec{\eta},
\quad
B_1(u^{s-1})  \vec{\bar{u}}^{s-1}
= \sum_{j=1}^m \theta_{K,j} \vec{\xi}_{j} = \vec{\theta}.
\label{vi009}
\end{align}

Then it is easy to tell that

\begin{align}
\|\vec{\eta}\|_{K,\infty} &\le
C \|A_{N2}(u^{s-1}) \vec{\bar{p}}^{\,s} \|_2
\le C \vertiii{A_{N2}(u^{s-1})}_2
\|\vec{\bar{p}}^{\,s} \|_{K,\infty} \nonumber \\
&\le C \|u^{s-1} \|_{L^2(\Omega)}
\|\vec{\bar{p}}^{\,s} \|_{K,\infty},
\end{align}

where (<ref>) is used in the last step.

\begin{align}
\|\vec{\theta}\|_{K,\infty} &\le
C \|B_1(u^{s-1})  \vec{\bar{u}}^{s-1} \|_2
\le C \vertiii{B_1(u^{s-1})}_2
\|\vec{\bar{u}}^{s-1} \|_2  \nonumber \\
&\le C  (1+k\|u^{s-1}\|_{L^2(\Omega)})
\|\vec{\bar{u}}^{s-1} \|_2.
\end{align}

From (<ref>) and (<ref>), we obtain, for $j=1, \cdots, m$,

\begin{align}
\vec{\bar{p}}^{\,s+1}_{K,j}
= (1-\lambda_j - k\eta_{K,j} ) \vec{\bar{p}}^{\,s}_{K,j}
- \theta_{K,j}.
\end{align}

This implies

\begin{align}
\|\vec{\bar{p}}^{\,s+1}\|_{K,\infty}
&\le (\rho(I-K) + k \|\vec{\eta}\|_{K,\infty} )
\|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+  \|\vec{\theta}\|_{K,\infty} \nonumber\\
(\rho(I-K) + k C \|u^{s-1} \|_{L^2(\Omega)} )
\|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+  C  (1+k\|u^{s-1}\|_{L^2(\Omega)})
\|\vec{\bar{u}}^{s-1} \|_2 \\
&= a_{s-1} \|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+ b_{s-1} \|\vec{\bar{u}}^{s-1} \|_2,
\label{kun603}
\end{align}

where for $s=-1, 0, \cdots$,

\begin{align}
a_{s} &= \rho(I-K) + k C \|u^{s} \|_{L^2(\Omega)}
\le \rho(I-K) + kC ( \|\vec{\bar{u}}^s \|_2 + \|u\|_{L^2(\Omega)}), \label{def_as}
\\
b_{s} &= C  (1+k\|u^{s}\|_{L^2(\Omega)})
\le C( 1+ k ( \|\vec{\bar{u}}^s \|_2 + \|u\|_{L^2(\Omega)})).
\label{def_bs}
\end{align}

Therefore, using these two quantities, (<ref>) and (<ref>) become, for $s=0, 1, \cdots$,

\begin{align}
\|\vec{\bar{p}}^{\,s+1}\|_{K,\infty}
&\le  a_{s-1} \|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+ b_{s-1} \|\vec{\bar{u}}^{s-1} \|_2,
\label{kun605}\\
\|\vec{\bar{u}}^s\|_2 &\le  k b_{s-1}
( \|\vec{\bar{p}}^{\,s}\|_{K,\infty}
+  \|\vec{\bar{u}}^{s-1}\|_2).
\label{kun606}
\end{align}

Step 2. prove (<ref>) by induction.
Fix $C$ as the maximum of one and all $C$'s that appear above (this enforces $C\ge 1$), which only depends on $h$ and bases of $V_h$ and $Q_h$.

\begin{align}
+ \|\vec{\bar{u}}^{-1}\|_2, \label{def_C1} \\
C_2&=\max(\|u^{-1}\|_{L^2(\Omega)} , CC_1+\|u\|_{L^2(\Omega} ), \label{def_C2} \\
P_1&=\rho(I-K) + C\sqrt{k} + kC(1+\sqrt{k})(\sqrt{k} C_1+\|u\|_{L^2(\Omega)}), \label{def_P1} \\
P_2&=\sqrt{k}(1+\sqrt{k}) C (1+k\sqrt{k})(\sqrt{k} C_1+\|u\|_{L^2(\Omega)}), \label{def_P2}  \\
P_3&=\sqrt{k} C (1+ k\| u^{-1}\|_{L^2(\Omega)}).
\label{def_P3}
\end{align}

If $\rho(I-K)<1$, then there exists $k_0>0$ such that when
$0<k<k_0$,   $\max(P_1, P_2, P_3)<1$ because $C, C_1, u, u^{-1}$ are fixed.
We fix a value $k$ such that

\begin{align}
0< k<\min \left(1, k_0,  \frac{0.5}{CC_2 \vertiii{A^{-1}}_2} \right).
\label{def_k}
\end{align}

The value of $\gamma$ is chosen satisfying

\begin{equation}
\max(P_1, P_2, P_3)<\gamma<1.
\label{def_gamma}
\end{equation}

In the induction, we not only need to prove (<ref>), but also  need to show $A_N(u^{s})$ is invertible for $s=-1, 0, \cdots$.

In the initial step of induction, we show $A_N(u^{-1})$ is invertible, $\|\vec{\bar{p}}^{\,0}\|_{K,\infty}\le C_1$,
$\|\vec{\bar{u}}^{\,0}\|_{2}\le \sqrt{k} C_1 \gamma$,
and $A_N(u^0)$ is invertible.
First, (<ref>) implies
$k<\frac{0.5}{CC_2 \vertiii{A^{-1}}_2}$, which in turn suggests
$k< \frac{0.5}{C \vertiii{A^{-1}}_2 \|u^{-1}\|_{L^2(\Omega)} }$ due to (<ref>). According to Lemma <ref> and (<ref>), $A_N(u^{-1})$ is invertible.

When $s=0$, the inequality for $\|\vec{\bar{p}}^{\,0}\|_{K,\infty}$ in (<ref>) is trivially true due to the definition of $C_1$ in (<ref>). From (<ref>),

\begin{align}
\|\vec{\bar{u}}^{\,0}\|_{2} &\le
k b_{-1}
( \|\vec{\bar{p}}^{\,0}\|_{K,\infty}
+  \|\vec{\bar{u}}^{-1}\|_2) \nonumber\\
&= k C(1+k \|u^{-1}\|_{L^2(\Omega)})
( \|\vec{\bar{p}}^{\,0}\|_{K,\infty}
+  \|\vec{\bar{u}}^{-1}\|_2)
=  \sqrt{k} C_1 P_3 \le \sqrt{k} C_1 \gamma.
\end{align}

Next, we show $A_N(u^0)$ is invertible.

\begin{align}
\|u^0\|_{L^2(\Omega)} &\le
\|\bar{u}^0\|_{L^2(\Omega)} + \|u\|_{L^2(\Omega)}
\le
C \|\vec{\bar{u}}^0\|_2 + \|u\|_{L^2(\Omega)}   \nonumber\\
&\le C \sqrt{k} C_1 \gamma + \|u\|_{L^2(\Omega)}
\le CC_1 + \|u\|_{L^2(\Omega)},
\label{formula_u0}
\end{align}

where $0<k<1$ and $0<\gamma<1$ are used in the last step.
According to (<ref>), (<ref>), (<ref>),

\begin{align}
k < \frac{0.5}{C C_2 \vertiii{A^{-1}}_2}
\le \frac{0.5 }{C \vertiii{A^{-1}}_2 (CC_1 + \|u\|_{L^2(\Omega)})}
\le \frac{0.5 }{C \vertiii{A^{-1}}_2  \|u^0\|_{L^2(\Omega)}}
= k_{u^0}.
\end{align}

Due to Lemma <ref> and (<ref>), $A_N(u^0)$ is invertible.

Second, assume the inequalities in (<ref>) hold and $A_N(u^{s})$ is invertible for $s\le s_0$ where $s_0\ge 0$.
Then, (<ref>) implies

\begin{align}
\|\vec{\bar{p}}^{\,s_0+1}\|_{K,\infty}
&\le  a_{s_0-1} \|\vec{\bar{p}}^{\,s_0}\|_{K,\infty}
+ b_{s_0-1} \|\vec{\bar{u}}^{s_0-1} \|_2 \nonumber\\
[ \rho(I-K) + kC (\|\vec{\bar{u}}^{s_0}\|_2
+  \|u\|_{L^2(\Omega)}) ]
\|\vec{\bar{p}}^{\,s_0}\|_{K,\infty} \label{kun541}\\
+ C [ 1+k  (\|\vec{\bar{u}}^{s_0}\|_2
+  \|u\|_{L^2(\Omega)}) ]
\|\vec{\bar{u}}^{s_0-1} \|_2  \label{kun542} \\
&\le [ \rho(I-K) + kC(\sqrt{k}C_1 \gamma^{s_0} +
\|u\|_{L^2(\Omega)} ) ] C_1 \gamma^{s_0} \label{kun543}\\
C[ 1 + k( \sqrt{k} C_1 \gamma^{s_0}+ \|u\|_{L^2(\Omega)} )]
\sqrt{k}C_1\gamma^{s_0} \label{kun544}\\
&=C_1\gamma^{s_0} \bigg[
\rho(I-k) + \sqrt{k} C + kC (1+\sqrt{k}) (\sqrt{k}C_1\gamma^{s_0}  + \|u\|_{L^2(\Omega)} )
\bigg]  \label{kun545} \\
&< C_1\gamma^{s_0} \bigg[
\rho(I-k) + \sqrt{k} C + kC (1+\sqrt{k}) (\sqrt{k}C_1 + \|u\|_{L^2(\Omega)} )
\bigg] \label{kun546}\\
&= C_1 \gamma^{s_0} P_1 \le  C_1 \gamma^{s_0+1},  \label{kun547}
\end{align}

where the inequality (<ref>) and (<ref>) are used in attaining (<ref>), $0<\gamma<1$ is applied from(<ref>) to (<ref>), the definition $P_1$ in (<ref>) is employed to go from (<ref>) to (<ref>), and the relation (<ref>) is utilized in (<ref>).

Furthermore, (<ref>) implies

\begin{align}
\|\vec{\bar{u}}^{s_0+1}\|_2 &\le kb_{s_0}
+ \|\vec{\bar{u}}^{s_0}\|_2) \nonumber \\
&\le kC [1+k(\|\vec{\bar{u}}^{s_0}\|_2
+ \|u\|_{L^2(\Omega)}  )] \cdot
+ \|\vec{\bar{u}}^{s_0}\|_2) \label{kun531}\\
&\le kC[ 1 + k(\sqrt{k}C_1 \gamma^{s_0+1} + \|u\|_{L^2(\Omega)} )] \cdot (C_1 \gamma^{s_0+1} + \sqrt{k} C_1 \gamma^{s_0+1} ) \label{kun532} \\
&= \sqrt{k} C_1 \gamma^{s_0+1}
\bigg[
\sqrt{k} (1+\sqrt{k})  C
( 1+ k(\sqrt{k}C_1 \gamma^{s_0+1} + \|u\|_{L^2(\Omega)} )
\bigg] \label{kun533} \\
&< \sqrt{k} C_1 \gamma^{s_0+1}
\bigg[
\sqrt{k} (1+\sqrt{k})  C
( 1+ k(\sqrt{k}C_1 + \|u\|_{L^2(\Omega)} )
\bigg] \label{kun534} \\
&= \sqrt{k} C_1 \gamma^{s_0+1} P_2
\le
\sqrt{k} C_1 \gamma^{s_0+2}, \label{kun535}
\end{align}

where the inequality (<ref>) is used in achieving (<ref>), $0<\gamma<1$ is used in obtaining (<ref>), the definition $P_2$ in (<ref>) is applied from (<ref>) to (<ref>), and the relation (<ref>) is used in (<ref>).

Finally, following the same process as we prove that $A_N(u^0)$ is invertible, we can prove $A_N(u^{s_0+1})$ is invertible.  The induction proof is completed.

§ STABILITY AND ERROR ANALYSIS OF CONVERGENT ITERATIVE PROJECTION METHOD

From here, we denote the solution at time step $n$ with  iteration number $s$ in the finite element space $V_h\times Q_h$  as $(u_h^{n,s}, p_h^{n,s})$.
Thus, the finite element formulation of the iterative scheme (<ref>)–(<ref>) is, at the time  $t_n$,  seeking ${u}^{n+1,s}_h\in V_h$, $p^{n+1,s}_h \in Q_h$, $\phi^{n+1,s}_h\in Q_h$, $s=0, 1, \cdots$, such that for any ${ v_h}\in V_h$, $q_h\in Q_h$,

\begin{align}
1.5 ( u_h^{n+1,s}, v_h)
&+k\nu (\nabla u_h^{n+1,s}, \nabla  v_h)
+ k (\text{NL}(\tilde{w}_h^{n+1,s},u_h^{n+1,s}), v_h)
- k (\nabla\cdot v_h, p^{n+1,s}_h)
\nonumber \\
%% equation 1, 2nd half:
(2u^n_h -0.5u^{n-1}_h, v_h)
+ k ( f^{n+1}, v_h),
\label{eqn_u}
\\
%% equation 2:
(\nabla \phi^{n+1,s}_h,  \nabla q_h )
(\nabla\cdot { u}^{n+1,s}_h, q_h),
\label{eqn_phi_FEM}
\\
%% equation 3:
(p^{n+1,s+1}_h, q_h)
( p^{n+1,s}_h + \alpha \phi_h^{n+1,s} -
\rho (\nabla\cdot u^{n+1,s}_h),  q_h ),
\label{p_update}
\end{align}

where $p_h^{n+1,0}=2p_h^n - p_h^{n-1}$.
In the semi-implicit convection, $\tilde{w}^{n+1,s}_h = 2u^n_h -u^{n-1}_h$.
In the fully implicit convection case, $\tilde{w}^{n+1,s}_h = u^{n+1,s-1}_h$ and when $s=0$, $u^{n+1,-1}_h = 2u^n_h -u^{n-1}_h$.
When the conditions in Theorem <ref> or Theorem <ref> are satisfied,   the iterative solution ($u_h^{n+1,s}$, $p_h^{n+1,s}$)  converges linearly to ($u_h^{n+1}$, $p_h^{n+1}$) of the following system

\begin{align}
1.5 ( u_h^{n+1}, v_h)
k\nu (\nabla u_h^{n+1}, \nabla  v_h)
+ k(\text{NL}(\tilde{w}_h^{n+1},u_h^{n+1}), v_h)
- k (\nabla\cdot v_h, p^{n+1}_h)
\nonumber \\
%% equation 1, 2nd half:
(2u^n_h -0.5u^{n-1}_h, v_h)
+ k ( f^{n+1}, v_h),
\quad \forall v_h\in V_h,
\label{eqn_u_skew_conv}\\
(\nabla \cdot u^{n+1}_h, q_h) &= 0,
\quad \forall q_h\in Q_h,
\label{eqn_p_skew_conv}
\end{align}

where $\tilde{w}^{n+1}_h=2u^n_h-u^{n-1}_h$ for the semi-imcplicit convection and $\tilde{w}^{n+1}_h=u^{n+1}_h$ for the fully implicit convection.

According to [17],  the stability result of the traditional rotational  projection method, of only one iterative projection in each time step,  is not proven for the scheme (<ref>), (<ref>), and (<ref>).
This is because the intermediate terms $u_h^{n+1,s}$ and $\phi_h^{n+1,s}$  make the analysis very complicated in the multi-step time discretization. Nevertheless, as for the limit scheme (<ref>) and  (<ref>),  the stability and error analyses are  established below.
In what follows, we use $\|\cdot\|_{L_x^2}$ to denote $\|\cdot\|_{L^2(\Omega)}$ for simplicity.

§.§ Stability analysis for the convergent iterative projection method

Suppose the convection term adopts SI-SKEW or FI-SKEW with both non-div-free and div-free FEMs.
If $0<k<1$, then the velocity solution of the scheme
(<ref>) and  (<ref>)  satisfies for a given final time $T>0$, at any time step $2\le n\le \lfloor \frac{T}{k} \rfloor$,

\begin{equation}
\|u_h^{n}\|_{L^2_x}^2 + \|2u_h^{n}-u_h^{n-1}\|_{L^2_x}^2
\le  \exp\Big(\frac{T}{1-k}\Big) \cdot  \Big( \|u_h^{1}\|_{L^2_x}^2 + \|2u_h^{1}-u_h^0\|_{L^2_x}^2
+ 4k \sum_{j=2}^{n} \|f^{j}\|_{L^2_x}^2 \Big).
\label{cat8910}
\end{equation}

It is easy to derive that for any $w,v\in V$,

\begin{align}
(\text{NL}(w,v),v) = 0.
\label{kun036}
\end{align}

Thus, if  SI-SKEW or FI-SKEW form is used in (<ref>), the convection term
$(\text{NL}(\tilde{w}_h^{n+1}, u_h^{n+1}), u^{n+1}_h)$
Letting $v_h=u_h^{n+1}$ in (<ref>) and
$q_h=p^{n+1}_h$ in (<ref>) eliminates the convection term  and $k(\nabla\cdot u^{n+1}_h,p^{n+1}_h)$,
and  we obtain

\begin{equation}
(u^{n+1}_h, 1.5u^{n+1}_h-2u^n_h+0.5 u^{n-1}_h) + k\nu \|\nabla u^{n+1}_h\|^2_{L^2_x}
= k (f^{n+1}, u^{n+1}_h).
\label{nan001}
\end{equation}

Applying the following algebraic identity for real numbers $a^{n+1}, a^n, a^{n-1}$ (e.g., (4.7) in [19]),

\begin{align}
a^{n+1} \cdot (1.5a^{n+1}- 2a^n+ 0.5 a^{n-1})
=\frac{1}{4} \big( &|a^{n+1}|^2 + |2a^{n+1}-a^n|^2 + |a^{n+1}-2a^n+a^{n-1}|^2  \nonumber \\
& - |a^n|^2 - |2a^n-a^{n-1}|^2 \big)
\label{Shen_identity}
\end{align}

on the first term in (<ref>), we obtain

\begin{align}
&\,\|u^{n+1}_h\|^2_{L^2_x} + \|2u^{n+1}_h - u^n_h\|^2_{L^2_x} + \|u^{n+1}_h -2u^n_h + u^{n-1}_h\|^2_{L^2_x}
+ 4k\nu \|\nabla u^{n+1}_h\|^2_{L^2_x} \nonumber \\
= &\,\|u^{n}_h\|^2_{L^2_x} + \|2u^{n}_h-u^{n-1}_h\|^2_{L^2_x} + 4k (f^{n+1}, u^{n+1}_h).
\end{align}

Dropping the third and fourth terms on the left and replacing the last term on the right by the upper bound in Young's inequality $|4k (f^{n+1}, u^{n+1}_h)| \le 4k \|f^{n+1}\|^2_{L^2_x} + k \|u^{n+1}\|^2_{L^2_x}$, we obtain

\begin{align}
(1-k) \|u^{n+1}_h\|^2_{L^2_x}
+ \|2u^{n+1}_h-u^n_h\|^2_{L^2_x}
\le \|u^{n}_h\|^2_{L^2_x}
+ \|2u^{n}_h - u^{n-1}_h\|^2_{L^2_x}
+ 4k \|f^{n+1}\|^2_{L^2_x}.
\label{cat8921}
\end{align}

Denoting $A^{n}=\|u^{n}_h\|^2_{L^2_x}
+ \|2u^{n}_h-u^{n-1}_h\|^2_{L^2_x}$, then (<ref>) yields

\begin{equation}
A^{n+1} \le \frac{1}{1-k} A^n+ \frac{4k}{1-k} \|f^{n+1}\|^2_{L^2_x}.
\end{equation}

Repeating this inequality backward leads to

\begin{align}
A^{n+1} \le \frac{1}{(1-k)^n} A^1 +  4k \sum_{j=2}^{n+1} \frac{1}{(1-k)^{n+2-j}} \|f^{j}\|^2_{L^2_x}.
\end{align}

Using a basic inequality $(1+x)^j\le \exp(jx)$ for $x>-1$ and $j\in\mathbb{N}$, we get

\begin{equation}
\frac{1}{(1-k)^j} = \Big( 1 + \frac{k}{1-k} \Big)^j \le \exp\Big(\frac{jk}{1-k}\Big)
\le \exp\Big(\frac{nk}{1-k}\Big),
\quad j=1,\cdots, n.
\end{equation}

\begin{align}
A^{n+1} \le  \exp\Big(\frac{nk}{1-k}\Big) \cdot \Big( A^1 +  4k \sum_{j=2}^{n+1} \|f^j\|^2_{L^2_x}
\Big).
\end{align}

At any time step $1<n\le \lfloor \frac{T}{k} \rfloor$,  one has $nk\le T$ and thus the estimate (<ref>) follows.

§.§ Error analysis for the convergent iterative projection method

In the error analysis, we adopt the symbols and approach from [23, 14]. The error analysis in  [23, 14]  does not have  time discretizations and focuses on the fully implicit convection case.
As in [14], we introduce the Stokes projection $\pi_s: V\to V^{div}_h$ where

\begin{equation}
V^{div}_h = \{v_h\in V_h: (\nabla\cdot v_h, q_h)=0, \ \forall\, q_h\in Q_h  \}
\end{equation}

\begin{equation}
(\nabla \pi_s (u), \nabla v_h) = (\nabla u, \nabla v_h), \quad
\forall\, v_h\in V^{div}_h.
\end{equation}

In addition, let the $L^2$-projection $\pi_Q: Q\to Q_h$ be

\begin{equation}
(\pi_Q(p), q) = (p,q), \quad \forall\, q\in Q_h.
\end{equation}

To simplify notations, we denote $H^m_x=H^m(\Omega)$, $m\in\mathbb{N}$, $L^\infty_x=L^\infty(\Omega)$, $L^2_t=L^2(0,T)$, $L^\infty_t=L^\infty(0,T)$, $L^\infty_t H^m_x=L^\infty(0,T; H^m(\Omega))$, $L^2_x L^\infty_t=L^2(\Omega; L^\infty(0,T))$, and so on.

Suppose the degrees of local finite element polynomials in $V_h$ and $Q_h$ are $r$ and $r-1$, respectively, where $r\ge 2$. There exist the following estimates  according to [14],

\begin{align}
\|u - \pi_s(u)\|_{H^m_x} &\le C h^{r+1-m} \|u\|_{H^{r+1}_x},
\quad m=0,1,
\label{John_001}\\
\|\pi_s(u)\|_{L^\infty_x} &\le C \|u\|_{H^2_x},\label{John_002}\\
\|\nabla \pi_s(u)\|_{L^\infty_x} &\le
C \|\nabla u\|_{L^\infty_x}, \label{John_003}\\
\|\nabla \pi_s(u)\|_{L^p_x} &\le
C \|\nabla u\|_{L^p_x},
\quad \quad\ \ \forall\, p\in[2,\infty), \label{John_004}\\
\|p-\pi_Q(p)\|_{H^m_x} &\le C h^{r-m} \|p\|_{H^r_x}, \quad m=0,1.
\label{John_005}
\end{align}

Furthermore, the followings hold

\begin{equation}
\|u\|_{L^\infty_x} \le ( \|u\|_{H^1_x}
\|u\|_{H^2_x} )^{1/2} \le \|u\|_{H^2_x},
\label{agmon}
\end{equation}

where the first one is Agmon's inequality.

Denote $s^n_h=\pi_s(u^n)$ and the error as

\begin{equation}
e^n=u^n - u^n_h = (u^n-s^n_h) - (u^n_h-s^n_h) \triangleq
\eta^n - \phi^n_h.
\label{err_defs}
\end{equation}

The following lemma is used in the proof of Theorem <ref>.

When $u_t\in L^\infty_xH^{r+1}_x$, then
$\|\eta^{n+1}-\eta^n\|_{L^2_x} \le Ck h^{r+1} \|u_t\|_{L^\infty_tH^{r+1}_x}$.

Denote $u(t)= u(t,x)$ as  a map from $[0,T]$ to $V$. Note

\begin{align}
\eta^{n+1} - \eta^n=& (u^{n+1}-\pi_s(u^{n+1})) - (u^n-\pi_s(u^n))
=(I-\pi_s)(u(t^{n+1}) - u(t^n))
\nonumber \\
=& \int_{t^n}^{t^{n+1}} \partial_t (I-\pi_s) u(t) dt
=\int_{t^n}^{t^{n+1}} (I-\pi_s) u_t(t) dt.
\nonumber
\end{align}

\begin{align}
\|\eta^{n+1}-\eta^n\|^2_{L^2_x} = &
\int_\Omega \bigg| \int_{t^n}^{t^{n+1}} (I-\pi_s) u_t(t) dt \bigg|^2 dx
\le \int_\Omega k  \int_{t^n}^{t^{n+1}} |(I-\pi_s) u_t(t)|^2 dt dx \label{cat1102} \\
=&\, k  \int_{t^n}^{t^{n+1}}  \int_\Omega  |(I-\pi_s) u_t(t)|^2 dx dt
\le\, k \int_{t^n}^{t^{n+1}}  \|(I-\pi_s) u_t(t)\|^2_{L^2_x}dt \label{cat1103}\\
\le&\, k^2 C h^{2r+2} \| u_t(t)\|^2_{L^\infty_tH^{r+1}_x}
\end{align}

where the Cauchy-Schwatz inequality is used in the last step of (<ref>), Fubini's theorem in the first step of (<ref>) to switch the integration order, and the estimate (<ref>) is applied in the last step.

The following definition is taken from  [14].

A scheme is  robust if there is no negative power of viscosity in the $L^2$ error bound of velocity, and not robust otherwise.

Assume the semi-implicit convection form, $\text{NL}(2u^n_h-u^{n-1}_h, u^{n+1}_h)$, is used in (<ref>).
Suppose $u, u_t\in L^\infty_t(H^{r+1}_x)$ with $r\ge 2$, $u_{tt}, u_{ttt}\in L^2_x(L^\infty_t)$ and $p\in L_t^\infty H^r_x$.  Let the time step size $k<1$ and the spatial mesh size $h<1$.  Then for any $2\le n \le \lfloor \frac{T}{k} \rfloor$, there exists a constant $C$ independent of the solution $(u,p)$ and the discretization parameters $k$ and $h$, such that the velocity solution of the scheme (<ref>) and  (<ref>)  satisfies,

(i) in SI-SKEW with non-div-free FEMs,

\begin{align}
&\|u^n - u^n_h\|^2_{L^2_x} +
\frac{k\nu}{2} \|\nabla (u^n-u^n_h)\|^2_{L^2_x}
\le  Ch^{2r} (h^2+ \nu k) \|u\|^2_{L^\infty_t H^{r+1}_x}
\nonumber\\
&\,+ \exp{\Big( \frac{(1+CM_1) T}{1-k}\Big)} \cdot
\Big(\|\phi^1_h\|^2_{L^2_x}
+ \|2\phi^1_h - \phi^0_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \|\nabla \phi^1\|^2_{L^2_x}
+ \frac{CM^2_2}{1+CM_1}
\Big).
\label{error_anal}
\end{align}

In this case, the scheme is not robust and has spatial accuracy order $r$.

(ii) in SI-SKEW with div-free FEMs,
the error estimate (<ref>) with $M_1$ replaced with $M^{div0}_1$ and $M_2$ replaced with $M^{div0}_2$. In this case, the scheme is robust with spatial accuracy order $r$.

The constants $M_1$, $M^{div0}_1$, $M_2$, $M^{div0}_2$ are defined as

\begin{align}
M_1 =&\, M^{div0}_1+  \frac{1}{\nu} \|u\|^2_{L^\infty_t H^{r+1}_x},\quad
M^{div0}_1 =\, \|u\|^2_{L^\infty_t H^{r+1}_x},
\quad
M_2 =\,
M^{div0}_2 + \frac{h^r}{\sqrt{\nu}} \|p\|_{L^\infty_t H^r_x},
\label{def_M1M2}\\
M^{div0}_2 =&\,
h^{r+1} \|u_t\|_{L^\infty_t H^{r+1}_x}
+ h^r \|u\|^2_{L^\infty_t H^{r+1}_x}
+  k^2 \big( \|u_{ttt}\|_{L^2_x L^\infty_t}
+ \|u\|_{L^\infty_t H^3_x}
\|u_{tt}\|_{L^2_x L^\infty_t}
\big).
\label{def_M2_div0}
\end{align}

In Case (i) of Theorem <ref>, , if all the norms of $u$ and $p$ are bounded and
$\phi^0_h=\phi^1_h=0$, then
when $\nu\to 0+$,

\begin{align}
\|u^u-u^n_h\|_{L_x^2} =&\,
O(h^{r+1}+h^r\sqrt{k\nu} )
+ \exp\Big(O\Big(\frac{1}{\nu} \Big)\Big)
\cdot O\big( h^{r+1}\sqrt{\nu} +
h^r + k^2\sqrt{\nu} ) \big),\\
\|\nabla(u^u-u^n_h)\|_{L_x^2} =&\,
O\left(\frac{h^{r+1}} {\sqrt{k\nu}} + h^r\right)
+ \exp\Big(O\Big(\frac{1}{\nu} \Big)\Big)
\cdot O\bigg( \frac{h^{r}} {\sqrt{k\nu}} + \frac{h^r}{\sqrt{k}}+ k^{1.5}  \bigg).
\end{align}

The spatial accuracy order of both the $L^2$ and $H^1$ errors of velocity is $r$. On the other hand,
when $h=0$, the  accuracy order in time step size for the velocity in $H^1$ norm is 1.5, the same as that for the rotational projection method [19] where the error is achieved without spatial discretizations.

The negative viscosity power in $M_1$ comes from $(\text{NL}(2\phi^n_h - \phi^{n-1}_h, s_h^{n+1}), \phi_h^{n+1})$, which produces the viscosity in (<ref>) from Young's inequality in order to suppress $\nabla\phi^{n+1}_h$. The negative viscosity power in $M_2$ is created from the pressure term
$(\nabla\cdot\phi_h^{n+1}, p^{n+1}-\pi_Q(p^{n+1}))$ by suppressing $\nabla\cdot\phi^{n+1}_h$ also using Young's inequality.

In Case (ii) of Theorem <ref>, if all the norms of $u$ and $p$ are bounded and
$\phi^0_h=\phi^1_h=0$, then
when $\nu\to 0+$,

\begin{align}
\|u^u-u^n_h\|_{L_x^2} =
O(h^{r} + k^2), \quad
\|\nabla(u^u-u^n_h)\|_{L_x^2} =
\frac{1}{\sqrt{k\nu}} O(h^{r} + k^2).
\end{align}

Since the $L^2$ error of velocity does not contain the negative power of $\nu$, this method is robust.

In the proof, we use $C$ to denote a generic constant which is independent of $u, p, k, h$, but may depend on the domain and the terminal time $T$. The value of the constant may vary line by line according to the context.

Step 1. Derivation of the error equation (<ref>).
We rewrite (<ref>) and (<ref>) as

\begin{align}
\Big( \frac{1.5  u_h^{n+1} - 2u^n_h + 0.5 u^{n-1}_h}{k}, v_h\Big)
&+ \nu (\nabla u_h^{n+1}, \nabla  v_h)
+ (\text{NL}(\tilde{w}^{n+1}_h, u^{n+1}_h),v_h)
- (\nabla\cdot v_h, p^{n+1}_h)
\nonumber \\
&= ( f^{n+1}, v_h), \qquad \forall\, v_h\in V_h,
\label{eqn_u_skew_02}\\
(\nabla \cdot u^{n+1}_h, q_h) &= 0,
\qquad \forall\, q_h\in Q_h,
\label{eqn_p_skew_02}
\end{align}

where $\tilde{w}^{n+1}_h=2u^n_h-u^{n-1}_h$ for SI convections and $u^{n+1}_h$ for FI convections.

Consider the Navier-Stokes equations in the following weak form

\begin{align}
\Big( \frac{1.5  u^{n+1} - 2u^n + 0.5 u^{n-1}}{k}, v\Big)
&+ \nu (\nabla u^{n+1}, \nabla  v)
+ (\text{NL}(u^{n+1}, u^{n+1}),v)
- (\nabla\cdot v, p^{n+1})
\nonumber \\
&= ( f^{n+1}, v) + (R^{n+1}(k),v),
\qquad \forall\, v\in V,
\label{eqn_NSu}\\
(\nabla \cdot u^{n+1}, q) &= 0,
\qquad \forall\, q\in Q,
\label{eqn_NSp}
\end{align}

\begin{equation}
(R^{n+1}(k), v) = \Big(
\frac{1.5 u^{n+1}-2u^n+0.5 u^{n-1}}{k} -u^{n+1}_t,v\Big)
\label{def_Rk}
\end{equation}

is the truncation term from the time discretization.

Subtracting (<ref>) from (<ref>) with $v=v_h$ in (<ref>), we obtain

\begin{align}
&\, \Big( \frac{1.5  e^{n+1} - 2e^n + 0.5 e^{n-1}}{k}, v_h\Big)
+ \nu (\nabla e^{n+1}, \nabla  v_h)
- (\nabla\cdot v_h, p^{n+1}- p^{n+1}_h)
\nonumber \\
=&\, B(v_h)
+ ( R^{n+1}(k), v_h),
\qquad\forall\, v_h\in V_h,
\label{eqn_NSu_02}
\end{align}

\begin{align}
B(v_h)=(\text{NL}(\tilde{w}^{n+1}_h, u^{n+1}_h),v_h)
- (\text{NL}(u^{n+1}, u^{n+1}),v_h).
\label{def_B}
\end{align}

Using the notations in (<ref>), the above equation becomes

\begin{align}
&\,\Big( \frac{1.5  \phi^{n+1}_h - 2\phi_h^n + 0.5 \phi_h^{n-1}}{k}, v_h\Big)
+ \nu (\nabla \phi_h^{n+1}, \nabla  v_h)
+ (\nabla\cdot v_h, p^{n+1}- p^{n+1}_h)
+ B(v_h)
\nonumber \\
\Big( \frac{1.5  \eta^{n+1} - 2\eta^n + 0.5 \eta^{n-1}}{k}, v_h\Big)
+ \nu (\nabla \eta^{n+1},\nabla v_h)
- ( R^{n+1}(k), v_h),
\qquad\forall\, v_h\in V_h.
\label{eqn_NSu_03}
\end{align}

Letting $v_h=\phi^{n+1}_h$ in (<ref>),
we get

\begin{align}
&\,\frac{1}{k} (1.5\phi^{n+1}_h-2\phi^n_h+0.5\phi^{n-1}_h,\phi^{n+1}_h)
+ \nu \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
+ B(\phi^{n+1}_h) \nonumber\\
&\quad + (\nabla\cdot \phi^{n+1}_h, p^{n+1}-\pi_Q (p^{n+1}) )
+ (\nabla\cdot \phi^{n+1}_h, \pi_Q(p^{n+1})-p^{n+1}_h)
\nonumber \\
=&\, \frac{1}{k} (1.5\eta^{n+1}-2\eta^n +0.5\eta^{n-1},\phi^{n+1}_h)
+ \nu (\nabla\eta^{n+1},\nabla\phi^{n+1}_h)
- (R^{n+1}(k),\phi^{n+1}_h).
\label{eqn_NSu_04}
\end{align}

Note that $(\nabla\cdot \phi^{n+1}_h, \pi_Q (p^{n+1})-p^{n+1}_h )=0$ because $\phi^{n+1}_h\in V^{div}_h$ and $\pi_Q (p^{n+1})-p^{n+1}_h\in Q_h$.

In addition, $ (\nabla\eta^{n+1},\nabla\phi^{n+1}_h)=0$
because $\eta^{n+1}=u^{n+1} - s^{n+1}_h$
and $(\nabla u^{n+1},\nabla\phi_h^{n+1}) = (\nabla s^{n+1}_h, \nabla\phi_h^{n+1})$ by Stokes projection with $\phi_h^{n+1}\in V^{div}_h$. Thus, (<ref>) becomes

\begin{align}
&\,\frac{1}{k} (1.5\phi^{n+1}_h-2\phi^n_h+0.5\phi^{n-1}_h,\phi^{n+1}_h)
+ \nu \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\nonumber\\
\frac{1}{k} (1.5\eta^{n+1}-2\eta^n +0.5\eta^{n-1},\phi^{n+1}_h)
-(\nabla\cdot \phi^{n+1}_h, p^{n+1}-\pi_Q (p^{n+1}) ) \notag\\
&\, -B(\phi^{n+1}_h) - (R^{n+1}(k),\phi^{n+1}_h).
\label{eqn_NSu_05}
\end{align}

Applying the identity (<ref>) on the first term  on the left of (<ref>) leads to

\begin{align}
&\,\frac{1}{k} \big(
\|\phi^{n+1}_h\|^2_{L^2_x}
+ \|2\phi^{n+1}_h-\phi^n_h\|^2_{L^2_x}
+ \|\phi^{n+1}_h - 2\phi^n_h + \phi^{n-1}_h\|^2_{L^2_x}
\big)
+ \nu \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\nonumber\\
\frac{1}{k}
\big(
\|\phi^{n}_h\|^2_{L^2_x}
+ \|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}
\big)
+ \frac{1}{k} (1.5\eta^{n+1}-2\eta^n +0.5\eta^{n-1},\phi^{n+1}_h)
\nonumber \\
- (\nabla\cdot \phi^{n+1}_h, p^{n+1}-\pi_Q (p^{n+1}) )
- (R^{n+1}(k),\phi^{n+1}_h)
-  B(\phi^{n+1}_h).
\label{eqn_NSu_06}
\end{align}

Step 2. Estimation of  common terms for SI and FI convections in  (<ref>).
This refers to the second, third, and fourth terms on the right of (<ref>).

As for the second term on the right of (<ref>), we have

\begin{align}
&\,\frac{1}{k} \big|(1.5\eta^{n+1}-2\eta^n +0.5\eta^{n-1},\phi^{n+1}_h) \big|
= \frac{1}{k}
|( 1.5(\eta^{n+1}-\eta^n)- 0.5(\eta^n-\eta^{n-1}),\phi^{n+1}) \big|
\nonumber \\
\le&\, \frac{1}{k} ( 1.5 \|\eta^{n+1}-\eta^n\|_{L^2_x}
+ 0.5 \|\eta^{n}-\eta^{n-1}\|_{L^2_x} )
\|\phi^{n+1}_h\|_{L^2_x}\label{cat2011} \\
\le &\, Ch^{r+1} \|u_t\|_{L^\infty_tH^{r+1}_x}
\|\phi^{n+1}_h\|_{L^2_x}
\le
\frac{1}{10}  \|\phi^{n+1}_h\|^2_{L^2_x}
+ Ch^{2r+2} \|u_t\|^2_{L^\infty_t H^{r+1}_x},
\label{mp001}
\end{align}

where Lemma <ref> is used in  (<ref>) and Young's inequality is used to earn (<ref>).

As for the third term on the right of (<ref>),
we can get

\begin{align}
\big|(\nabla\cdot \phi^{n+1}_h,  p^{n+1}-\pi_Q (p^{n+1}) )
\big| &\le \|\nabla\cdot \phi^{n+1}_h \|_{L^2_x} \| p^{n+1} - \pi_Q(p^{n+1}) \|_{L^2_x}
\\
&\le \|\nabla \phi^{n+1}_h \|_{L^2_x} Ch^r
\|p \|_{L^\infty_t H^r_x}
\label{cat1015} \\
&\le \frac{\nu}{4} \|\nabla\phi^{n+1}_h\|^2_{L^2_x} +
\frac{Ch^{2r}\|p\|^2_{L^\infty_t H^{r}_x}}{\nu},
\label{mp002}
\end{align}

where the inequality (<ref>) and the estimate (<ref>) with $m=0$ are applied to get (<ref>) and Young's inequality is used in the last step. One way to avoid negative powers of $\nu$ is throwing derivatives to the pressure, thus we can get

\begin{align}
&\big|(\nabla\cdot \phi^{n+1}_h,  p^{n+1}-\pi_Q (p^{n+1}) )
\big|
= \big|(\phi^{n+1}_h,  \nabla(p^{n+1}-\pi_Q (p^{n+1})) )
\big|  \nonumber \\
\le& \|\phi^{n+1}_h \|_{L^2_x} \| \nabla(p^{n+1} - \pi_Q(p^{n+1})) \|_{L^2_x}
\le \|\phi^{n+1}_h \|_{L^2_x} Ch^{r-1}
\|p \|_{L^\infty_t H^r_x}
\label{cat1015b} \\
\le& \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x}
+ Ch^{2r-2}\|p\|^2_{L^\infty_t H^{r}_x},
\label{mp002b}
\end{align}

where the estimate (<ref>) is used in (<ref>).

Note if a divergence-free FE pair is used, then $\nabla\cdot \phi^{n+1}_h=0$ in $L^2$ sense and
$(\nabla\cdot \phi^{n+1}_h,  p^{n+1}-\pi_Q (p^{n+1}) ) =0$.

As for the fourth term on the right of (<ref>),
because  $\frac{1.5u^{n+1}-2u^n+0.5 u^{n-1}}{k}$is an approximation of $u^{n+1}_t$  with  second order truncation error, there exists a positive constant $C_1$ such that for almost every $x\in\Omega$,

\begin{align}
\Big( \frac{1.5u^{n+1}-2u^n+0.5 u^{n-1}}{k} - u^{n+1}_t \Big) (x)
&= C_1 k^2 u_{ttt}(\xi_1(x),x),
\end{align}

where $\xi_1(x)\in (t^{n-1}, t^{n+1})$.

\begin{align}
&\, \big| (R^{n+1}(k),\phi^{n+1}_h) \big|
\le\, \big| (C_1  k^2 u_{ttt}(\xi_1(x),x), \phi^{n+1}_h(x))\big|   \nonumber\\
\le&\, Ck^2 \|u_{ttt}(\xi_1(x),x)\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x}
\le\, Ck^2 \|u_{ttt}\|_{L^2_x L^\infty_t}
\|\phi^{n+1}_h\|_{L^2_x}.
\label{cat2016a}
\end{align}

To reach (<ref>), we have used
$\|\nabla u^{n+1}\|_{L^\infty_x}  \le \|u\|_{L^\infty_t H^3_x}$
due to (<ref>), and
the following process

\begin{align}
\|u_{ttt}(\xi_1(x),x)\|_{L^2_x}
\int_\Omega |u_{ttt}(\xi_1(x),x)|^2 \mathrm{d}x
\Big)^{1/2}
\nonumber \\
\Big(
\int_\Omega \esssup_{t\in [0,T]} |u_{ttt}(t,x)|^2 \mathrm{d}x
\Big)^{1/2}
= \|u_{ttt}\|_{L^2_xL^\infty_t}.
\label{LxtoLt}
\end{align}

Using Young's inequality on (<ref>) yields

\begin{align}
\big| (R^{n+1}(k),\phi^{n+1}_h) \big|
\le
\frac{1}{10}\|\phi^{n+1}_h\|^2_{L^2_x}
+ Ck^4  \|u_{ttt}\|^2_{L^2_x L^\infty_t}.
\label{mp003}
\end{align}

Step 3. Estimation of SI convection term
$B(\phi^{n+1}_h)$ in  (<ref>).
In the SI convection case,  $B(\phi^{n+1}_h) =
(\text{NL}(2u^n_h-u^{n-1}_h, u^{n+1}_h), \phi^{n+1}_h)
-(\text{NL}(u^{n+1}, u^{n+1}), \phi^{n+1}_h)$.
We replace $u^j_h$ by $u^j-\eta^j + \phi^j_h$ for $j=n+1, n, n-1$ and obtain

\begin{align}
B(\phi^{n+1}_h) =&\, (\text{NL}((2u^n-u^{n-1}) - (2\eta^n-\eta^{n-1}) + (2\phi^n_h-\phi^{n-1}_h),
u^{n+1} - \eta^{n+1} + \phi^{n+1}_h), \phi^{n+1}_h)
\nonumber\\
- ( \text{NL}(2u^u - u^{n-1}, u^{n+1}), \phi^{n+1}_h)
+(\text{NL}(2u^u-u^{n-1}-u^{n+1}, u^{n+1}),\phi^{n+1}_h)
\\
=&\, -(\text{NL}(2u^n-u^{n-1},\eta^{n+1}), \phi^{n+1}_h)
-(\text{NL}(2\eta^n-\eta^{n-1},s^{n+1}_h), \phi^{n+1}_h) \nonumber \\
&\,+(\text{NL}(2s^n_h-s^{n-1}_h,\phi^{n+1}_h), \phi^{n+1}_h)
+(\text{NL}(2\phi^n_h-\phi^{n-1}_h,s^{n+1}_h),\phi^{n+1}_h) \nonumber \\
+(\text{NL}(2u^n - u^{n-1} - u^{n+1}, u^{n+1}),\phi^{n+1}_h).
\label{cat0112}
\end{align}

Based on (<ref>), the estimate of $B(\phi^{n+1}_h)$ is split into the following 6 parts.

(3A) Estimate of
$(\text{NL}(2u^n-u^{n-1},\eta^{n+1}), \phi^{n+1}_h)$.
Because $u^n$ and $u^{n-1}$ are exact solutions,
$\nabla\cdot(2u^n-u^{n-1})=0$ almost everywhere in $\Omega$.
Thus $\text{NL}(2u^n-u^{n-1},\eta^{n+1})
=((2u^n-u^{n-1})\cdot \nabla) \eta^{n+1}$.

\begin{align}
& | (\text{NL}(2u^n-u^{n-1},\eta^{n+1}), \phi^{n+1}_h) | \nonumber \\
\le&   \|2u^n-u^{n-1}\|_{L^\infty_x}
\|\nabla \eta^{n+1}\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x} \\
\le & 3 \|u\|_{L^\infty_tH^2_x}
Ch^r  \|u\|_{L^\infty_tH^{r+1}_x}
\|\phi^{n+1}_h\|_{L^2_x}
\\
\le &
\frac{1}{10}   \|\phi^{n+1}_h\|^2_{L^2_x}
+ Ch^{2r}  \|u\|^4_{L_t^\infty H^{r+1}_x},
\label{mp003a}
\end{align}

where (<ref>) is used to obtain
$\|2u^n-u^{n-1}\|_{L^\infty_x} \le 3\|u\|_{L^\infty_t H^2_x}$,
the estimate (<ref>) to get  $ \|\nabla \eta^{n+1}\|_{L^2_x}\le C h^r \|u\|_{L^\infty_t H^{r+1}_x}$,  and $r\ge 2$ and  Young's inequality are used in the last step.

Estimate of $(\text{NL}(2\eta^n-\eta^{n-1},s^{n+1}_h), \phi^{n+1}_h)$.
In the SI-SKEW convection form,

\begin{align}
&\, \big| (\text{NL}(2\eta^n-\eta^{n-1}, s_h^{n+1}), \phi^{n+1}_h)
\big| \nonumber \\
=&\, \big| (((2\eta^n-\eta^{n-1})\cdot\nabla)s_h^{n+1},\phi^{n+1}_h)
+ \frac{1}{2} ((\nabla\cdot (2\eta^n-\eta^{n-1})) s_h^{n+1},\phi^{n+1}_h) \big| \label{cat1104} \\
\le &\, \|\nabla s_h^{n+1}\|_{L^\infty_x}
\|2\eta^n-\eta^{n-1}\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x} \nonumber\\
&\qquad\qquad\qquad\qquad\qquad + \frac{1}{2} \|s^{n+1}_h\|_{L^\infty_x}
\|\nabla\cdot(2\eta^n-\eta^{n-1})\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x}  \\
\le&\, \big(C h^{r+1} \|u\|_{L^\infty_t H^{3}_x}
\|u\|_{L^\infty_t H^{r+1}_x}
+ Ch^{r}  \|u\|_{L^\infty_t H^{2}_x}
\|u\|_{L^\infty_t H^{r+1}_x}
\big)
\|\phi^{n+1}_h\|_{L^2_x}
\label{cat0311}\\
\le&\, Ch^r \|u\|^2_{L^\infty_t H^{r+1}_x}
\|\phi^{n+1}_h\|_{L^2_x}
\label{cat1011}\\
\le&\, \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x}
+  Ch^{2r}  \|u\|^{4}_{L^\infty_t H^{r+1}_x}.
\label{mp003b}
\end{align}

To achieve (<ref>), we have used $ \|\nabla s_h^{n+1}\|_{L^\infty_x}
\le C\|\nabla u\|_{L^\infty_t L^\infty_x}
\le C\|u\|_{L^\infty_t H^3_x}$
due to (<ref>) and (<ref>),
$ \|2\eta^n-\eta^{n-1}\|_{L^2_x}
\le Ch^{r+1} \|u\|_{L^\infty_t H^{r+1}_x}$
and $\|\nabla\cdot(2\eta^n-\eta^{n-1})\|_{L^2_x}\le Ch^{r} \|u\|_{L^\infty_t H^{r+1}_x}$
due to (<ref>),
and $\|s^{n+1}_h\|_{L^\infty_x}\le
C \|u\|_{L^\infty_t H^2_x}$
due to (<ref>).
To get (<ref>), we have used
$\|u\|_{L^\infty_t H^2_x}
\le \|u\|_{L^\infty_t H^3_x}
\le \|u\|_{L^\infty_t H^{r+1}_x}$ when $r\ge2$ and $h<1$.
In the step (<ref>), Young's inequality is applied.

In div-free FEMs, the second term in (<ref>) vanishes and the final upper bound
in(<ref>) becomes
$\frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x}
+  Ch^{2r+2}  \|u\|^{4}_{L^\infty_t H^{r+1}_x}$.

(3C) $(\text{NL}(2s^n_h-s^{n-1}_h,\phi^{n+1}_h), \phi^{n+1}_h)=0$.

(3D) Estimate of $(\text{NL}(2\phi^n_h-\phi^{n-1}_h,s^{n+1}_h), \phi^{n+1}_h)$.
We need to avoid the differentiation of $(2\phi^n_h-\phi_h^{n-1})$ because there are no terms containing the derivatives of $(2\phi_h^{n+1}-\phi_h^{n})$ on the left of the formula (<ref>), thus unable to generate recursive relations. For this purpose, we rewrite
$(\text{NL}(2\phi^n_h-\phi^{n-1}_h,s^{n+1}_h), \phi^{n+1}_h)$ in the following form by using integration by parts,

\begin{align}
&\,\big| (\text{NL}(2\phi^n_h-\phi^{n-1}_h,s^{n+1}_h), \phi^{n+1}_h) \big|\nonumber \\
=&\, \frac{1}{2} \big|
(((2\phi^n_h-\phi^{n-1}_h)\cdot\nabla) s^{n+1}_h,\phi^{n+1}_h)
- (((2\phi^n_h-\phi^{n-1}_h)\cdot\nabla) \phi^{n+1}_h, s^{n+1}_h) \big| \label{cat1012}\\
\le&\, \frac{1}{2}\|\nabla s^{n+1}_h\|_{L^\infty_x}
\|2\phi^n_h-\phi^{n-1}_h\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x} \nonumber \\
&\qquad\qquad\qquad\qquad\qquad + \frac{1}{2}
\|s^{n+1}_h\|_{L^\infty_x}
\|2\phi^n_h-\phi^{n-1}_h\|_{L^2_x}
\|\nabla \phi^{n+1}_h\|_{L^2_x}
\label{cat1013}
\\
\le&\, C \|\nabla u\|_{L^\infty_t L^\infty_x}
\|2\phi^n_h-\phi^{n-1}_h\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x}
\nonumber \\
&\qquad\qquad\qquad\qquad\qquad + C \|u\|_{L^\infty_t H^2_x} \|2\phi^n_h-\phi^{n-1}_h\|_{L^2_x}
\|\nabla \phi^{n+1}_h\|_{L^2_x}
\label{cat2015}\\
\le&\, \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x}
+ C \|u\|^2_{L^\infty_t H^{r+1}_x}
\|2\phi^n_h-\phi^{n-1}_h\|^2_{L^2_x}
\nonumber \\
&\qquad\qquad\qquad\qquad\qquad + \frac{\nu}{4} \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
+ \frac{C \| u\|^2_{L^\infty_t H^{r+1}_x}}{\nu}
\|2\phi^n_h-\phi^{n-1}_h\|^2_{L^2_x}.
\label{mp003d}
\end{align}

From (<ref>) to (<ref>) then to (<ref>),  the estimates (<ref>) and (<ref>) are employed to obtain $\|\nabla s^{n+1}_h\|_{L^\infty_x}\le C\|\nabla u\|_{L^\infty_t L^\infty_x} \le C \|u\|_{L^\infty_t H^3_x} \le C \|u\|_{L^\infty_t H^{r+1}_x}$ when $r\ge 2$,
and (<ref>) is used to obtain $\|s^{n+1}_h\|_{L^\infty_x}\le \|u\|_{L^\infty_t H^2_x} \le
\|u\|_{L^\infty_t H^{r+1}_x}$. Young's inequality is applied from (<ref>) to (<ref>).

Note when a divergence-free finite element space pair is used,
$(\text{NL}(2\phi^n_h-\phi^{n-1}_h,s^{n+1}_h), \phi^{n+1}_h)
= (((2\phi^n_h-\phi^{n-1}_h)\cdot\nabla) s^{n+1}_h,\phi^{n+1}_h)$, whose upper bound would contain only the first two terms in (<ref>), thus without the viscosity terms.

(3E) $(\text{NL}(2\phi^n_h-\phi^{n-1}_h,\phi^{n+1}_h), \phi^{n+1}_h)=0$.

(3F) Estimate of $(\text{NL}(2u^n-u^{n-1}-u^{n+1}, u^{n+1}), \phi^{n+1}_h )$.
Since $u^{n}$, $u^{n-1}$ and $u^{m+1}$ are exact solutions,
$\nabla\cdot( 2u^n-u^{n-1}-u^{n+1})=0$ almost everywhere in $\Omega$.  Thus, $(\text{NL}(2u^n-u^{n-1}-u^{n+1}, u^{n+1}), \phi^{n+1}_h ) = (((2u^n-u^{n-1}-u^{n+1})\cdot \nabla) u^{n+1}, \phi^{n+1}_h)$.
Because $2u^n-u^{n-1}$ is a second order accurate approximation of $u^{n+1}$, there exists a positive constants  $C_2$ such that for almost every $x\in\Omega$,

\begin{align}
( 2u^n-u^{n-1}-u^{n+1}) (x)
&= C_2 k^2 u_{tt}(\xi_2(x),x),
\end{align}

where $\xi_2(x)\in (t^{n-1}, t^{n+1})$.

\begin{align}
&\, \big| (\text{NL}(2u^n-u^{n-1}-u^{n+1}, u^{n+1}), \phi^{n+1}_h ) \big|
\le \big|
(( C_2 k^2 u_{tt}(\xi_2(x),x) \cdot\nabla ) u^{n+1}, \phi^{n+1}_h)
\big|  \nonumber \\
\le&\,  C_2 k^2 \|\nabla u^{n+1}\|_{L^\infty_x}
\|u_{tt}(\xi_2(x),x)\|_{L^2_x}
\|\phi^{n+1}_h\|_{L^2_x}
\le Ck^2 \|u\|_{L^\infty_t H^3_x}
\|u_{tt}\|_{L^2_x L^\infty_t}
\|\phi^{n+1}_h\|_{L^2_x}
\nonumber \\
\le&\,
\frac{1}{10}\|\phi^{n+1}_h\|^2_{L^2_x}
+ Ck^4  \|u\|^2_{L^\infty_t H^3_x}  \|u_{tt}\|^2_{L^2_x L^\infty_t}.
\label{mp003f}
\end{align}

To reach (<ref>), we have used
$\|\nabla u^{n+1}\|_{L^\infty_x}  \le \|u\|_{L^\infty_t H^3_x}$
due to (<ref>),
\le \|u_{tt}\|_{L^2_x L^\infty_t}$ similar to (<ref>), and Young's inequality.

Step 4. Finalization of error analysis from (<ref>) and upper bounds obtained in Step 2.
Dropping $\|\phi^{n+1}_h - 2\phi^n_h + \phi^{n-1}_h\|^2_{L^2_x}$ on the left of (<ref>) and replacing the last four terms by their upper bounds in (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), we obtain

\begin{align}
&\,\frac{1}{k} \big(
\|\phi^{n+1}_h\|^2_{L^2_x}
+ \|2\phi^{n+1}_h-\phi^n_h\|^2_{L^2_x}
\big)
+ \nu \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\nonumber\\
\le&\,
\frac{1}{k}
\big(
\|\phi^{n}_h\|^2_{L^2_x}
+ \|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}
\big)
+ \|\phi^{n+1}_h\|^2_{L^2_x}
+ \frac{\nu}{2} \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\nonumber \\
C \Big(1+ \frac{1}{\nu} \Big)
\|u\|^2_{L^\infty_t H^{r+1}_x}
\|2\phi^n_h-\phi^{n-1}_h\|^2_{L^2_x}
\nonumber\\
C \Big( h^{2r+2} \|u_t\|^2_{L^\infty_t H^{r+1}_x}
+ \frac{h^{2r} \|p\|^2_{L^\infty_t H^{r}_x}}{\nu}
+ h^{2r}   \|u\|^{4}_{L^\infty_t H^{r+1}_x}
\nonumber \\
&\qquad\quad + k^4  \|u_{ttt}\|^2_{L^2_x L^\infty_t} + k^4 \|u\|^2_{L^\infty_t H^3_x}
\|u_{tt}\|^2_{L^2_x L^\infty_t}
\Big).
\label{eqn_NSu_07}
\end{align}

According to the inequality $\sum_{i=1}^m x_i^2 \le \left( \sum_{i=1}^m x_i \right)^2$ for any positive integer $m$ and positive real numbers $x_i, i=1, \cdots, m$,  the sum of the five terms in the last pair of brackets on the right of (<ref>) is bounded above by $M_2^2$,  where $M_2$ is defined in (<ref>).

Then, by multiplying $k$ on both sides and rearranging terms,  we can rewrite (<ref>) as

\begin{align}
(1-k) \|\phi^{n+1}_h\|^2_{L^2_x}
+ \|2\phi^{n+1}_h-\phi^n_h\|^2_{L^2_x}
\big)
+ \frac{k \nu}{2} \|\nabla\phi^{n+1}_h\|^2_{L^2_x}
\nonumber\\
\le&\,
\big(
\|\phi^{n}_h\|^2_{L^2_x}
+ (1+ kCM_1) \|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}
\big)
+ kC M_2^2,
\label{eqn_NSu_09}
\end{align}

where $M_1$ is defined in (<ref>).

In the case of divergence free FEMs, many terms in (<ref>) are changed based on the above analysis, but the same estimate  (<ref>) holds with  $M_1$  and $M_2$  replaced with  $M^{div0}_1$ defined in (<ref>) and  $M^{div0}_2$ in (<ref>), respectively.

We further denote

\begin{equation}
a^n\triangleq \|\phi^{n}_h\|^2_{L^2_x}
+ \|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}
+ \frac{k\nu}{2}
\|\nabla\phi^{n+1}_h\|^2_{L^2_x}.
\label{def_an}
\end{equation}

It is easy to check that when $k<1$,  the left side of (<ref>) is bounded below by $(1-k ) a^{n+1}$, and
the right side of (<ref>)
is bounded above by
$(1+kCM_1) a^n + kCM_2^2 $.
Then (<ref>)  can be thrown into a recursive relation $a^{n+1}\le \beta a^n + \alpha$,

\begin{align}
\alpha = \frac{kC M^2_2}{1-k},
\qquad
\beta = \frac{ 1 + kCM_1 }{1-k}.
\end{align}

This relation ultimately generates
$a^{n} \le e^{n (\beta-1)} \big( a^1 + \frac{\alpha}{\beta-1} \big)$, which corresponds to

\begin{align}
a^{n} \le
\exp\Big(
\frac{nk (1+CM_1) } { 1-k }
\Big)
\cdot
\Big(
a^1 + \frac{C M^2_2}{1+CM_1 }
\Big).
\end{align}

By using $nk\le T$ for $n\le \lfloor \frac{T}{k} \rfloor$
and  dropping
$\|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}$ in $a^n$, we obtain

\begin{align}
\|\phi^n_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \|\nabla \phi^n_h\|^2_{L^2_x}
\le
\exp{ \Big( \frac{(1+CM_1) T}{1-k}  \Big) }
\cdot
\Big(a^1 + \frac{ CM^2_2}{1+CM_1}
\Big).
\label{error_anal_01}
\end{align}

As for the full error $u^n - u^n_h=(u^n-s^n_h) + \phi^n_h$, we have

\begin{align}
&\, \|u^n-u^n_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \| \nabla(u^n-u^n_h)\|^2_{L^2_x}
\nonumber \\
\le&\,
\|\phi^n_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \|\nabla \phi^n_h\|^2_{L^2_x}
+ \|u^n- s^n_h\|^2_{L^2_x}
+  \frac{k\nu}{2} \| \nabla(u^n-s^n_h)\|^2_{L^2_x}
\\
\le&\,
\|\phi^n_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \|\nabla \phi^n_h\|^2_{L^2_x}
+ Ch^{2r+2} \|u\|^2_{L^\infty_t H^{r+1}_x}
+ C\nu k h^{2r} \|u\|^2_{L^\infty_t H^{r+1}_x}
\label{error_anal_02}
\end{align}

where the inequality (<ref>) is applied in the last step.
Finally, the error estimate (<ref>) is obtained by combining (<ref>) and (<ref>).

The error analysis of the fully implicit convection form is similar to the semi-implicit convection case and the result is summarized below.

Assume the fully implicit convection form $\text{NL}(u^{n+1}_h, u^{n+1}_h)$ is used in(<ref>).
Suppose $u, u_t\in L^\infty_t(H^{r+1}_x)$ with $r\ge 2$, $u_{tt}, u_{ttt}\in L^2_x(L^\infty_t)$ and $p\in L_t^\infty H^r_x$.  Let the time step size $k<1$ and the spatial mesh size $h<1$.  Then for any $2\le n \le \lfloor \frac{T}{k} \rfloor$, there exists a constant $C$ independent of the solution $(u,p)$ and the discretization parameters $k$ and $h$, such that the velocity solution of the scheme (<ref>) and  (<ref>)  satisfies,

(i) under the FI-SKEW convection with non-div-free FEMs,

\begin{align}
&\,\|u^n - u^n_h\|^2_{L^2_x} +
\frac{k\nu}{2} \|\nabla (u^n-u^n_h)\|^2_{L^2_x}
\le Ch^{2r} (h^2+ \nu k) \|u\|^2_{L^\infty_t H^{r+1}_x}\nonumber\\
&\,+\exp{\Big( \frac{T M_3}{1-k  M_3}\Big)} \cdot
\Big(\|\phi^1_h\|^2_{L^2_x}
+ \|2\phi^1_h - \phi^0_h\|^2_{L^2_x}
+ \frac{k\nu}{2} \|\nabla \phi^1\|^2_{L^2_x}
+ \frac{CM^2_4}{M_3}
\Big).
\label{error_anal_FI}
\end{align}

This error estimate is non-robust and has  spatial accuracy order $r$.

under the FI-SKEW  with div-free FEMs, the estimate  (<ref>)  with $M_3$ and $M_4$ replaced with $M^{div0}_3$ and $M^{div0}_4$, respectively. This error estimate is robust with spatial accuracy order $r$.

The constants $M_3$, $M^{div0}_3$, $M_4$, and $M^{div0}_4$ are defined below,

\begin{align}
M_3 =& M^{div0}_3+  \frac{1}{\nu} \|u\|^2_{L^\infty_t H^{r+1}_x},
\quad  M^{div0}_3 =\, 1+ C \|u\|_{L^\infty_t H^{r+1}_x},
\label{def_M3} \\
M_4 =&\,
M^{div0}_4 + \frac{h^r}{\sqrt{\nu}} \|p\|_{L^\infty_t H^r_x},
M^{div0}_4 =
h^{r+1} \|u_t\|_{L^\infty_t H^{r+1}_x}
+ h^r \|u\|^2_{L^\infty_t H^{r+1}_x}
+  k^2 \|u_{ttt}\|_{L^2_x L^\infty_t}.
\label{def_M4M5}
\end{align}

Step 1 and Step 2 are identical to those in the proof of Theorem <ref>.

Step 3. In the FI convection case,  $B(\phi^{n+1}_h) =
(\text{NL}(u^{n+1}_h, u^{n+1}_h), \phi^{n+1}_h)
-(\text{NL}(u^{n+1}, u^{n+1}), \phi^{n+1}_h)$.
We replace $u^j_h$ by $u^j-\eta^j + \phi^j_h$ for $j=n+1, n, n-1$ and obtain

\begin{align}
=&\, -(\text{NL}(u^{n+1},\eta^{n+1}), \phi^{n+1}_h)
-(\text{NL}(\eta^{n+1},s^{n+1}_h), \phi^{n+1}_h)
+(\text{NL}(s^{n+1}_h,\phi^{n+1}_h), \phi^{n+1}_h)
\nonumber \\
\label{hai001}
\end{align}

Below is the estimate of the 5 terms on the right of (<ref>).

(3A) Similar to (3A) in the proof of Theorem <ref>, one gets

\begin{align}
|\text{NL}(u^{n+1},\eta^{n+1}, \phi^{n+1}_h) |
\le&\, \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x} + Ch^{2r} \|u\|^4_{L^2_xH^{r+1}_x},
\label{hai002}\\
|\text{NL}(u^{n+1},\eta^{n+1}, \phi^{n+1}_h) |
\le&\, \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x} + Ch^{2r+2} \|u\|^4_{L^2_xH^{r+1}_x}.
\label{hai003}
\end{align}

(3B) Similar to (3B) in the proof of Theorem <ref>, one gets

\begin{align}
|\text{NL}(\eta^{n+1},s_h^{n+1}, \phi^{n+1}_h) |
\le&\, \frac{1}{10} \|\phi^{n+1}_h\|^2_{L^2_x} + Ch^{2r} \|u\|^4_{L^2_xH^{r+1}_x}.
\label{hai004}
\end{align}

(3C) $(\text{NL}(s^{n+1}_h,\phi^{n+1}_h), \phi^{n+1}_h) =0$.

(3D) Similar to (3D) in the proof of of Theorem <ref>, one gets for non-div-free FEMs

\begin{align}
|\text{NL}(\phi^{n+1}_h,s_h^{n+1}, \phi^{n+1}_h) |
\le&\, C \bigg(  \|u\|_{L^2_xH^{3}_x} + \frac{C\|u\|^2_{L^2_xH^{r+1}_x}}{\nu} \bigg) \|\phi^{n+1}_h\|^2_{L^2_x}
+ \frac{\nu}{4} \|\nabla\phi^{n+1}_h\|^2_{L^2_x}.
\label{hai007}
\end{align}

For div-free FEMs, the viscosity terms on the right of (<ref>) vanish.

(3E) $|\text{NL}(\phi^{n+1}_h, \phi^{n+1}_h, \phi^{n+1}_h) =0$.

Step 4.
After collecting all estimates from Step 2 and 3, we obtain for FI-SKEW with non-divergence free FEMs

\begin{align}
&\,\frac{1}{k} \big(
\|\phi^{n+1}_h\|^2_{L^2_x}
+ \|2\phi^{n+1}_h-\phi^n_h\|^2_{L^2_x}
\big)
+ \nu \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\le
\frac{1}{k}
\big(
\|\phi^{n}_h\|^2_{L^2_x}
+ \|2\phi^{n}_h-\phi^{n-1}_h\|^2_{L^2_x}
\big) \nonumber\\
&+ \frac{\nu}{2} \|\nabla \phi^{n+1}_h\|^2_{L^2_x}
\Big(1+ C \Big(\|u\|_{L^\infty_t H^{r+1}_x}
+ \frac{1}{\nu} \|u\|^2_{L^\infty_t H^{r+1}_x} \Big) \Big)
\|\phi^{n+1}_h\|^2_{L^2_x}
\nonumber\\
C \Big( h^{2r+2} \|u_t\|^2_{L^\infty_t H^{r+1}_x}
+ \frac{h^{2r} \|p\|^2_{L^\infty_t H^{r}_x}}{\nu}
+ h^{2r}   \|u\|^{4}_{L^\infty_t H^{r+1}_x}
+ k^4  \|u_{ttt}\|^2_{L^2_x L^\infty_t}
\Big).
\label{hai009}
\end{align}

Afterwards, using the same process as in Step 4 in the proof of Theorem <ref>, the error estimate (<ref>) is achieved. Specifically,

(i) under FI-SKEW with non-div-free FEMs, the error upper bound (<ref>) has the same accuracy order as shown in  Remark <ref>, which is non-robust and  of spatial order $r$.

(ii) under FI-SKEW with div-free FEMs, the error upper bound (<ref>) with $M_3$ replaced with $M^{div0}_3$ and $M_4$ replaced with $M^{div0}_4$ has the same accuracy order as shown in  Remark <ref>, which is robust and of spatial order $r$.

§ NUMERICAL TESTS

In all the numerical tests, the spatial domain is a unit cube and the regular tetrahedral meshes are used.  Denote the number of subdivisions in each direction as $N$.

The Reynolds number is defined as

\begin{equation}
Re=\frac{L U}{\nu},
\end{equation}

where the length scale $L=1$. The velocity scale $U$ is taken as the maximum initial velocity magnitude. The kinematic viscosity $\nu$ is modulated to give various values of $Re$.

The spatial discretization is chosen as the Taylor-Hood pair of  finite element spaces $P_2/P_1$, that is, $V_h$ is composed of continuous piecewise quadratic polynomials and $Q_h$ continuous piecewise linear polynomials. This finite element pair is not divergence free according to Definition <ref>.
The numerical integration uses a 24-point quadrature rule in [21], which is exact for 6th degree polynomials in a tetrahedron.

To solve the  non-symmetric matrix equations (<ref>), the ILUT preconditioned GMRES method from Yousef Saad's SPARSKIT [25, 24] (https://www-users.cs.umn.edu/s̃aad/software/SPARSKIT/) is adopted.
All the numerical simulations are performed in Michigan State University's High Performance Computing Center (HPCC).

§.§ Stopping criterion  and incompressibility measures

The stopping criterion of the iteration (<ref>) and (<ref>), (<ref>) is set as

\begin{align}
( \|p^{n+1,s+1}_h - p^{n+1,s}_h\|_{\max} <\epsilon
\text{ and }
\|p^{n+1,s+1}_h - p^{n+1,s}_h\|_{L^2_x} <\epsilon )
\text{ or  }
s\ge \text{IterMax}.
\label{stoppingcriterion}
\end{align}

It may be expensive to achieve the full  convergence in every time step  and thus in practice we use an inexact approach where $\epsilon=10^{-2}$ or $10^{-7}$ and $\text{IterMax}=50$.

To track how the divergence free condition is satisfied by the numerical velocity through iterations, we use two measures. The first is a weak incompressibility measure,

\begin{equation}
|\text{weak div}|_{\text{max}} =
|\pi_h(\nabla\cdot u_h^s)|_{\text{max}}
= \max_{q \text{ is a basis function of }Q_h}
|(\nabla\cdot {u}^s_h, q)|,
\label{weak_divfreemeasure}
\end{equation}

where $\pi_h$ is the $L^2$ projection from $V_h$ to $Q_h$, and the basis functions of $Q_h$ are chosen as the continuous, piecewise linear functions which are of value one at one mesh point and zero on other mesh points. This measure has been used by Fortin and Glowinski [13] (Chapter 2) to monitor the convergence of some iterative algorithms for the Stokes  equations.
From the iterative relations (<ref>) and (<ref>), it is found that

\begin{equation}
(\vec{\bar{p}}^{s+1} - \vec{\bar{p}}^{s})_i
= -\frac{1}{k} (\alpha G^{-1} + \rho k M^{-1}) \cdot (\nabla\cdot \bar{u}^s_h, q_i), \quad \forall\, i=1, \cdots, m.
\label{kun1022}
\end{equation}

Since the matrix $(\alpha G^{-1} + \rho k M^{-1})$ is invertible, the weak incompressibility measure provides an equivalent stopping criterion  to $(p^{s+1}-p^s)$ for the iterative scheme.
Therefore, when the iterations converge, the weak incompressibility measure is zero for the velocity limit.

The second measure is the strong incompressibility measure which is defined as the pointwise maximum magnitude of the  strong divergence of velocity,

\begin{equation}
|\text{strong div}|_{\text{max}}
= \max_{{\bf  x}\in \Omega}
|\nabla\cdot {u}^s_h({x})|.
\label{eqn_strong_measure}
\end{equation}

It is easy to see that  these two measures are equivalent in the div-free FEMs.
In the finite dimensional space $\nabla\cdot V_h\subset L^2(\Omega)$,
$\max_{x\in\Omega} |\nabla\cdot u^s_h(x)|$ is equivalent to $\|\nabla\cdot u^s_h\|_{L^2(\Omega)}$.
Using $\nabla\cdot u=0$ for the exact solution $u$ and the inequality (<ref>),
$\|\nabla\cdot u^s_h\|_{L^2(\Omega)} =
\|\nabla\cdot (u^s_h-u)\|_{L^2(\Omega)}\le
\|\nabla (u^s_h-u)\|_{L^2(\Omega)}$.
This correlation suggests the strong incompressibility measure could serve as a  marker of the spuriousness of the numerical solution.

§.§ Numerical tests of iteration convergence at a single time step

In Problem 1,  the exact solution is

\begin{align}
u_1&=\cos(t) g(x)  g'(y) g'(z), \hspace{.38 in} u_2=\cos(t) g'(x)  g(y) g'(z),\\
u_3&=-2 \cos(t) g'(x)  g'(y) g(z), \hspace{.25 in} p=\cos(3t) \sin(2\pi x) \sin(2\pi y) z^3,
\end{align}

where $g(x)=10x^2(1-x)^2$. The computational domain is $\Omega=[0,1]^3$.
Thus,  the boundary value of velocity is zero
and $\max{|u|}\approx 4.6$. In this problem, the meshes are always uniform of resolution $h$.  The Reynolds number of this problem is $Re=\frac{4.6}{\nu}$.

The first group of tests examines the convergence of the projection iterations in one single time step of the scheme (<ref>), (<ref>), and (<ref>)  on Problem 1.   The target is to find the velocity and pressure at $t_2$, where their values at $t_0$ and $t_1$ are chosen as the exact solutions.
All the simulations in this section use the SI-SKEW convection form, which have the similar behaviors as those in the FI-SKEW form. The results for Problem 2 mentioned in Section <ref> is similar and thus not shown here.

The mesh resolution is chosen as $N=80$ and the time step size is $k=10^{-3}$.
The computational results shown in Figure <ref> agree with the normal mode analysis. That is, the Uzawa iterations  reduce the weak incompressibility measure significantly only when the viscosity is large ($\nu\ge 10^{-2}$), and have almost no effects when the viscosity is small ($\nu\le 10^{-4}$). In contrast, the standard projection iterations lead to the fast decay of the weak measure only when the viscosity is small ($\nu\le 10^{-2}$). But the rotational projection iterations give the fast decay of the weak measure for all the viscosity values being studied. A common feature of all these simulations is that the strong incompressibility measure tends to a nonzero steady state,  which implies the velocity limit in this convergence process is not pointwise divergence free.
This is because the Taylor-Hood FEMs are not divergence free.

Weak and strong measures of velocity over iterations of Problem 1 in the time step $t=2\times 10^{-3}$ with SI-SKEW convection, $N=80$, $k=10^{-3}$. First row: weak measure. Second row: strong measure.

§.§ About optimal parameter values of $\alpha$ and $\rho$

The normal mode analysis of Section <ref> shows that the choice of $\alpha=1.5$ and $\rho=\nu$ results in super convergence when the numerical solution is smooth and the convection is absent. When the numerical velocity has only $H^1$ regularity and the convection is present, it is hard to analyze the optimal values of these two parameters. Thus, we use numerical tests to get some clues.

Firstly, we study the effect of parameter values on the iteration convergence at the initial time step $t=0.002$ of Problem 1 with the SI-SKEW convection form and the mesh resolution $N=80$. When $\alpha=1.5$, the convergence is faster as $\rho$ grows from 0 to $2\nu$ when $\nu=1$ and $\nu=0.001$, as shown in Figure <ref>[ab], with an exception at $\nu=1$, $\rho=2\nu$.  A direct comparison between $\nu=1$ and $\nu=0.001$ cases indicates the optimal parameters depend less sensitively on $\rho$ when $\nu$ is small.
On the other hand, when $\nu=0.001$ and $\rho=2\nu$, the convergence increases as $\alpha$ goes from 1.5 to 2.4 or 2.6 (Figure <ref>[c]) and decreases as $\alpha$ continues from 2.6 to 3.0. Thus, the optimal value should be close to 2.5.

Effects of $\alpha$ and $\rho$ on iteration convergence in Problem 1 in the time step $t=2\times 10^{-3}$ with SI-SKEW convection, $N=80$, $k=10^{-3}$.
[a]: $\nu=1$ and $\alpha=1.5$.
[b]: $\nu=10^{-3}$ and $\alpha=1.5$.
[c]: $\nu=10^{-3}$ and $\rho=2\nu$.

Secondly, we use the iterative projection method with FI-SKEW convection to solve Problem 1 from $t=0$ to $t=0.5$ with $Re=4600$, $N=60$, and the stopping criterion is $|\text{weak div}|_{\max}<0.01$. We test the speed of the iterative projection method with some different values of $\alpha\in (1, 3)$ and $\rho\in [\nu, 2\nu]$ and the results are shown in Table <ref>. All these simulations produce the same error as in the case $N=60$ in  Table <ref>. Note in Table <ref>, the simulations of the same $\alpha$ value but different $\rho$ values have the identical solution process, which means the $\rho$ choice is not important here. The optimal value of $\alpha$ seems to be close to 2.5 in this particular example, which is consistent with the result from the initial time step test.

Effects of parameters $\alpha$ and $\rho$ on iterations of FI-SKEW scheme of evolution of Problem 1

$\alpha$     $\rho$     average iterations

1.5     $\nu$, $1.9\nu$, $2\nu$     9.0

1.9     $\nu$, $1.9\nu$, $2\nu$     6.2

2.0     $\nu$, $1.9\nu$, $2\nu$     6.1

2.5     $\nu$, $1.9\nu$, $2\nu$     5.0

2.8     $\nu$, $1.9\nu$, $2\nu$     6.1

The similar tests are done for Problem 2 mentioned in Section <ref>, where the optimal value occurs around $\alpha=2.9$ and $\rho=2\nu$ when $\nu=10^{-5}$ (data not shown). Based on these numerical tests, we find that the optimal values of parameters to attain the fastest iteration convergence are in the range of $(\alpha,\rho) \in (1.5,3)\times (\nu,2\nu)$ and the specific value depends on the specific problem. It is recommended to test the optimal values at the initial time step with a coarse mesh before running a long time evolution.

§.§ Numerical tests on accuracy of iterative projection method of Problem 1

This group of tests investigates the stability and accuracy of  the full scheme (<ref>), (<ref>), and (<ref>) when it is applied to  Problem 1. The simulations run from $t=0$ to $t=0.5$  with various Reynolds numbers starting from $Re=920$, where the convergence is evaluated at time $t=0.5$.  All the simulations in Section <ref>, Section <ref>, and Section <ref> use $\alpha=1.5$, $\rho=\nu$, $k=0.001$.

§.§.§ $Re=920$:  one iteration attains  accurate solutions

In the first case, we choose $Re=920$ and the SI-SKEW convection form. When only one iteration is used, that is, the  conventional rotational projection method is employed, the numerical solutions show the optimal convergence (see Table <ref>): second order accuracy of  both the velocity (in $H^1$ norm) and  the pressure (in $L^2$ norm). Therefore, in this case, more than one iteration is not needed. Notably, both the weak and strong incompressibility measures are small (less than 1) as shown in Figure <ref>.

Errors at $t=0.5$ when $Re=920$,  iteration=1.

N      $\|{ u}_h - { u}\|_{H^1}$     	rate     	$|p_h-p|_{L^2}$	    rate

60	     $0.20\times 10^{-1}$	   	1.86	     $0.15\times 10^{-4}$	   	2.01

70	     $0.15\times 10^{-1}$	   	2.32	     $0.11\times 10^{-4}$	   	2.11

80	     $0.11\times 10^{-1}$     	1.89	     $0.83\times 10^{-5}$	   	2.08

90	     $0.88\times 10^{-2}$	   	2.04	     $0.65\times 10^{-5}$	   	2.12

100      $0.71\times 10^{-2}$	   		     $0.52\times 10^{-5}$

Weak [a] and strong [b] incompressibility measures for Problem 1 with $Re=920$, iteration=1.

§.§.§ $Re=2300$: iterations reduce/remove spurious solutions

In this subsubsection, all the simulations use the SI-SKEW convection.

(1) Only one projection per time step.

When  $Re=2300$, only one projection in every time step leads to large errors when $N=60, 70, 80$, but delivers accurate solutions when $N=90, 100$ (see Table <ref>). These coincide with sizes of the weak and strong incompressibility measures as shown in Figure <ref>: when $N=60, 70, 80$, the measures grow exponentially in time and then stay at high values, but when $N=90, 100$, the measures are small (less than 0.1).

Errors at $t=0.5$ when $Re=2300$, iteration=1.

N      $\|{ u}_h - { u}\|_{H^1}$     	rate     	$|p_h-p|_{L^2}$	    rate

60         $0.12\times 10^{+3}$       11        $0.31\times 10^{+1}$     	 6

70         $0.42\times 10^{+2}$       33        $0.12\times 10^{+1}$     	35

80         $0.51\times 10^{0}$       33        $0.10\times 10^{-1}$ 	    	62

90         $0.94\times 10^{-2}$       2.14      $0.66\times 10^{-5}$       1.90

100        $0.75\times 10^{-2}$                $0.54\times 10^{-5}$

Weak [a] and strong [b] incompressibility measures of Problem 1 when $Re=2300$, iteration=1.

(2) Using fixed numbers of  iterations per time step.

To see the effects of iterative projections, we  focus on the mesh $N=60$ with a fixed number of iterations per time step:  iteration=1, 2, 4, or 6.  Based on Table <ref>, the errors decrease when more iterations are used and when iteration=6, the errors are optimal in the given time and spatial discretizations.
The plottings of the numerical solutions in Figure <ref> reveal that there are strongly unphysical oscillations  when only one iteration is used, but the oscillations reduce when more iterations are added and completely vanish when iteration=6. Correspondingly, both the weak and strong incompressibility measures increase to large values (more than 1000) with respect to time when iteration=1,2,4, but small (less than 0.1) when iteration=6, as illustrated in Figure <ref>. Especially the strong incompressibility measure shares the same tendency as the $H^1$ error of the velocity in the fixed $V_h$ space ($h=1/60$) when the iteration number increases, which is an evidence to Remark <ref> that the strong incompressibility measure is a marker of spuriousness.

Errors at $t=0.5$ when  $Re=2300$, $N=60$,  and a fixed number of iterations, “iter”,  is applied in each time step.

$Re$     N     iter      $\|{ u}_h - { u}\|_{H^1}$      	$|p_h-p|_{L^2}$

2300     60     1     $0.12\times 10^{+3}$        $0.31\times 10^{+1}$

2300     60     2     $0.43\times 10^{+2}$        $0.13\times 10^{+1}$

2300     60     4     $0.29\times 10^{+1}$        $0.26\times 10^{-1}$

2300     60     6     $0.23\times 10^{-1}$        $0.15\times 10^{-4}$

Numerical solution $u_3$ on the plane $z=\frac{59}{60}$ of Problem 1 at t=0.5 when $Re=2300$ and $N=60$. Iteration=1 in [a], 2 in [b], 4 in [c], 6 in [d].

Weak [a] and strong [b] incompressibility measures when $Re=2300$, $N=60$, and iteration=1, 2, 4, 6.

(3) Using stopping criterion in formula  (<ref>)

One apparent issue of using a fixed number of iterations in every time step is that the smallest number of iterations sufficient for removing spurious modes is unknown before simulations. Therefore,  we switch the stopping criterion (<ref>) with $\epsilon=10^{-2}$  and re-run the simulations of Problem 1 with $Re=2300$. The optimal convergence is recovered as shown in Table <ref>. Both the weak and strong incompressibility measures are small, according to Figure <ref>.

Errors at $t=0.5$ when $Re=2300$ and the stopping criterion is (<ref>) with $\epsilon=10^{-2}$. “average iter” is the average iteration in all the time steps.

N     $\|{ u}_h - { u}\|_{H^1}$      	$|p_h-p|_{L^2}$     average iter

60      $0.23\times 10^{-1}$     $0.15\times 10^{-4}$      5.2

70      $0.16\times 10^{-1}$     $0.11\times 10^{-4}$      3.5

80      $0.12\times 10^{-1}$     $0.83\times 10^{-5}$      1.7

90      $0.94\times 10^{-2}$     $0.66\times 10^{-5}$      1

100     $0.75\times 10^{-2}$     $0.53\times 10^{-5}$      1

Weak [a] and strong [b] incompressibility measures of Problem 1 when $Re=2300$ and the stopping criterion is (<ref>) with $\epsilon=10^{-2}$.

§.§.§ $Re=4600$: semi-implicit vs fully-implicit schemes

(1) SI-SKEW-based iterative projection method.

When the Reynolds number increases to $Re=4600$,  we first employ the SI-SKEW convection and set the stopping criterion as (<ref>) with $\epsilon=10^{-7}$. The numerical solutions show a strange behavior. The standard exact solution of $u_1$ of Problem 1 is zero on the plane $z=0.5$ at all time, but the numerical solution is nonzero, as shown in Figure  <ref>. This solution is spurious  but the magnitude of oscillation is much smaller than what are presented in Figure <ref> when iteration=1. In addition, when the mesh is refined from $N=60$ to $N=100$, the numerical solution converges to that at $N=100$ for both the velocity and pressure in $L^2$ norm, as shown in Table <ref>.
The velocity error with respect to the exact solution and the strong incompressibility measure are shown in Figure <ref>[ab]. Both surge suddenly from less than one to some large values at certain time steps between $t=0.2$ and $t=0.4$ for all these simulations.  This further confirms the role of the strong measure as the marker of the spuriousness of numerical solutions.

Numerical solutions of $u_1$ on the $z=0.5$ plane of Problem 1  at $t=0.25, 0.3, 0.5$. Parameters: $Re=4600$, $N=100$, and $k=10^{-3}$. The standard exact solution $u_1$ is zero on this plane. The SI-SKEW convection is used.

Difference between the numerical solution with mesh resolution $N<100$ and that of $N=100$. $Re=4600$. The SI-SKEW convection is used. Stopping criterion is (<ref>) with $\epsilon=10^{-7}$.

N     $\|u_{N} - u_{100}\|_{L^\infty(0,0.5;L^\infty\Omega)}$
$\|u_{N} - u_{100}\|_{L^2(0,0.5;L^2\Omega)}$
$||p_{N} - p_{100}||_{L^2(0,0.5;L^2\Omega)}$
average iter

60         $0.25769\times 10^{+1}$
$0.19696\times 10^{0}$
$0.17290\times 10^{0}$     15.0

70         $0.25764\times 10^{+1}$
$0.15761\times 10^{0}$
$0.13233\times 10^{0}$     17.7

80         $0.26592\times 10^{+1}$
$0.11724\times 10^{0}$
$0.95676\times 10^{-1}$     21.4

90         $0.23662\times 10^{+1}$
$0.71958\times 10^{-1}$
$0.57587\times 10^{-1}$     27.9

96         $0.29156\times 10^{+1}$
$0.44540\times 10^{-1}$
$0.28873\times 10^{-1}$     25.0

97         $0.19185\times 10^{+1}$
$0.41189\times 10^{-1}$
$0.26961\times 10^{-1}$     25.6

98       $0.19823\times 10^{+1}$
$0.40630\times 10^{-1}$
$0.27062\times 10^{-1}$     26.4

99         $0.21280\times 10^{+1}$
$0.35669\times 10^{-1}$
$0.21571\times 10^{-1}$     25.5

[a,b]: velocity error [a] and strong incompressibility measure [b] of SI-SKEW with $Re=4600$. [c]: strong incompressibility measure of FI-SKEW when $Re=4600$.
The stopping criterion is (<ref>) with $\epsilon=10^{-2}$ in [abc].

The direct reason for this false numerical solution for the iterative SI-SKEW scheme is that at a certain time step, the weak incompressibility measure cannot converge to zero given sufficient iterations although this measure can be very small.  For example, when $N=100$, the weak incompressibility measure is stuck at $1.1067\times 10^{-7}$ at iteration steps from 22 to 50 at time step $t=0.219$, the moment when the strong incompressibility measure suddenly rises to a high value in Figure <ref>[b]. Therefore, the iterative  projection method fails to converge at this time step.

(2) FI-SKEW-based iterative projection method

When the FI-SKEW form is adopted and  the stopping criterion is (<ref>) with $\epsilon=10^{-2}$,  the correct solution is recovered as shown in Table <ref>. This suggests the fully implicit scheme is more robust than the semi-implicit scheme under the same setting and a possible reason is that the fully implicit scheme preserves better the nonlinear structure of the original Navier-Stokes equations.

A side note is that although the acceleration technique mentioned in Appendix <ref> can speed up the iteration convergence in many simulations including the cases $N=60, 80, 90$ in Table <ref>,  it is not guaranteed to always produce convergence. For instance, when $N=70$, the iteration using this acceleration blows up at time $t=0.489$. More studies are required to uncover the conditions to apply the acceleration techniques.

Errors when $Re=4600$ and the stopping criterion is (<ref>) with $\epsilon=10^{-2}$. The convection is FI-SKEW.

N     $||{ u}_h - { u}||_{H^1}$      rate     	$|p_h-p|_{L^2}$	    rate      average iter

60      $0.29\times 10^{-1}$     2.4      $0.18\times 10^{-4}$     2.6     9.0

70      $0.20\times 10^{-1}$     2.2      $0.12\times 10^{-4}$     1.4     6.5

80      $0.15\times 10^{-1}$     2.6      $0.10\times 10^{-4}$     3.3     5.1

90      $0.10\times 10^{-1}$     2.1      $0.68\times 10^{-5}$     0.6     3.9

100     $0.88\times 10^{-2}$              $0.64\times 10^{-5}$             2.4

§.§.§ Some higher Reynolds numbers

This set of simulations applies the iterative projection method with the FI-SKEW convection term to handle some high Reynolds numbers $Re=5\times 10^3$, $10^4$, $2\times 10^4$, $4\times 10^4$ of Problem 1,  and the stopping criterion is (<ref>) with $\epsilon=10^{-2}$. The computational results are shown in Figure <ref>.
Both the $L^2$ and $H^1$ errors of velocity are optimal and increase when the Reynolds number is larger, which is consistent with the non-robust error estimates of this non-div-free FEM. The average iterations per time step also increase with the Reynolds numbers and are below 9 in these simulations.

Simulation results for Problem 1 with FI-SKEW convection and the stopping criterion is (<ref>) with $\epsilon=10^{-2}$. [ab]: $L^2$ and $H^1$ errors of velocity. [c]: average iteration numbers per time step.

§.§.§ Comparison with iterative IMEX projection scheme

Following [3],  we use a BDF2-IMEX scheme, where the only difference with the method proposed in this work ((<ref>), (<ref>), and (<ref>))  is that the equation (<ref>) is replaced by

\begin{align}
1.5 ( u_h^{n+1,s}, v_h)
&+ k\nu (\nabla u_h^{n+1,s}, \nabla  v_h)
- k (\nabla\cdot v_h, p^{n+1,s}_h) \notag\\
&=(2u^n_h -0.5u^{n-1}_h, v_h)
+ k ( f^{n+1}, v_h) -k ((w^{n+1}_h\cdot\nabla) w^{n+1}_h, v_h),
\label{eqn_u_IMEX}
\end{align}

where $w^{n+1}_h=2 u^n_h - u^{n-1}_h$.
That is, the diffusion term $\nu (\nabla u_h^{n+1,s}, \nabla  v_h)$ is implicit but the convection term $((w^{n+1}_h\cdot\nabla) w^{n+1}_h, v_h)$ is completely explicit.
Using the same proof for Theorem <ref>, with $f$ replaced with $f-(w^{n+1}_h\cdot\nabla)w^{n+1}_h$ in (<ref>), it is easy to see that when $s\to\infty$, the scheme converges to the following system,

\begin{align}
1.5 ( u_h^{n+1}, v_h)
&+ k\nu (\nabla u_h^{n+1}, \nabla  v_h)
- k (\nabla\cdot v_h, p^{n+1}_h)
=(2u^n_h -0.5u^{n-1}_h, v_h)\notag\\
&+ k ( f^{n+1}, v_h) - k ((w^{n+1}_h\cdot\nabla) w^{n+1}_h, v_h),
\quad \forall v_h\in V_h,
\label{IMEX_eqn_u_skew_conv}\\
(\nabla \cdot u^{n+1}_h, q_h) &= 0,
\quad \forall q_h\in Q_h.
\label{IMEX_eqn_p_skew_conv}
\end{align}

For this scheme, it is impossible to prove an energy estimate similar to Theorem <ref>. Indeed, letting  $v_h=u^{n+1}$ in (<ref>) and $q_h=p^{n+1}_h$ in (<ref>), we obtain

\begin{align}
&\,\|u^{n+1}_h\|^2_{L^2_x} + \|2u^{n+1}_h - u^n_h\|^2_{L^2_x} + \|u^{n+1}_h -2u^n_h + u^{n-1}_h\|^2_{L^2_x}
+ 4k\nu \|\nabla u^{n+1}_h\|^2_{L^2_x} \nonumber \\
=&\, \|u^{n}_h\|^2_{L^2_x} + \|2u^{n}_h-u^{n-1}_h\|^2_{L^2_x} + 4k (f^{n+1}, u^{n+1}_h)
- 4k ((w^{n+1}_h\cdot\nabla) w^{n+1}_h, u^{n+1}_h).
\label{cat3001}
\end{align}

Following a process similar to that used to derive (<ref>) and (<ref>), we can get

\begin{align}
| 4k ((w^{n+1}_h\cdot\nabla) w^{n+1}_h, u^{n+1}_h)|
\le&\, 4k \|w^{n+1}_h\|_{L^4_x} \|u^{n+1}_h\|_{L^4_x}
\|\nabla w^{n+1}_h\|_{L^2_x} \notag\\
\le&\, 8k \|w^{n+1}_h\|^{\frac14}_{L^2_x}
\|\nabla w^{n+1}_h\|^{\frac34}_{L^2_x}
\|u^{n+1}_h\|^{\frac14}_{L^2_x}
\|\nabla u^{n+1}_h\|^{\frac34}_{L^2_x}
\|\nabla w^{n+1}_h \|_{L^2_x}.
\label{cat3000}
\end{align}

By the Poincaré inequality,
$\|w^{n+1}_h\|_{L^2_x} \le C\|\nabla w^{n+1}_h\|_{L^2_x}$
$\|u^{n+1}_h\|_{L^2_x} \le C\|\nabla u^{n+1}_h\|_{L^2_x}$. Then (<ref>) becomes

\begin{eqnarray}
| 4k ((w^{n+1}_h\cdot\nabla) w^{n+1}_h, u^{n+1}_h)|
\le kC \|\nabla (2u^n_h-u^{n-1}_h)\|^2_{L^2_x}
\| \nabla u^{n+1}_h \|_{L^2_x}.
\label{cat3002}
\end{eqnarray}

This nonlinear term makes it impossible to derive a recursive relation from (<ref>).

Similarly, the error estimate for (<ref>) and (<ref>), by following the same process as in the proof of Theorem <ref>, generates one term
$(((2\phi^n_h-\phi^{n-1}_h)\cdot\nabla) (2\phi^n_h-\phi^{n-1}_h), \phi^{n+1}_h)$. Using the same argument as in (<ref>) and (<ref>), we can obtain its upper bound as $C \|\nabla (2\phi^n_h-\phi^{n-1}_h)\|^2_{L^2_x} \|\nabla\phi^{n+1}_h\|_{L^2_x}$. With the same reason for the stability, it is also impossible to derive a recursive relation for the error analysis.

The numerical simulations further suggest this iterative IMEX projection scheme is not so stable as the proposed scheme
((<ref>), (<ref>), and (<ref>)). When $Re=920$ and $k=0.001$, this scheme blows up around $t=0.1$ (see Figure <ref>).
In contrast, the iterative SI-SKEW scheme obtains the optimal convergence with just one iteration in each time step according to Section <ref>.

Weak and strong incompressibility measures when iterative IMEX method is used in Problem 1 with $Re=920$ and $k=0.001$. The stopping criterion is (<ref>) with $\epsilon=10^{-7}$.

§.§ Problem 2: iterative projection method in lid driven cavity flow

Problem 2  is a classic lid driven cavity problem.  The domain is $\Omega=(-0.5,0.5)^3$  and the velocity is ${u} = (0,1,0)$ on the wall $x=-0.5$ and zero on the other walls. The initial velocity is set as zero everywhere except on the sliding wall.  The meshes use Gauss-Lobatto points (e.g. [2]), where
$x_i = 0.5 \cos(\frac{i\pi}{N})$,
$y_j = 0.5\cos(\frac{j\pi}{N})$,
$z_l = \frac{l\pi}{N}-0.5$,
$i,j,l=0,\cdots, N$,
where $N$ is the number of subdivisions in each direction.
The Reynolds number of this problem is $Re=\frac{1}{\nu}$.

We run the simulations with the SI-SKEW convection form,  $\nu=0.001$ or $Re=1000$,  $\alpha=1.5$, $\rho=\nu$, and $k=0.001$.
The simulations are run until the velocity reaches a steady state. Since the exact solution is not available, the solution at $t_1$ is taken the same as $t_0=0$. The solution at $t_2$ is obtained by running the scheme until it converges.

When $N=40$, the conventional rotational projection method shows very chaotic velocity and pressure fields at $t=1$, as shown in Figure <ref>. But when two projection iterations are used in each time step, the velocity field looks very smooth (Figure <ref>).

Velocity[a] and pressure[b] on the x-y plane of $z=0$, at time $t=1$. $Re=1000$, iteration=1, max speed=22.79.

Velocity[a] and ressure[b] on the x-y plane of  $z=0$  at time $t=1$. $Re=1000$, iteration=2, max speed=1.

The accuracy of the numerical solution is confirmed by the  agreement with the results of [2]  as shown in Figure <ref> when the mesh is  refined from $N=20$ to $N=40$ and iteration=2 in each time step.  In [2], the standard projection scheme with BDF2 scheme in time and spectral method in space is used. More importantly, they construct a special solution $u_c$ addressing the edge singularity to accompany the numerical solution. Therefore, without introducing the prescribed singular solution, the proposed iterative  projection method obtains the correct solution with just two projections in each time step.

Comparison of normalized velocities, $v/2$ and $u/2$ (the whole velocity vector is $(u,v,w)$) with those in [2],  on the center lines $(x, 0, 0)$ and $(0, y, 0)$ on the
x-y plane $z=0$ in the cubic cavity flow.
The square symbols are extracted from figures of
[2].  The solid lines are from the iterative projection method of this work.
[a]: $N=20$. [b]: $N=40$.
The Reynolds number $Re = 1000$. The time is $t=20$.

§ SUMMARY

This work introduces an iterative projection method for the  Navier-Stokes equations, whose methodology is to iterate the classic projection method in every time step with a proper convection form to guarantee stability. The benefits are tremendous. First, when fully convergent, it produces the weakly divergence free velocity (and almost everywhere divergence free velocity on divergence free finite element methods). Second, the stability and error estimates are theoretically established when the projection iterations are convergent in every time step. Third, this method can handle high Reynolds numbers.  When the finite element spaces are divergence free, this method is robust, that is, the error estimate of velocity in $L^2$ norm is independent of the viscosity or Reynolds numbers. When the projection iterations are not fully convergent, the numerical experiments show that the method can still preserve stability and accuracy for high Reynolds numbers. All these properties are absent in the classic projection method.  At the same time, this new approach keeps an essential property of the classic projection method:  the splitting of the velocity and pressure fields, which results in the  splitting of a large saddle point system into small elliptic problems.

A  condition to warrant the iteration convergence with two parameters is derived, as shown in Theorem <ref> and
Theorem <ref>. A range of the parameter values to achieve the fast iteration convergence is given in Section <ref>.
Two incompressibility measures are proposed to track the performance of numerical solutions. The weak measure monitors whether the iterations are convergent or not and the strong measure indicates whether the numerical solution is spurious or not.  These two measures are equivalent for the div-free FEMs but not for the non-div-free FEMs. Under the non-div-free FEMs, if the weak measure is small but strong measure is too large, the numerical solution would be non-physical.
In this work, both the semi-implicit and the fully implicit skew-symmetric convection forms are studied and tested.  The numerical simulations show that the fully implicit scheme is more stable and accurate and efficient than the semi-implicit form.