Unconditional long-time stability of a velocity–vorticity method for the 2D Navier–Stokes equations

Abstract

We prove unconditional long-time stability for a particular velocity–vorticity discretization of the 2D Navier–Stokes equations. The scheme begins with a formulation that uses the Lamb vector to couple the usual velocity–pressure system to the vorticity dynamics equation, and then discretizes with the finite element method in space and implicit–explicit BDF2 in time, with the vorticity equation decoupling at each time step. We prove the method’s vorticity and velocity are both long-time stable in the $L^2$ and $H^1$ norms, without any timestep restriction. Moreover, our analysis avoids the use of Gronwall-type estimates, which leads us to stability bounds with only polynomial (instead of exponential) dependence on the Reynolds number. Numerical experiments are given that demonstrate the effectiveness of the method.

1 Introduction

The paper addresses long-time stability of numerical methods for the two-dimensional Navier–Stokes system describing the motion of incompressible Newtonian fluids:

$$\begin{aligned} \frac{\partial {u}}{\partial t} - \nu \Delta {u} + ({u} \cdot \nabla ) {u} + \nabla p&= {f}, \nonumber \\ \text {div} \; {u}&= 0, \end{aligned}$$

(1.1)

where ${u}={u}({x},t)$ denotes a velocity vector field, $p=p({x},t)$ is the pressure, and $f=f({x},t)$ represents (given) external forcing. The solution to (1.1) is well-known (see [5]) to be smooth for all time in the periodic setting, that is, the domain $\Omega$ is a 2D torus $\mathbb {T}^2$, all functions have mean zero over the torus, and the forcing term f is smooth. Moreover, the solution of (1.1) is long-time stable, in the sense that the norms $\Vert u\Vert _{L^2(\Omega )}$ and $\Vert u\Vert _{H^1(\Omega )}$ are bounded uniformly in time for $f\in L^\infty (\mathbb {R}_+,L^2(\Omega ))$ and initial value $u_0\in H^1(\Omega )$, $\int _\Omega u_0=0$. The long-time stability is a key property of (1.1) if one is interested in simulation of a large time scale phenomenon or recovering long term statistics, as commonly the case for simulation of flows with large Reynolds’ numbers, weather prediction, or climate modeling. Therefore, it is of practical interest to design numerical methods for (1.1) which inherit this important property. It is also interesting to explore to what extent popular numerical approaches to (1.1) are long-time stable.

The topic of long-time stability and error control for numerical methods for the Navier–Stokes equations is not new in the literature. Heywood and Rannacher [13, 14] proved uniform in time stability and error estimate in the energy norm for a Crank–Nicolson Galerkin method applied to 3D Navier–Stokes system, assuming the solution of the initial boundary value problem is stable. Simo and Armero [24] examined the long-time stability in the energy norm of several time integration algorithms, including coupled schemes and fractional step/projection methods. More recent studies include the papers [1, 8, 25–27]. The work of Tone and Wirosoetisno [25, 26] proved uniform in time bounds on $\Vert u(t_n)\Vert _{L^2(\Omega )}$ and $\Vert \nabla u(t_n)\Vert _{L^2(\Omega )}$ for implicit Euler and Crank–Nicolson methods. These bounds are subject to restrictions on time step in terms of $\nu$ and a spatial discretization parameter. Badia et al. showed in [1] that $\nabla u \in L^\infty (0,\infty ;L^2(\Omega ))$ for a solution to spatially discretized equations (1.1). First and second order semi-explicit time discretization methods for (1.1) written in vorticity–stream function formulation were studied by Wang and co-workers [8, 27]. Both papers consider spectral discretization in space, and prove long-time stability bounds for the enstrophy and the $H^1$-norm of the vorticity, again all subject to a time step restriction of the form $\Delta t\le c\,Re^{-1}$. Thus, despite progress, the current understanding of the long-time behavior of numerical methods for (1.1) is far from being full: only a few studies address uniform in time error estimates for vorticity or velocity gradient, time step restrictions are common in the analyses, and semi-discrete methods are often treated rather than full discretizations. Moreover, to our knowledge, all proofs of long-time numerical stability bounds for vorticity and the gradient of velocity, invoke a variant of the discrete Gronwall lemma, which results in the dependence of the bounds on the Reynolds number of the form $O(\exp (c^2 Re))$ or even $O(\exp (c^2 Re^2))$. Although being time independent, such bounds are not very practical for higher Reynolds number flows; see [16] for a discussion and an effort to improve numerical stability and error estimates dependence on Re number, but only locally in time.

In this paper, we prove unconditional long-time stability of a fully discrete numerical method for (1.1): for $f\in L^\infty (0,\infty ;H^1(\Omega ))$ we prove uniform in time estimates for the kinematic energy, enstrophy, as well as the $L^2$ norms of velocity gradient and vorticity gradient of a discrete system. A finite element method is used for the spatial discretization, and both first and second order time stepping semi-implicit (linear at each time step) schemes are studied. The stability bounds are unconditional, i.e., no time step restrictions are imposed. Furthermore, our analysis does not rely on any Gronwall type estimate, which allows us to avoid exponential dependence of stability bounds on the Reynolds number. In the present analysis, the dependence is polynomial. Our analysis reveals that the polynomials degree can be significantly lowered at the expense of logarithmic dependence on the spatial mesh size.

The results of the paper systematically exploit the relationship between the vorticity and velocity of the Navier–Stokes system by considering the vorticity dynamics equation and writing the inertia in the momentum equation in the form of Lamb vector. For $w=\nabla \times u$ and $P=\frac{1}{2}|u|^2+p$, we reformulate (1.1) as:

$$\begin{aligned} \frac{\partial {u}}{\partial t} - \nu \Delta {u} + w\times {u} + \nabla P= & {} {f}, \nonumber \\ \text {div} \; {u}= & {} 0,\\ \frac{\partial {w}}{\partial t} - \nu \Delta {w} + ({u} \cdot \nabla ) {w}= & {} \nabla \times {f},\nonumber \end{aligned}$$

(1.2)

where $w\times u:= \left[ -u_2 w, \, u_1 w \right] ^T$. Vorticity plays a fundamental role in fluid dynamics, and studying properties of (1.1) through the vorticity equation is a well established approach in the Navier–Stokes theory, see, e.g., [6, 18]. It is also not uncommon in numerical analysis to design numerical methods based on the vorticity equation, e.g., [7, 10]. For numerical methods, standard closures for the vorticity equations are obtained either in vorticity–stream function variables or with the help of the vector Poisson equation, $\Delta u=-\nabla \times w$. However, recent papers [17, 21] have demonstrated numerical advantages of complementing the vorticity equation with the velocity dynamic equation as in (1.2). Thus, Eq. (1.2) will be the departure point in the present analysis.

The rest of the paper is organized as follows. Section 2 gathers necessary definitions and preliminary results for the analysis that follows. In Sect. 3, we introduce a first order time stepping method and prove its long-time stability with respect to the velocity and vorticity $H^1$ norms. Section 4 introduces a second order method based on BDF2 time discretization. We extend the long-time stability results for this method by taking care of some extra technical details. Since the numerical scheme is non-standard, we also provide with our analysis a series of numerical experiments for a 2D flow past a bluff object. The results of the experiments are presented in Sect. 5.2, and they illustrate the long-time stability and the performance of the method.

We finish the introduction with the following remark. Most of our stability analysis is restricted to the 2D case and, due to the current lack of understanding of the long time behavior of 3D Navier–Stokes solutions, we cannot say to what an extend the results remain valid in 3D. However, the numerical approach studied here has a straightforward extension to 3D, and relying on a past experience, we believe that numerical methods which are physically consistent and computationally efficient for 2D problems are commonly found to be also advantageous for solving 3D Navier–Stokes equations.

2 Notation and preliminaries

We consider a domain $\Omega =(0,2\pi )^2\subset \mathbb {R}^2$, and we restrict this study to the case of periodic boundary conditions. We note that our stability analysis also holds for the case of full Dirichlet velocity and vorticity boundary conditions.

We use the notation $(\cdot ,\cdot )$ and $\Vert \cdot \Vert$ for the $L^2(\Omega )$ inner product and norm, respectively. All other norms will be clearly labeled with subscripts.

The natural velocity and pressure spaces in the periodic setting for the Navier–Stokes equations are

$$\begin{aligned} X&:= H^1_{\#}(\Omega )^2 = \left\{ v\in H^1_{loc}(\mathbb {R})^2,\ v \text{ is } 2\pi \text{-periodic } \text{ in } \text{ each } \text{ direction }, \ \int _{\Omega } v \ dx=0\right\} , \\ Q&:= L^2_{\#}(\Omega ) = \left\{ q\in L^2_{loc}(\mathbb {R})^2,\ q \text{ is } 2\pi \text{-periodic } \text{ in } \text{ each } \text{ direction }, \ \int _{\Omega } q \ dx=0 \right\} . \end{aligned}$$

In two dimensions, vorticity is considered as a scalar, and we define vorticity space as

$$\begin{aligned} Y := H^1_{\#}(\Omega ) = \left\{ v\in H^1_{loc}(\mathbb {R}),\ v \text{ is } 2\pi \text{-periodic } \text{ in } \text{ each } \text{ direction }, \ \int _{\Omega } v \ dx=0 \right\} . \end{aligned}$$

For the discrete setting, we assume $\tau _h$ is a regular, conforming triangulation of $\Omega$ which is compatible with periodic boundary conditions. Let $(X_h,Q_h)\subset (X,Q)$ be inf-sup stable velocity–pressure finite element spaces, $Y_h\subset Y$ be the discrete vorticity space, all defined as piecewise polynomials on $\tau _h$.

The discretely divergence-free subspace will be denoted by

$$\begin{aligned} V_h := \{ v_h \in X_h,\ (\nabla \cdot v_h,q_h)=0\quad \forall q_h\in Q_h \}. \end{aligned}$$

The dual space of $V_h$ is denoted by $V_h^*$ with norm $\Vert \cdot \Vert _{V_h^*}$.

We will utilize in our analysis discrete analogues of the Laplacian operator. Define $\Delta _h$ to be the discrete Laplacian operator on $Y_h$: Given $\phi \in H^1(\Omega )$, $\Delta _h \phi \in Y_h$ satisfies

$$\begin{aligned} (\Delta _h \phi ,v_h) = -(\nabla \phi ,\nabla v_h) \quad \forall v_h\in Y_h. \end{aligned}$$

Define $A_h$ to be a discretely divergence-free Laplace operator, often referred to as a Stokes operator by: given $\phi \in H^1(\Omega )$, $A_h\phi \in V_h$ satisfies

$$\begin{aligned} (A_h \phi ,v_h) = (\nabla \phi ,\nabla v_h) \quad \forall v_h\in V_h, \end{aligned}$$

or equivalently,

$$\begin{aligned} (A_h \phi ,v_h) - (l_h,\nabla \cdot v_h) + (\nabla \cdot A_h\phi ,q_h) = (\nabla \phi _h,\nabla v_h) \quad \forall (v_h,q_h) \in (X_h,Q_h), \end{aligned}$$

where $l_h$ is a Lagrange multiplier used so the divergence constraint does not overdetermine the system. Restricted to $X_h$ and $V_h$, respectively, the linear operators $-\Delta _h$ and $A_h$ are self-adjoint and positive definite. In this case, $(-\Delta _h)^{\frac{1}{2}}:X_h\rightarrow X_h$ and $A_h^{\frac{1}{2}}:V_h\rightarrow V_h$ are well defined and will be used in the paper.

The Poincare inequality will be used heavily throughout: there exists $\lambda$, dependent only on $\Omega$, satisfying

$$\begin{aligned} \Vert \phi \Vert \le \lambda \Vert \nabla \phi \Vert \quad \forall \phi \in X. \end{aligned}$$

An immediate consequence on the Poincare inequality and the definition of discrete Stokes and Laplace operators is that the following bounds hold

$$\begin{aligned} \Vert \nabla v_h \Vert&\le \lambda \Vert A_h v_h \Vert \quad \forall v_h \in V_h, \\ \Vert \nabla z_h \Vert&\le \lambda \Vert \Delta _h z_h \Vert \quad \forall z_h \in Y_h. \end{aligned}$$

We recall the following discrete Agmon inequalities, which are also consequences of discrete Gagliardo-Nirenberg estimates, see [12, p. 298]:

$$\begin{aligned} \Vert v_h \Vert _{L^{\infty }}&\le C \Vert v_h \Vert ^{1/2} \Vert A_h v_h \Vert ^{1/2} \quad \forall v_h\in V_h, \end{aligned}$$

(2.1)

$$\begin{aligned} \Vert z_h \Vert _{L^{\infty }}&\le C \Vert z_h \Vert ^{1/2} \Vert \Delta _h z_h \Vert ^{1/2} \quad \forall z_h\in Y_h, \end{aligned}$$

(2.2)

where C is independent of h. The discrete Sobolev inequality (proven in [9]),

$$\begin{aligned} \Vert \nabla \phi _h \Vert _{L^4} \le C \Vert \nabla \phi _h \Vert ^{1/2} \Vert \Delta _h \phi _h \Vert ^{1/2} \quad \forall \phi _h\in X_h, \end{aligned}$$

(2.3)

again with C independent of h, allows us to prove the following lemma.

Lemma 2.1

For every $z_h\in Y_h$, there exists a constant C, independent of h, satisfying

$$\begin{aligned} \Vert \nabla z_h \Vert _{L^3} \le C \Vert z_h \Vert ^{1/3} \Vert \Delta _h z_h \Vert ^{2/3} \quad \forall z_h\in Y_h. \end{aligned}$$

(2.4)

Proof

By Hölder’s inequality,

$$\begin{aligned} \Vert \nabla z_h \Vert ^3_{L^3} \le \Vert \nabla z_h \Vert \Vert \nabla z_h \Vert ^2_{L^4}, \end{aligned}$$

and thus using (2.3) provides the bound

$$\begin{aligned} \Vert \nabla z_h \Vert ^3_{L^3} \le C \Vert \nabla z_h \Vert ^2 \Vert \Delta _h z_h \Vert . \end{aligned}$$

Since $\Vert \nabla z_h \Vert ^2 = (\nabla z_h,\nabla z_h)=-(z_h,\Delta _h z_h) \le \Vert z_h \Vert \Vert \Delta _h z_h \Vert$, the estimate becomes

$$\begin{aligned} \Vert \nabla z_h \Vert ^3_{L^3} \le C \Vert z_h \Vert \Vert \Delta _h z_h \Vert ^2. \end{aligned}$$

Taking cube roots of both sides completes the proof. $\square$

Define the skew-symmetric trilinear operator $b^*:X_h\times Y_h \times Y_h \rightarrow \mathbb {R}$ by

$$\begin{aligned} b^*(u,w,\chi )=(u\cdot \nabla w,\chi ) + \frac{1}{2} ((\nabla \cdot u)w,\chi ). \end{aligned}$$

We will exploit the property that $b^*(u,w,w)=0$ in our analysis of the vorticity equation.

3 Backward Euler

We first consider long-time stability of the velocity–vorticity scheme with finite element spatial discretization and backward Euler temporal discretization. The algorithm decouples the vorticity equation by using a first order approximation of the vorticity in the momentum equation, and reads as follows.

Algorithm 3.1

Given the forcing f and initial velocity $u_0$, set $u_h^0$ to be the interpolant of $u_0$, and $w_h^0$ the interpolant of the curl of $u_0$. Select a timestep $\Delta t>0$, and for $n=0,1,2,\ldots$

Step 1 Find $(u_h^{n+1},P_h^{n+1})\in (X_h,Q_h)$ satisfy for every $(v_h,q_h)\in (X_h,Q_h)$,

$$\begin{aligned}&\frac{1}{\Delta t}( u_h^{n+1} - u_h^n ,v_h) + ( w_h^n \times u_h^{n+1},v_h) - (P_h^{n+1},\nabla \cdot v_h) \nonumber \\&\quad + \nu (\nabla u_h^{n+1},\nabla v_h) = (f^{n+1},v_h). \end{aligned}$$

(3.1)

$$\begin{aligned}&(\nabla \cdot u_h^{n+1},q_h) = 0, \end{aligned}$$

(3.2)

Step 2 Find $w_h^{n+1}\in Y_h$ satisfy for every $\chi _h\in Y_h$,

$$\begin{aligned} \frac{1}{\Delta t}( w_h^{n+1} - w_h^n,\chi _h) + b^*( u_h^{n+1} , w_h^{n+1},\chi _h) + \nu (\nabla w_h^{n+1},\nabla \chi _h) = (\nabla \times f^{n+1},\chi _h). \end{aligned}$$

(3.3)

We note that $P_h^{n}$ represents Bernoulli pressure, and thus is intended to approximate $\frac{1}{2} | u(t^{n}) |^2 + p(t^{n})$. To recover a zero-mean approximation to the kinematic pressure, one can rescale $p_h^{n}:=P_h^n - \frac{1}{2} |u_h^n|^2$ accordingly.

We will prove long-time $L^2$ and $H^1$ stability of both the velocity and the vorticity. We begin with the $L^2$ results.

Theorem 3.1

(Long-time $L^2$ stability of velocity and vorticity) Suppose $f\in L^{\infty }(0,\infty ;L^2(\Omega ))$, and $u_0 \in H^1(\Omega )$. Denote $\alpha :=(1+\nu \lambda ^{-2}\Delta t)$. For any $\Delta t>0$, we have that solutions of Algorithm 3.1 satisfy for every positive integer n,

$$\begin{aligned}&\Vert u_h^{n} \Vert ^2+\frac{\nu \Delta t}{2} \sum _{k=0}^{n-1}\left( \frac{1}{\alpha }\right) ^{n-k}\Vert \nabla u_h^{k+1}\Vert ^2 \le \left( \frac{1}{\alpha }\right) ^{n} \Vert u_h^0 \Vert ^2 + \frac{2\alpha \lambda ^{2}}{\nu ^2}\,\Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2 =: C_0^2, \end{aligned}$$

(3.4)

$$\begin{aligned}&\Vert w_h^{n} \Vert ^2+\frac{\nu \Delta t}{2} \sum _{k=0}^{n-1}\left( \frac{1}{\alpha }\right) ^{n-k}\Vert \nabla w_h^{k+1}\Vert ^2 \nonumber \\&\quad \le \left( \frac{1}{\alpha }\right) ^{n} \Vert w_h^0 \Vert ^2 + \frac{2\alpha \lambda ^{2}}{\nu ^2 }\,\Vert f \Vert _{L^{\infty }(0,\infty ;L^2(\Omega ))}^2 =: C_1^2, \end{aligned}$$

(3.5)

Remark 3.1

The constants $C_0$ and $C_1$ are independent of n and therefore hold for arbitrarily large n. These bounds can be considered as dependent only on the data (since time step sizes are inherently bounded above), and moreover, for sufficiently large n, the bounds are independent of the initial condition.

Proof

Take $v_h=2\Delta t u_h^{n+1}$, $q_h=P_h^{n+1}$, and $\chi _h=2 \Delta t w_h^{n+1}$, which vanishes the nonlinear and pressure terms, and leaves

$$\begin{aligned} \Vert u_h^{n+1} \Vert ^2 - \Vert u_h^n \Vert ^2 + \Vert u_h^{n+1} - u_h^n \Vert ^2 + 2\Delta t \nu \Vert \nabla u_h^{n+1} \Vert ^2&= 2\Delta t (f^{n+1},u_h^{n+1}), \\ \Vert w_h^{n+1} \Vert ^2 - \Vert w_h^n \Vert ^2 + \Vert w_h^{n+1} - w_h^n \Vert ^2 + 2\Delta t \nu \Vert \nabla w_h^{n+1} \Vert ^2&= 2\Delta t (\nabla \times f^{n+1},w_h^{n+1}). \end{aligned}$$

We majorize the forcing terms after integrating by parts in $(\nabla \times f^{n+1},w_h^{n+1}),$ applying Young’s inequality, and dropping positive terms on the left hand sides to get

$$\begin{aligned} \Vert u_h^{n+1} \Vert ^2 - \Vert u_h^n \Vert ^2 + \frac{3}{2}\nu \Delta t \Vert \nabla u_h^{n+1} \Vert ^2&\le 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert _{V_h^*}^2, \\ \Vert w_h^{n+1} \Vert ^2 - \Vert w_h^n \Vert ^2 + \frac{3}{2}\nu \Delta t \Vert \nabla w_h^{n+1} \Vert ^2&\le 2\nu ^{-1} \Delta t \Vert f^{n+1} \Vert ^2. \end{aligned}$$

From here, the velocity and vorticity estimates follow identically, except that the norm on the forcing term is different, and thus we restrict the remainder of the proof to only the velocity. Applying the Poincare inequality to lower bound the viscous term yields

$$\begin{aligned} (1+\nu \lambda ^{-2}\Delta t) \Vert u_h^{n+1} \Vert ^2 + \frac{\nu }{2}\Delta t \Vert \nabla u_h^{n+1} \Vert ^2&\le \Vert u_h^n \Vert ^2 + 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert _{V_h^*}^2. \end{aligned}$$

Now fix an integer $N>0$ and divide the above inequality by $\alpha ^{N-n}$ to obtain

$$\begin{aligned}&\left( \frac{1}{\alpha }\right) ^{N-n-1}\Vert u_h^{n+1} \Vert ^2 + \left( \frac{1}{\alpha }\right) ^{N-n} \frac{\nu }{2}\Delta t \Vert \nabla u_h^{n+1} \Vert ^2 \\&\quad \le \left( \frac{1}{\alpha }\right) ^{N-n}\Vert u_h^n \Vert ^2 + \left( \frac{1}{\alpha }\right) ^{N-n} 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert _{V_h^*}^2. \end{aligned}$$

Summing up for $n=0,\dots ,N-1$ and reducing, we get

$$\begin{aligned}&\Vert u_h^{N} \Vert ^2 +\frac{\nu \Delta t}{2} \sum _{n=0}^{N-1}\left( \frac{1}{\alpha }\right) ^{N-n}\Vert \nabla u_h^{n+1}\Vert ^2 \\&\quad \le \left( \frac{1}{\alpha }\right) ^{N} \Vert u_h^0 \Vert ^2 + 2\nu ^{-1}\Delta t \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2 \sum _{n=0}^{N-1} \left( \frac{1}{\alpha } \right) ^{N-n} \\&\quad \le \left( \frac{1}{\alpha }\right) ^{N} \Vert u_h^0 \Vert ^2 + 2\nu ^{-1}\Delta t \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2 \frac{ \alpha }{\alpha - 1 }. \end{aligned}$$

Substituting for $\alpha$ proves the velocity result. Applying the same steps for vorticity produces estimate (3.5), which finishes the proof of the theorem. $\square$

Theorem 3.2

(Long-time $H^1$ stability of velocity) Suppose $f\in L^{\infty }(0,\infty ;L^2(\Omega ))$, and $u_0 \in H^1(\Omega )$. Denote $\alpha :=(1+\nu \lambda ^{-2}\Delta t)$. For any $\Delta t>0$, the solutions of Algorithm 3.1 satisfy for every positive integer n,

$$\begin{aligned} \Vert \nabla u_h^{n} \Vert ^2 \le \left( \frac{1}{\alpha }\right) ^{n} \Vert \nabla u_h^0 \Vert ^2 + ( 2\nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2)}^2 + C \nu ^{-3}C_1^4 C_0^2) \frac{\alpha \lambda ^{2}}{\nu } =: C_2^2. \end{aligned}$$

(3.6)

and

$$\begin{aligned} \Vert \nabla u_h^{n} \Vert ^2 \le \left( \frac{1}{\alpha }\right) ^{n} \Vert \nabla u_h^0 \Vert ^2 + 2\frac{\alpha \lambda ^{2}}{\nu ^2 } \Vert f \Vert _{L^{\infty }(0,\infty ;L^2)}^2 + C(1+|\ln h|)\nu ^{-2} C_1^2 C_0^2 =: \widetilde{C}_2^2. \end{aligned}$$

(3.7)

where C is a generic constant, which depends on Sobolev’s embedding inequalities optimal constants and constants from Agmon’s type inequalities (2.1)–(2.4).

Remark 3.2

The theorem above proves that the long-time velocity solution is bounded in the $H^1$ norm only by the problem data, and similar to the $L^2$ bound, it is independent of the initial condition when n is sufficiently large.

With respect to the dependence on Re, the estimate (3.6) gives $\Vert \nabla u_h^{n} \Vert \le O(Re^5)$, while estimate (3.7) gives $\Vert \nabla u_h^{n} \Vert \le O((1+|\ln h|)^{\frac{1}{2}} Re^3)$.

Proof

Take $v_h=2\Delta t A_h u_h^{n+1}$ in (3.1) to obtain

$$\begin{aligned}&\Vert \nabla u_h^{n+1} \Vert ^2 - \Vert \nabla u_h^n \Vert ^2 + \frac{3}{2}\nu \Delta t \Vert A_h u_h^{n+1} \Vert ^2 \le 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert ^2 \\&\quad +\, 2\Delta t | ( w_h^n \times u_h^{n+1},A_h u_h^{n+1} ) | . \end{aligned}$$

For the last term on the right-hand side, we majorize it first using Holder’s inequality, the discrete Agmon inequality (2.1), Young’s inequality, and Theorem 3.1 to find

$$\begin{aligned} | ( w_h^n \times u_h^{n+1},A_h u_h^{n+1} ) |&\le \Vert w_h^n \Vert \Vert u_h^{n+1} \Vert _{L^{\infty }} \Vert A_h u_h^{n+1} \Vert \\&\le C \Vert w_h^n \Vert \Vert u_h^{n+1} \Vert ^{1/2} \Vert A_h u_h^{n+1} \Vert ^{3/2} \\&\le C\nu ^{-3} \Vert w_h^n \Vert ^4 \Vert u_h^{n+1} \Vert ^{2} + \frac{\nu }{2} \Vert A_h u_h^{n+1} \Vert ^{2} \\&\le C\nu ^{-3}C_1^4 C_0^{2} + \frac{\nu }{2} \Vert A_h u_h^{n+1} \Vert ^{2}. \end{aligned}$$

Combining these last two inequalities produces

$$\begin{aligned} \Vert \nabla u_h^{n+1} \Vert ^2 - \Vert \nabla u_h^n \Vert ^2 + \nu \Delta t \Vert A_h u_h^{n+1} \Vert ^2 \le 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert ^2 + C\Delta t \nu ^{-3} C_1^4 C_0^2, \end{aligned}$$

and thanks to Poincare, we obtain

$$\begin{aligned} (1+ \nu \lambda ^{-2}\Delta t )\Vert \nabla u_h^{n+1} \Vert ^2&\le \Vert \nabla u_h^n \Vert ^2 + \Delta t ( 2\nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2)}^2 + C \nu ^{-3} C_1^4 C_0^2 ). \end{aligned}$$

Recalling the notation $\alpha = (1+ \nu \lambda ^{-2}\Delta t )$, this relation can be written as

$$\begin{aligned} \Vert \nabla u_h^{n+1} \Vert ^2 \le \frac{1}{\alpha } \Vert \nabla u_h^n \Vert ^2 +\frac{1}{\alpha } \Delta t\,( 2\nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2)}^2 + C \nu ^{-3} C_1^4 C_0^2 ). \end{aligned}$$

(3.8)

Recursive substitution and an estimate for the partial sum of a geometric progression lead us to (3.6).

Alternatively, we can employ the finite element inverse inequality $\Vert {u}_h\Vert _{L^\infty (\Omega )}\le C(1+|\ln h|)^{\frac{1}{2}}\Vert \nabla {u}_h\Vert$, valid in 2D (see p. 124 in [3]), and estimate the nonlinear terms in the different way:

$$\begin{aligned} | ( w_h^n \times u_h^{n+1},A_h u_h^{n+1} ) |&\le \Vert w_h^n \Vert \Vert u_h^{n+1} \Vert _{L^{\infty }} \Vert A_h u_h^{n+1} \Vert \\&\le C (1+|\ln h|)^{\frac{1}{2}}\Vert w_h^n \Vert \Vert \nabla u_h^{n+1} \Vert \Vert A_h u_h^{n+1} \Vert \\&\le C(1+|\ln h|)\nu ^{-1} C_1^2 \Vert \nabla u_h^{n+1} \Vert ^{2} + \frac{\nu }{2} \Vert A_h u_h^{n+1} \Vert ^{2}. \end{aligned}$$

Similar arguments that produced (3.8) give

$$\begin{aligned} \Vert \nabla u_h^{n+1} \Vert ^2 \le \frac{1}{\alpha } \Vert \nabla u_h^n \Vert ^2 +\frac{1}{\alpha } \Delta t\,( 2\nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2)}^2 + C(1+|\ln h|) \nu ^{-1} C_1^2 \Vert \nabla u_h^{n+1} \Vert ^{2} ). \end{aligned}$$

Doing recursive substitution and employing (3.4) to estimate the resulting sum $\sum _{k=1}^{n+1}\alpha ^{k-n-1} \Vert \nabla u_h^{k} \Vert ^{2}$ leads to (3.7). $\square$

Theorem 3.3

(Long-time $H^1$ stability of vorticity) Suppose $f\in L^{\infty }(0,\infty ;H^1(\Omega ))$, and $u_0 \in H^2(\Omega )$. Let $\alpha :=(1+\nu \lambda ^{-2}\Delta t)$. For any $\Delta t>0$, solutions of Algorithm 3.1 satisfy for every positive integer n,

$$\begin{aligned} \Vert \nabla w_h^{n} \Vert ^2 \le \left( \frac{1}{\alpha }\right) ^{n} \Vert \nabla w_h^0 \Vert ^2 + ( 2\nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;H^1(\Omega ))} +\nu ^{-5} C_2^6 C_1^2 + C\nu ^{-3} C_2^4 C_1^2) \frac{\alpha \lambda ^{2}}{\nu }, \end{aligned}$$

(3.9)

and

$$\begin{aligned} \Vert \nabla w_h^{n} \Vert ^2 \le \left( \frac{1}{\alpha }\right) ^{n} \Vert \nabla w_h^0 \Vert ^2 + \frac{2\alpha \lambda ^{2}}{\nu ^2 } \Vert f \Vert _{L^{\infty }(0,\infty ;H^1(\Omega ))} + C (1+|\ln h|)\nu ^{-2} \widetilde{C}_2^2 C_1^2. \end{aligned}$$

(3.10)

Remark 3.3

The theorem above proves that the long-time vorticity solution is bounded in the $H^1$ norm only by the problem data, and similar to the $L^2$ bound, it is independent of the initial condition when n is sufficiently large.

Using $C_0=O(Re),\ C_1=O(Re)$ from Theorem 3.1, and $C_2=O(Re^5)$ from Theorem 3.2, the estimate (3.9) gives $\Vert \nabla w_h^{n} \Vert \le O(Re^{19})$. Estimate (3.10) gives $\Vert \nabla w_h^{n} \Vert \le O({ (1+|\ln h|)} Re^5)$, using additionally that $\widetilde{C}_2=O((1+|\ln h|)^{\frac{1}{2}} Re^3)$ from Theorem 3.2.

Proof

Take $\chi _h=2\Delta t \Delta _h w_h^{n+1}$ in (3.3), and majorize the forcing term using Cauchy–Schwarz and Young’s inequalities to obtain

$$\begin{aligned}&\Vert \nabla w_h^{n+1} \Vert ^2 - \Vert \nabla w_h^n \Vert ^2 + \frac{3}{2}\nu \Delta t \Vert \Delta _h w_h^{n+1} \Vert ^2 \le 2\nu ^{-1}\Delta t \Vert f^{n+1} \Vert _{H^1}^2 \\&\quad +\,2\Delta t | b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}) |. \end{aligned}$$

We bound the nonlinear term using Hölder inequality, Sobolev embeddings, discrete Agmon (2.2) and discrete Sobolev inequality (2.4), and Theorems 3.1 and 3.2 to reveal

$$\begin{aligned}&| b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}) | \\&\quad \le | (u_h^{n+1} \cdot \nabla w_h^{n+1},\Delta _h w_h^{n+1}) |+ \frac{1}{2} | ( (\nabla \cdot u_h^{n+1}) w_h^{n+1},\Delta _h w_h^{n+1}) | \\&\quad \le \Vert u_h^{n+1}\Vert _{L^6} \Vert \nabla w_h^{n+1} \Vert _{L^3} \Vert \Delta _h w_h^{n+1} \Vert + \frac{1}{2} \Vert \nabla u_h^{n+1} \Vert \Vert w_h^{n+1} \Vert _{L^{\infty }} \Vert \Delta _h w_h^{n+1} \Vert \\&\quad \le C C_2 \Vert w_h^{n+1} \Vert ^{1/3} \Vert \Delta _h w_h^{n+1} \Vert ^{5/3} + C C_2 \Vert w_h^{n+1} \Vert ^{1/2} \Vert \Delta _h w_h^{n+1} \Vert ^{3/2} \\&\quad \le C C_2 C_1^{1/3} \Vert \Delta _h w_h^{n+1} \Vert ^{5/3} + C C_2 C_1^{1/2} \Vert \Delta _h w_h^{n+1} \Vert ^{3/2}. \end{aligned}$$

The generalized Young’s inequality now provides the bound

$$\begin{aligned} | b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}) | \le C\nu ^{-5} C_2^6 C_1^2 + C\nu ^{-3} C_2^4 C_1^2 + \frac{\nu }{4} \Vert \Delta _h w_h^{n+1} \Vert ^2. \end{aligned}$$

Combining the estimates above yields

$$\begin{aligned}&\Vert \nabla w_h^{n+1} \Vert ^2 + \nu \Delta t \Vert \Delta _h w_h^{n+1} \Vert ^2 \le \Vert \nabla w_h^n \Vert ^2 + C\Delta t ( \nu ^{-1} \Vert f^{n+1} \Vert _{H^1} \\&\quad +\, \nu ^{-5} C_2^6 C_1^2 + C\nu ^{-3} C_2^4 C_1^2), \end{aligned}$$

and after applying Poincare we get

$$\begin{aligned} ( 1+ \lambda ^{-2}\nu \Delta t ) \Vert \nabla w_h^{n+1} \Vert ^2 \le \Vert \nabla w_h^n \Vert ^2 + C\Delta t ( \nu ^{-1} \Vert f^{n+1} \Vert _{H^1} + \nu ^{-5} C_2^6 C_1^2 + C\nu ^{-3} C_2^4 C_1^2). \end{aligned}$$

The remainder of the proof of (3.9) follows analogous to the $H^1$ case for velocity.

Alternatively, we may bound the nonlinear terms as follows:

$$\begin{aligned}&| b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}) | \\&\quad \le | (u_h^{n+1} \cdot \nabla w_h^{n+1},\Delta _h w_h^{n+1}) |+ \frac{1}{2} | ( (\nabla \cdot u_h^{n+1}) w_h^{n+1},\Delta _h w_h^{n+1}) | \\&\quad \le \Vert u_h^{n+1}\Vert _{L^\infty } \Vert \nabla w_h^{n+1} \Vert \Vert \Delta _h w_h^{n+1} \Vert + \frac{1}{2} \Vert \nabla u_h^{n+1} \Vert \Vert w_h^{n+1} \Vert _{L^{\infty }} \Vert \Delta _h w_h^{n+1} \Vert \\&\quad \le C(1+|\ln h|)^{\frac{1}{2}}\Vert \nabla u_h^{n+1}\Vert \Vert \nabla w_h^{n+1} \Vert \Vert \Delta _h w_h^{n+1} \Vert \\&\quad \le C(1+|\ln h|)^{\frac{1}{2}} \widetilde{C}_2 \Vert \nabla w_h^{n+1} \Vert \Vert \Delta _h w_h^{n+1} \Vert \\&\quad \le C \nu ^{-1} (1+|\ln h|) \widetilde{C}_2^2 \Vert \nabla w_h^{n+1} \Vert ^2 + \frac{\nu }{4} \Vert \Delta _h w_h^{n+1} \Vert ^2. \end{aligned}$$

To complete the proof of (3.10) we proceed as above and employ estimate (3.5) for the weighted sum of $\Vert \nabla w_h^{n+1} \Vert ^2$ norms. $\square$

4 Second-order method

We consider next a velocity–vorticity scheme with BDF2 timestepping. The scheme decouples the update of velocity and vorticity on each time step. Similar to the backward Euler case, we shall prove that the velocity and vorticity are both unconditionally long-time stable in both the $L^2$ and $H^1$ norms, and the scalings of the stability estimates with Re are the same as those from the backward Euler analysis. However, the analysis is somewhat more technical here, and a special norm is used to handle the time derivative terms.

Algorithm 4.1

Given the forcing f and initial velocity $u_0$, set $u_h^{-1}=u_h^0$ to be the interpolant of $u_0$, and $w_h^{-1}=w_h^0$ the interpolant of the curl of $u_0$. Select a timestep $\Delta t>0$, and for n=0,1,2,...

Step 1 Find $(u_h^{n+1},P_h^{n+1})\in (X_h,Q_h)$ satisfy for every $(v_h,q_h)\in (X_h,Q_h)$,

$$\begin{aligned}&\frac{1}{2\Delta t}( 3u_h^{n+1} - 4u_h^n + u_h^{n-1},v_h) + ( (2w_h^n - w_h^{n-1}) \times u_h^{n+1},v_h) \nonumber \\&\quad \quad - (P_h^{n+1},\nabla \cdot v_h) + \nu (\nabla u_h^{n+1},\nabla v_h) = (f^{n+1},v_h), \end{aligned}$$

(4.1)

$$\begin{aligned}&(\nabla \cdot u_h^{n+1},q_h) = 0. \end{aligned}$$

(4.2)

Step 2 Find $w_h^{n+1}\in Y_h$ satisfy for every $\chi _h\in Y_h$,

$$\begin{aligned}&\frac{1}{2\Delta t}( 3w_h^{n+1} - 4w_h^n + w_h^{n-1},v_h) + b^*( u_h^{n+1} , w_h^{n+1},\chi _h) \\&\quad + \nu (\nabla w_h^{n+1},\nabla \chi _h) = (\nabla \times f^{n+1},\chi _h). \end{aligned}$$

For the $2\times 2$ matrix

$$\begin{aligned} G := \left( \begin{array}{cc} 1/2 &{}\quad -1 \\ -1 &{}\quad 5/2 \end{array} \right) , \end{aligned}$$

we introduce the G-norm $\Vert \chi \Vert _{G} = (\chi ,G\chi )$, $\chi$ is vector valued. The G-norm is widely used in BDF2 analysis, see e.g. [4, 11]. The following property of the G-norm is well-known [11], however for completeness we will provide a short proof.

Lemma 4.1

Set $\chi ^0 = [v^0,\ v^1]^T$ and $\chi ^1 = [v^1,\ v^2]^T$. Then

$$\begin{aligned} \left( \frac{3}{2} v^2 -2v^1+ \frac{1}{2} v^0,v^2 \right) = \frac{1}{2} ( \Vert \chi ^1 \Vert _G^2 - \Vert \chi ^0 \Vert _G^2 ) + \frac{1}{4}\Vert v^2 - 2v^1 + v^0 \Vert ^2. \end{aligned}$$

Proof

Noting the algebraic identity $a(3a-4b+c)=\frac{1}{2}( (a^2-b^2) + (2a-b)^2-(2b-c)^2 + (a-2b+c)^2)$, we can write

$$\begin{aligned} \frac{1}{2}( 3v^2 - 4v^1 + v^0,v^2 )= & {} \frac{1}{2} \left( \frac{\Vert v^2 \Vert ^2 - \Vert v^1\Vert ^2}{2} + \frac{ \Vert 2v^2 - v^1 \Vert ^2 - \Vert 2v^1 - v^0\Vert ^2}{2} \right) \\&+\, \frac{1}{4} \Vert v^2 - 2v^1 + v^0 \Vert ^2 \\= & {} \frac{1}{2} \left( \frac{\Vert v^2 \Vert ^2 + \Vert 2v^2 - v^1 \Vert ^2}{2} - \frac{ \Vert 2v^1 - v^0\Vert ^2 + \Vert v^1\Vert ^2 }{2} \right) \\&+\, \frac{1}{4}\Vert v^2 - 2v^1 + v^0 \Vert ^2. \end{aligned}$$

The term $\Vert \chi ^0 \Vert _G^2$ can be decomposed as

$$\begin{aligned} \Vert \chi ^0 \Vert _G^2=(\chi ^0,G\chi ^0) = \left( \left( \begin{array}{c} v^0 \\ v^1 \end{array} \right) , \left( \begin{array}{c@{\quad }c} 1/2 &{} -1 \\ -1 &{} 5/2 \end{array} \right) \left( \begin{array}{c} v^0 \\ v^1 \end{array} \right) \right) , \end{aligned}$$

and with some arithmetic we get that

$$\begin{aligned} \Vert \chi ^0 \Vert _G^2 = \int _{\Omega } \left( v^0 \left( \frac{1}{2} v^0 - v^1\right) + v^1\left( -v^0 + \frac{5}{2} v^1\right) \right) \ dx = \frac{1}{2} \Vert v^1 \Vert ^2 + \frac{1}{2} \Vert 2v^1 - v^0 \Vert ^2. \end{aligned}$$

A similar decomposition of $\Vert \chi ^1\Vert$ completes the proof. $\square$

It is also known that the G norm is equivalent to the $L^2(\Omega )$ norm in the sense of there existing $C_l$ and $C_u$ such that

$$\begin{aligned} C_l \Vert \chi \Vert _G \le \Vert \chi \Vert \le C_u \Vert \chi \Vert _G. \end{aligned}$$

Use of the G-norm and this norm equivalence will allow for a smoother analysis.

We begin our analysis with the long-time $L^2$ stability of velocity and vorticity.

Theorem 4.1

(Long-time $L^2$ stability of velocity and vorticity) Let $f\in L^{\infty }(0,\infty ;V_h^{*})$ and $u_0\in H^1(\Omega )$. Then for any $\Delta t>0$, set $\alpha := \min \{ 1/2,\ \frac{\nu \Delta t C_l}{4} \}$, and it holds that solutions of Algorithm 4.1 satisfy for every positive integer n,

$$\begin{aligned}&C_l^2 ( \Vert u_h^{n} \Vert ^2 + \Vert u_h^{n-1} \Vert ^2 ) + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n} \Vert ^2 \le \left( \frac{1}{1+\alpha } \right) ^n \left( 2C_u \Vert u_h^0 \Vert ^2+ \frac{\nu \Delta t}{4}\Vert \nabla u_h^0 \Vert ^2 \right) \nonumber \\&\quad + \max \left( 2\Delta t, \frac{4}{\nu C_l^2} \right) \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2 =: C_4. \end{aligned}$$

(4.3)

If additionally $f\in L^{\infty }(0,\infty ;L^2(\Omega ))$ and $w_h^0 \in H^1(\Omega )$, then for any $\Delta t>0$, solutions of Algorithm 4.1 satisfy for every positive integer n,

$$\begin{aligned}&C_l^2 ( \Vert w_h^{n} \Vert ^2 + \Vert w_h^{n-1} \Vert ^2 ) + \frac{\nu \Delta t}{4} \Vert \nabla w_h^{n} \Vert ^2 \le \left( \frac{1}{1+\alpha } \right) ^n \left( 2C_u \Vert w_h^0 \Vert ^2 + \frac{\nu \Delta t}{4}\Vert \nabla w_h^0 \Vert ^2 \right) \nonumber \\&\quad + \max \left( 2\Delta t, \frac{4}{\nu C_l^2} \right) \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2(\Omega ))}^2 =: C_5. \end{aligned}$$

(4.4)

Remark 4.1

A more technical analysis can be made, similar to the backward Euler case, that includes the terms $\frac{\nu \Delta t}{2} \sum _{k=0}^{n-1} (1 + a)^{k-n} \Vert \nabla u_h^{k+1} \Vert ^2$ and $\frac{\nu \Delta t}{2} \sum _{k=0}^{n-1} \alpha ^{k-n} \Vert \nabla w_h^{k+1} \Vert ^2$ on the left hand sides of (4.3) and (4.4), respectively.

Proof

Choose $v_h=2\Delta t u_h^{n+1}$ in (4.1), which vanishes the nonlinear and pressure terms, and then upper bound the forcing term just as in the backward Euler case to get

$$\begin{aligned} \Vert \chi ^{n+1} \Vert _G^2 - \Vert \chi ^n \Vert _G^2 + \frac{1}{2} \Vert u_h^{n+1} - 2u_h^n + u_h^{n-1} \Vert ^2 + \nu \Delta t \Vert \nabla u_h^{n+1} \Vert ^2 \le \nu ^{-1} \Delta t \Vert f^{n+1} \Vert _{V_h^*}^2, \end{aligned}$$

(4.5)

where $\chi ^{n+1} = [u_h^{n},\ u_h^{n+1}]^T$ and $\chi ^{n} = [u_h^{n-1},\ u_h^{n}]^T$. Dropping the third term on the left-hand side, and adding $\frac{\nu \Delta t}{4}\Vert \nabla u_h^n \Vert ^2$ to both sides produces

$$\begin{aligned}&\left( \Vert \chi ^{n+1} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n+1} \Vert ^2 \right) + \frac{\nu \Delta t}{4}( \Vert \nabla u_h^{n+1} \Vert ^2 + \Vert \nabla u_h^n \Vert ^2 ) + \frac{\nu \Delta t}{2} \Vert \nabla u_h^{n+1} \Vert ^2 \nonumber \\&\quad \le \left( \Vert \chi ^n \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert \nabla u_h^n \Vert ^2 \right) + \nu ^{-1} \Delta t \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2. \end{aligned}$$

(4.6)

Using the Poincare inequality and then the equivalence of the G-norm with the $L^2$ norm, we have that

$$\begin{aligned}&\frac{\nu \Delta t}{4}( \Vert \nabla u_h^{n+1} \Vert ^2 + \Vert \nabla u_h^n \Vert ^2 ) \ge \frac{\nu \lambda ^{-2}}{\Delta t}{4}( \Vert u_h^{n+1} \Vert ^2 + \Vert u_h^n \Vert ^2 ) \\&\quad = \frac{\nu \Delta t}{4} \Vert \chi ^{n+1} \Vert ^2 \ge \frac{\nu \Delta t C_l^2}{4} \Vert \chi ^{n+1} \Vert _G^2, \end{aligned}$$

and thus setting $\alpha := \min \big \{ 1/2, \frac{\nu \Delta t C_l^2}{4} \big \}$, it holds that

$$\begin{aligned} \frac{\nu \Delta t}{4}( \Vert \nabla u_h^{n+1} \Vert ^2 + \Vert \nabla u_h^n \Vert ^2 ) + \frac{\nu \Delta t}{2} \Vert \nabla u_h^{n+1} \Vert ^2 \ge \alpha \left( \Vert \chi ^{n+1} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n+1} \Vert ^2 \right) . \end{aligned}$$

(4.7)

Combining (4.7) and (4.6) yields

$$\begin{aligned}&(1+ \alpha ) \left( \Vert \chi ^{n+1} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n+1} \Vert ^2 \right) \le \left( \Vert \chi ^n \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert \nabla u_h^n \Vert ^2 \right) \\&\quad +\, \nu ^{-1} \Delta t \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2, \end{aligned}$$

which immediately implies that

$$\begin{aligned}&\Vert \chi ^{n} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n} \Vert ^2 \le \left( \frac{1}{1+\alpha } \right) ^n \left( \Vert \chi ^0 \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert \nabla u_h^0 \Vert ^2 \right) \nonumber \\&\quad + \left( \frac{1}{1+\alpha } + \cdots +\left( \frac{1}{1+\alpha } \right) ^{n} \right) \nu ^{-1} \Delta t \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2. \end{aligned}$$

(4.8)

Since $\alpha >0$,

$$\begin{aligned} \frac{1}{1+\alpha } + \cdots + \left( \frac{1}{1+\alpha } \right) ^{n} = \frac{ \frac{1}{1+\alpha } - \left( \frac{1}{1+\alpha } \right) ^{n+1} }{1 - \frac{1}{1+\alpha } } \le \frac{ \frac{1}{1+\alpha } }{ \frac{\alpha }{1+\alpha } } = \frac{1}{\alpha } = \max \left\{ 2, \frac{4}{\nu \Delta t C_l^2} \right\} , \end{aligned}$$

and thus

$$\begin{aligned}&\Vert \chi ^{n} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert \nabla u_h^{n} \Vert ^2 \le \left( \frac{1}{1+\alpha } \right) ^n \left( \Vert \chi ^0 \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert \nabla u_h^0 \Vert ^2 \right) \nonumber \\&\quad + \max \left\{ 2\Delta t, \frac{4}{\nu C_l^2} \right\} \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;V_h^*)}^2. \end{aligned}$$

(4.9)

Now using the equivalence of norms for the G norm and $L^2$ norm of $\chi$ completes the velocity proof.

The proof for vorticity follows identically, modulo a higher order norm on the forcing, after taking the test function to be $w_h^{n+1}$. $\square$

We prove next the unconditional long-time $H^1$ stability of velocity.

Theorem 4.2

(Long-time $H^1$ stability of velocity) Let $f\in L^{\infty }(0,\infty ;L^2(\Omega ))$, $u_0\in H^2(\Omega )$, and set $\alpha := \min \{ 1/2, \frac{\nu \Delta t C_l^2}{4} \}$. Then for any $\Delta t>0$, solutions of Algorithm 4.1 satisfy for every positive integer n,

$$\begin{aligned}&C_l^2 ( \Vert \nabla u_h^{n} \Vert ^2 + \Vert \nabla u_h^{n-1} \Vert ^2 ) + \frac{\nu \Delta t}{4} \Vert A_h u_h^{n} \Vert ^2 \nonumber \\&\quad \le \left( \frac{1}{1+\alpha } \right) ^n \left( 2C_u \Vert \nabla u_h^0 \Vert ^2 + \frac{\nu \Delta t}{4}\Vert A_h u_h^0 \Vert ^2 \right) \nonumber \\&\quad \quad + \max \left( 2\Delta t, \frac{4}{\nu C_l^2} \right) ( \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2(\Omega ))}^2 + C\nu ^{-3} C_5^4 C_4^2 ) =: C_6.\nonumber \\ \end{aligned}$$

(4.10)

Remark 4.2

Similar to the backward Euler case, the long-time $H^1$ stability bound for velocity gives $\Vert \nabla u_h^{n} \Vert \le O(Re^5)$. If we instead bounded the nonlinear term as in the backward Euler case via

$$\begin{aligned} | 2\Delta t ( (2w_h^n - w_h^{n-1}) \times u_h^{n+1},A_h u_h^{n+1}) | \le C\Delta t(1+| \ln h |) \nu ^{-1} C_5^2 C_4^2 + \frac{\nu \Delta t}{2} \Vert A_h u_h^{n+1} \Vert ^2, \end{aligned}$$

then we can get instead $\Vert \nabla u_h^{n} \Vert \le O( (1+| \ln h |)^{1/2} Re^3 )$.

Proof

Choose $v_h=2\Delta t A_h u_h^{n+1}$ in (4.1), which vanishes the pressure terms, and then upper bound the forcing term to get

$$\begin{aligned}&( \Vert \chi ^{n+1} \Vert _G^2 - \Vert \chi ^n \Vert _G^2 ) + \frac{1}{2} \Vert \nabla ( u_h^{n+1} - 2u_h^n + u_h^{n-1} ) \Vert ^2 + \nu \Delta t \Vert A_h u_h^{n+1} \Vert ^2 \nonumber \\&\quad \le \nu ^{-1} \Delta t \Vert f^{n+1} \Vert ^2 - 2\Delta t ( (2w_h^n - w_h^{n-1}) \times u_h^{n+1},A_h u_h^{n+1}) , \end{aligned}$$

(4.11)

where $\chi ^{n+1} = [A_h^{1/2} u_h^{n},\ A_h^{1/2} u_h^{n+1}]^T$ and $\chi ^{n} = [A_h^{1/2} u_h^{n-1},\ A_h^{1/2} u_h^{n}]^T$. The last term on the right hand side is estimated using the same technique as in the backward Euler case from Sect. 3, and then applying the $L^2$ stability estimates (which is from Theorem 4.1 in this case):

$$\begin{aligned}&| 2\Delta t ( (2w_h^n - w_h^{n-1}) \times u_h^{n+1},A_h u_h^{n+1}) | \nonumber \\&\quad \le C \Delta t \nu ^{-3} ( \Vert w_h^{n} \Vert ^4 + \Vert w_h^{n-1} \Vert ^4 ) \Vert u_h^{n+1} \Vert ^2 + \frac{\nu \Delta t}{2} \Vert A_h u_h^{n+1} \Vert ^2 \nonumber \\&\quad \le C \Delta t \nu ^{-3} C_5^4 C_4^2 + \frac{\nu \Delta t}{2} \Vert A_h u_h^{n+1} \Vert ^2. \end{aligned}$$

(4.12)

Combining this with (4.11), dropping the second term on the left-hand side, and adding $\frac{\nu \Delta t}{4}\Vert A_h u_h^n \Vert ^2$ to both sides produces

$$\begin{aligned}&\left( \Vert \chi ^{n+1} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert A_h u_h^{n+1} \Vert ^2 \right) + \frac{\nu \Delta t}{4}( \Vert A_h u_h^{n+1} \Vert ^2 + \Vert A_h u_h^n \Vert ^2 ) + \frac{\nu \Delta t}{2} \Vert A_h u_h^{n+1} \Vert ^2 \nonumber \\&\quad \le \left( \Vert \chi ^n \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert A_h u_h^n \Vert ^2 \right) + \Delta t ( \nu ^{-1} \Vert f\Vert _{L^{\infty }(0,\infty ;L^2(\Omega ))}^2 + C\nu ^{-3} C_5^4 C_4^2 ).\nonumber \\ \end{aligned}$$

(4.13)

From here, setting $\alpha := \min \big \{ 1/2, \frac{\nu \Delta t C_l^2}{4} \big \}$ and taking analogous steps as in the proof of the long-time $L^2$ estimate (starting from (4.6)) provides us with

$$\begin{aligned}&\Vert \chi ^{n} \Vert _G^2 + \frac{\nu \Delta t}{4} \Vert A_h u_h^{n} \Vert ^2 \nonumber \\&\quad \le \left( \frac{1}{1+\alpha } \right) ^n \left( \Vert \chi ^0 \Vert _G^2 + \frac{\nu \Delta t}{4}\Vert A_h u_h^0 \Vert ^2 \right) + \max \left\{ 2\Delta t, \frac{4}{\nu C_l^2} \right\} \nonumber \\&\quad \quad \times ( \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;L^2(\Omega ))}^2+ C \nu ^{-3} C_5^4 C_4^2 ). \end{aligned}$$

(4.14)

Finally, applying the norm equivalence for the G-norm finishes the proof. $\square$

We can now prove the unconditional long-time $H^1$ stability of the vorticity.

Theorem 4.3

(Long-time $H^1$ stability of vorticity) Let $f\in L^{\infty }(0,\infty ;H^1(\Omega ))$, $u_0\in H^3(\Omega )$, and set $\alpha := \min \{ 1/2, \frac{\nu \Delta t C_l^2}{4} \}$. Then for any $\Delta t>0$, solutions of Algorithm 4.1 satisfy for every positive integer n,

$$\begin{aligned}&C_l^2 ( \Vert \nabla w_h^{n} \Vert ^2 + \Vert \nabla w_h^{n-1} \Vert ^2 ) + \frac{\nu \Delta t}{4} \Vert \Delta _h w_h^{n} \Vert ^2 \nonumber \\&\quad \le \left( \frac{1}{1+\alpha } \right) ^n \left( 2C_u \Vert \nabla w_h^0 \Vert ^2 + \frac{\nu \Delta t}{4}\Vert \Delta _h w_h^0 \Vert ^2 \right) \nonumber \\&\quad \quad + \max \left( 2\Delta t, \frac{4}{\nu C_l^2} \right) ( \nu ^{-1} \Vert f \Vert _{L^{\infty }(0,\infty ;H^1(\Omega ))}^2 + C \nu ^{-5} C_6^6 C_5^2 + \nu ^{-3}C_6^4 C_5^2 ) =: C_7.\nonumber \\ \end{aligned}$$

(4.15)

Remark 4.3

Similar to the backward Euler case, we find that $\Vert \nabla w_h^n \Vert \le O(Re^{19})$. However, different estimates of the nonlinear terms (i.e., using an inverse inequality as in the backward Euler case) can be used to find $\Vert \nabla w_h^n \Vert \le O((1+| \ln h |^{3/2} Re^{5})$.

Proof

Begin by choosing $\chi _h=-2\Delta t \Delta _h w_h^{n+1}$ to get

$$\begin{aligned}&( \Vert \chi ^{n+1} \Vert _G^2 - \Vert \chi ^n \Vert _G^2 ) + \frac{1}{2} \Vert \nabla ( w_h^{n+1} - 2w_h^n + w_h^{n-1} ) \Vert ^2 + \nu \Delta t \Vert \Delta _h w_h^{n+1} \Vert ^2 \nonumber \\&\quad \le \nu ^{-1} \Delta t \Vert \nabla \times f^{n+1} \Vert ^2 + 2\Delta t b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}), \end{aligned}$$

(4.16)

where $\chi ^{n+1} = [(-\Delta _h)^{1/2} w_h^{n},\ (-\Delta _h)^{1/2} w_h^{n+1}]^T$ and $\chi ^{n} = [(-\Delta _h)^{1/2} w_h^{n-1}, (-\Delta _h)^{1/2} w_h^{n}]^T$. Upper bounding the nonlinear term exactly as in the backward Euler case, and then using the long-time estimates proven above for the BDF2 scheme gives

$$\begin{aligned} | 2\Delta t b^*(u_h^{n+1},w_h^{n+1},\Delta _h w_h^{n+1}) |&\le C \nu ^{-5} C_6^6 C_5^2 + C\nu ^{-3}C_6^4 C_5^2 + \frac{\nu \Delta t}{4} \Vert \Delta _h w_h^{n+1} \Vert ^2, \end{aligned}$$

and thus using this and dropping the second term on the left side of (4.16) yields

$$\begin{aligned}&( \Vert \chi ^{n+1} \Vert _G^2 - \Vert \chi ^n \Vert _G^2 ) + \frac{\nu \Delta t}{2} \Vert \Delta _h w_h^{n+1} \Vert ^2 \le \nu ^{-1} \Delta t \Vert \nabla \times f^{n+1} \Vert ^2 \nonumber \\&\quad +\,C\Delta t( \nu ^{-5} C_6^6 C_5^2 + \nu ^{-3}C_6^4 C_5^2). \end{aligned}$$

(4.17)

From here, the same techniques as for the long-time $H^1$ stability of velocity can be used to complete the proof, modulo a higher norm on the forcing term. $\square$

5 Numerical experiments

We run several numerical experiments in order to test the long-time stability of Algorithm 4.1, which is the BDF2 timestepping algorithm for the proposed velocity–vorticity method. However, as our interest is in practical applications, we do not consider a test problem with periodic boundary conditions; instead, we consider 2D channel flow past a flat plate, which uses a Dirichlet velocity inflow, no-slip velocity on the walls, and a zero-traction outflow condition. Thus we must appropriately modify Algorithm 4.1 so that physical boundary conditions for the velocity and vorticity can be applied.

Fig. 1

Setup for the flow past a normal flat plate

As a numerical illustration of the long-term numerical stability, we compute the flow past a normal flat plate following [22, 23], see Fig. 1. We take as the domain $\Omega =[-7, 20] \times [-10, 10]$, with a hole of size $0.125 \times 1$ (representing the flat plate) removed from 7 units into the channel from the left, vertically centered. The inflow velocity is $u_{in}=\langle 1, 0 \rangle ^T$, and no-slip velocity is enforced on the walls and plate. Direct numerical simulations for this experiment are done for various Reynolds numbers Re, which can be considered here as $Re=\nu ^{-1}$, since the length of the plate is 1, and the inflow velocity has average magnitude 1. This is relatively simple, but interesting problem, which resembles the flow past other bluff objects. The flow undergoes a first Hopf bifurcation from steady to unsteady at a relatively low Reynolds numbers between 30 and 35 [22] and a second transition, also known as spatial transition from two-dimensional to three-dimensional, occurs around $Re=200$ [19]. We will test the velocity–vorticity algorithm and its long-time stability for $Re=100$ and $Re=125$.

The mathematical formulation of the problem has a constant in time non-homogeneous inflow boundary condition and zero source term. We deem this setting somewhat similar to the one analyzed in the paper [periodic boundary conditions and $L^{\infty }(0,\infty ;H^1(\Omega ))$-bounded right hand side], but more practically relevant.

5.1 Velocity–vorticity formulation with boundary conditions

Denote the domain by $\Omega$, with boundary $\partial \Omega =\Gamma _{in} \cup \Gamma _{out} \cup \Gamma _w$ split into inflow $\Gamma _{in}$, outflow $\Gamma _{out}$, and walls (of channel and plate) $\Gamma _w$. Denote by $\tau _h$ a regular, conforming triangulation of $\Omega$. The trial and test spaces for velocity functions are defined by

$$\begin{aligned} X_h^0&:= \{ v_h \in C^0(\Omega )^2 \cap P_2(\tau _h)^2,\ \ v_h |_{\Gamma _{in} \cup \Gamma _w}=0 \}, \\ X_h^g&:= \{ v_h \in C^0(\Omega )^2 \cap P_2(\tau _h)^2,\ v_h |_{\Gamma _{in} \cup \Gamma _w}=g \}, \end{aligned}$$

with $g=\langle 1, 0 \rangle ^T$ at the inflow, $g=\langle 0, 0 \rangle ^T$ on the walls, and with $P_2(\tau _h)$ denoting the space of globally continuous functions which are quadratic on each triangle. The discrete pressure space is taken to be

$$\begin{aligned} Q_h = \{ q \in C^0(\Omega ) \cap P_1(\tau _h) \}, \end{aligned}$$

and the zero traction boundary condition will be enforced weakly in the formulation. Note that $(X_h^0,Q_h)$ is the Taylor-Hood velocity–pressure element, which is known to be inf-sup stable [3]. The vorticity trial and test spaces are equal, since we take the vorticity at the inflow to be 0. The outflow condition for vorticity is a homogeneous Neumann condition, which is enforced weakly by the formulation. The appropriate vorticity boundary condition on $\Gamma _w$ is

$$\begin{aligned} (\nabla \times w) \times n=-(\nabla p)\times n \end{aligned}$$

where n in a normal vector on $\Gamma _w$, see [20]. In the finite element formulation this is a natural boundary condition, resulting in the presence of the following term, cf. [20],

$$\begin{aligned} \int _{\Gamma _w} p_h(\nabla \times {\chi }_h)\cdot n\,ds -\int _{\partial \Gamma _w} p_h\, {\chi }_h dl \quad \forall {\chi }_h\in W_h. \end{aligned}$$

The term is added to the formulation with the known pressure from Step 1. Thus the vorticity space is

$$\begin{aligned} W_h := \{ w_h \in C^0(\Omega ) \cap P_2(\tau _h), {w}_h |_{\Gamma _{in}}=0 \}. \end{aligned}$$

A second modification is made to the algorithm to avoid using the Bernoulli pressure, since there is an outflow boundary. Here, we use the identity from [2],

$$\begin{aligned} (\nabla \times u) \times u + \nabla \left( p + \frac{1}{2} |u|^2 \right) = \frac{1}{2} (\nabla \times u) \times u + \nabla p + D(u)u, \end{aligned}$$

where $D(u) = \frac{1}{2} \left( \nabla u + (\nabla u)^T \right)$ is the rate of deformation tensor. Thus the ‘do-nothing’ conditions we use on the outflow boundary correspond to enforcing

$$\begin{aligned} -\nu \frac{\partial u}{\partial n}+ p\,n=0,\quad \frac{\partial w}{\partial n}=0\quad \text {on}~\Gamma _{out} \end{aligned}$$

in the strong formulation.

Since there is no forcing in this test problem, we set $f=0$, and thus now Steps 1 and 2 of Algorithm 4.1 can now be written as they are computed:

Step 1 Find $({u}_h^{n+1},p_h^{n+1})\in (X_h^g,Q_h)$ satisfying

$$\begin{aligned}&\frac{1}{2\Delta t}(3u_h^{n+1} - 4u_h^n + u_h^{n-1},{v}_h) + \frac{1}{2} ((2{w}_h^{n}-{w}_h^{n-1}) \times {u}_h^{n+1},{v}_h)\quad \\&\quad + ( D(u_h^{n+1}) (2u_h^n - u_h^{n-1}),v_h ) - (p_h^{n+1},\nabla \cdot {v}_h) + \nu (\nabla {u}_h^{n+1},\nabla {v}_h) = 0\quad \forall {v}_h \in X_h^0 \\&(\nabla \cdot {u}_h^{n+1},q_h) = 0 \quad \forall q_h \in Q_h. \end{aligned}$$

Step 2 Find $w_h^{n+1} \in W_h$ satisfying

$$\begin{aligned}&\frac{1}{2\Delta t}(3{w}_h^{n+1} - 4{w}_h^n + w_h^{n-1},{\chi }_h) + ({u}_h^{n+1} \cdot \nabla {w}_h^{n+1},{\chi }_h) + \nu (\nabla {w}_h^{n+1},\nabla {\chi }_h) \nonumber \\&\quad = - \int _{\Gamma _w} p_h^{n+1}(\nabla \times {\chi }_h)\cdot n\,ds +\int _{\partial \Gamma _w} p_h^{n+1}\, {\chi }_h dl \quad \forall {\chi }_h\in W_h. \end{aligned}$$

(5.1)

We note that since globally continuous pressure elements are used, the right hand side of (5.1) can be equivalently written as

$$\begin{aligned} - \int _{\Gamma _w} p_h^{n+1}(\nabla \times {\chi }_h)\cdot n\,ds +\int _{\partial \Gamma _w} p_h^{n+1}\, {\chi }_h dl = \int _{\Gamma _w} (\nabla p_h^{n+1} \times n) \cdot {\chi }_h \,ds. \end{aligned}$$

5.2 Channel flow past a flat plate at $Re = 100$ and $Re = 125$

The BDF2 velocity–vorticity scheme was computed for both $Re=100$ and $Re=125$ ($\nu =Re^{-1}$), using three Delaunay generated triangular meshes which provided 79,509 total degrees of freedom (dof), 116,045 dof, and 159,055 dof with the $(P_2^2,P_1,P_2)$ velocity–pressure–vorticity elements. The simulations started the flow from rest $(u_h^0=0)$, and were run to an endtime $T=200$. For each mesh, several timestep choices were made, starting with $\Delta t=0.04$, and then cutting $\Delta t$ in half until convergence (i.e., successive solutions’ statistics matched). For both $Re=100$ and $Re=125$, the smallest $\Delta t$ was 0.01.

Quantities of interest for this problem is the long-time average of the drag coefficient $C_d$, and the Strouhal number. The Strouhal number was calculated as in [22, 23], using the fast Fourier transform of the transverse velocity at (4.0, 0.0) from $T=120$ to $T=200$. The drag coefficients are defined at each $t^n$ to be

$$\begin{aligned} C_d(t^m)&= \frac{2}{\rho L U_{max}^2} \int _{S} \left( \rho \nu \frac{\partial u_{t_S}(t^m)}{\partial n} n_y - p_h^m n_x \right) \ dS, \end{aligned}$$

where S is the plate, $\mathbf{n}=\langle n_x,n_y \rangle$ is the outward normal vector to S pointing into the domain, ${u}_{t_S}(t^m)$ is the tangential velocity of $u_h^m$, the density $\rho =1$, the max velocity at the inlet $U_{max}=1$, and $L=1$ is the length of the plate. The integral is calculated by transforming it into a global integral, which is believed to be more accurate [15]. The results for time averaged $C_d$ and the Strouhal numbers from the simulations for each Re, and for each mesh (with $\Delta t=0.01$), are shown in Table 1, along with reference values taken from [23]. We observe that the 116 K dof mesh and the 159 K dof meshes agree well with the reference values at $Re=100$ and $Re=125$. It appears we have achieved (or are close to) grid-convergence, and we note that for the Strouhal number, since the FFT was used with 8000 timesteps, 0.177 was the closest discrete frequency value to 0.174, and 0.189 was the next biggest discrete value compared to 0.183. We also plot the time-averaged vorticity in Fig. 2, and instantaneous velocity (as speed contours) in Fig. 3; both plots match the reference plots given in [23].

Table 1 Shown above are Strouhal numbers and long-time average drag coefficients for solutions on varying meshes, for $Re=100$ and $Re=125$

Also of interest is the stability of computed solutions in the $\Vert u_h^n \Vert _{L^2},\ \Vert u_h^n \Vert _{H^1}, \Vert w_h^n \Vert _{L^2},\ \Vert w_h^n \Vert _{H^1}$ norms versus time $t^n$, since we proved in Sect. 4 that these norms are all long-time stable (at least, in the periodic setting), independent of the timestep $\Delta t$ and mesh width h. Plots of these norms versus time are shown for $Re=100$ in Fig. 4 and for $Re=125$ in Fig. 5 for varying timesteps. Each norm appears to be long-time stable. Moreover, we do not observe the very large scaling of any of the norms with Re. Although $\Vert \nabla w_h^n\Vert \approx O(500)$ is an order of magnitude larger than $\Vert w_h^n\Vert$, it is still a very reasonable size and nowhere near $O(Re^{19})$ or even $O((1+|\ln h|)^{\frac{3}{2}} Re^{5})$.

Fig. 2

Shown above are plots of the time-averaged vorticity contours

Fig. 3

Shown above are plots of the speed contours of the velocity solutions at $T=200$

Fig. 4

Shown above are plots of the $Re=100$ solution norms versus time, found using Mesh 3 (the finest mesh)

Fig. 5

Shown above are plots of the $Re=125$ solution norms versus time, found using Mesh 3 (the finest mesh)

6 Conclusions and future directions

We have proven unconditional long-time stability of a scheme based on a velocity–vorticity formulation, and a finite-element-in-space BDF2-in-time IMEX discretization for the 2D Navier–Stokes equations. Long-time stability was proven in both the $L^2$ and $H^1$ norms for both velocity and vorticity, and the estimates hold for any $\Delta t>0$. The scheme is non-standard, and so we tested it on a benchmark problem on flow past a flat plate; it performed very well.

It would be interesting to study Algorithm 4.1, and variations thereof, for 3D flows. The difference in 3D is that the vortex stretching term $-(w\cdot \nabla u)$ appears in the vorticity equation. Since the 2D algorithm is proven herein to be unconditionally long-time stable, any instability in the 3D algorithm can be immediately attributed to the vortex stretching term and/or its numerical treatment. Isolating this behavior may give insight into better stabilization methods for higher Reynolds number flows in 3D.