A split-step finite-element method for incompressible Navier-Stokes equations with high-order accuracy up-to the boundary

TL;DR

This paper introduces a split-step finite-element method for incompressible Navier-Stokes equations that achieves high-order accuracy up to the boundary, improving efficiency and boundary-layer error handling.

Contribution

The paper presents a novel split-step finite-element algorithm with a new pressure boundary condition that enhances accuracy and flexibility for solving incompressible Navier-Stokes equations.

Findings

Achieves second-order accuracy for velocity and pressure up to the boundary.

Demonstrates robustness and efficiency with equal-order Lagrange finite elements.

Successfully applied to flow past a cylinder with higher-order elements.

Abstract

An efficient and accurate finite-element algorithm is described for the numerical solution of the incompressible Navier-Stokes (INS) equations. The new algorithm that solves the INS equations in a velocity-pressure reformulation is based on a split-step scheme in conjunction with the standard finite-element method. The split-step scheme employed for the temporal discretization of our algorithm completely separates the pressure updates from the solution of velocity variables. When the pressure equation is formed explicitly, the algorithm avoids solving a saddle-point problem; therefore, our algorithm has more flexibility in choosing finite-element spaces. For efficiency and robustness, Lagrange finite elements of equal order for both velocity and pressure are used. Motivated by a post-processing technique that calculates derivatives of a finite element solution with super-convergent error estimates, an alternative numerical boundary condition is proposed for the pressure equation at the discrete level. The new numerical pressure boundary condition that can be regarded as a better implementation of the compatibility boundary condition improves the boundary-layer errors of the pressure solution. Normal-mode analysis is performed using a simplified model problem on a uniform mesh to demonstrate the numerical properties of our methods. Convergence study using $\mathbb{P}_1$ elements confirms the analytical results and demonstrates that our algorithm with the new numerical boundary condition achieves the optimal second-order accuracy for both velocity and pressure up-to the boundary. Benchmark problems are also computed and carefully compared with existing studies. Finally, as an example to illustrate that our approach can be easily adapted for higher-order finite elements, we solve the classical flow-past-a-cylinder problem using $\mathbb{P}_n$ finite elements with $n\geq 1$.

Tables

Table 1. Table 1: Maximum values of the drag and lift coefficients and the pressure difference at final time t = 8 𝑡 8 t=8 .

	$C_{d,\max }$	$t\left(C_{d,\max }\right)$	$C_{l,\max }$	$t\left(C_{l,\max }\right)$	$\Delta p\left(8\right)$
TN (${\mathbb{P}}_1,{𝒢}_4$)	$2.9070$	$3.9343$	$0.4901$	$5.7189$	$-0.1115$
TN (${\mathbb{P}}_2,{𝒢}_2$)	$2.9527$	$3.9376$	$0.4871$	$5.6933$	$-0.1115$
TN (${\mathbb{P}}_4,{𝒢}_1$)	$2.9371$	$3.9365$	$0.4797$	$5.6885$	$-0.1111$
WABE (${\mathbb{P}}_1,{𝒢}_4$)	$2.9036$	$3.9341$	$0.4663$	$5.7274$	$-0.1111$
WABE (${\mathbb{P}}_2,{𝒢}_2$)	$2.9415$	$3.9339$	$0.4663$	$5.6989$	$-0.1126$
WABE (${\mathbb{P}}_4,{𝒢}_1$)	$2.9363$	$3.9365$	$0.4738$	$5.6890$	$-0.1103$
Reference Values
V. John [50]	$2.9509$	$3.9362$	$0.4779$	$5.6931$	$-0.1116$
Liu et al. [24]	$2.9541$	$3.9364$	$0.4791$	$5.6928$	$-0.1116$
Schäfer et al. [49]	$[2.930,2.970]$		$[0.470,0.490]$		$[-0.115,-0.105]$

Equations

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\mu \Delta \mathbf{u}+\mathbf{F}$$

$$\Delta p=-\rho \nabla \mathbf{u}:{\left(\nabla \mathbf{u}\right)}^T+\nabla \cdot \mathbf{F}+\alpha (\mathbf{x})\nabla \cdot \mathbf{u}$$

$$\mathcal{B}(\mathbf{u},p)=\mathbf{g}(\mathbf{x},t)$$

$$\nabla \cdot \mathbf{u}=0$$

$$\mathbf{u}(\mathbf{x},0)=\mathbf{f}(\mathbf{x})$$

$$\nabla \mathbf{u}:(\nabla \mathbf{u}{)}^T\equiv \sum _{i=1}^d\sum _{j=1}^d\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}.$$

$$\frac{\partial p}{\partial n}=\mathbf{n}\cdot \left(-\rho \frac{\partial \mathbf{u}}{\partial t}-\rho \mathbf{u}\cdot \nabla \mathbf{u}-\mu \nabla \times \nabla \times \mathbf{u}+\mathbf{F}\right),$$

$$\frac{\partial p}{\partial n}=\mathbf{n}\cdot \left(-\rho \frac{\partial \mathbf{u}}{\partial t}-\rho \mathbf{u}\cdot \nabla \mathbf{u}+\mu \Delta \mathbf{u}+\mathbf{F}\right),$$

$$\Delta \mathbf{u}=\nabla (\nabla \cdot \mathbf{u})-\nabla \times \nabla \times \mathbf{u}.$$

$$\mathbf{u}{|}_{\partial \Omega }=\mathbf{g}(\mathbf{x},t).$$

$$\rho \frac{{\mathbf{u}}^p-{\mathbf{u}}^n}{\Delta t}=\frac{3}{2}\left(L{\mathbf{u}}^n+{\mathbf{F}}^n\right)-\frac{1}{2}\left(L{\mathbf{u}}^{n-1}+{\mathbf{F}}^{n-1}\right)\text{with}L\mathbf{u}=-\rho \mathbf{u}\cdot \nabla \mathbf{u}-\nabla p+\mu \Delta \mathbf{u}.$$

$$\Delta p^p=-\rho \nabla {\mathbf{u}}^p:{\left(\nabla {\mathbf{u}}^p\right)}^T+\nabla \cdot {\mathbf{F}}^{n+1}+\alpha (\mathbf{x})\nabla \cdot {\mathbf{u}}^p.$$

$$\rho \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}=\frac{1}{2}\left(L{\mathbf{u}}^n+{\mathbf{F}}^n\right)+\frac{1}{2}\left(L{\mathbf{u}}^p+{\mathbf{F}}^{n+1}\right).$$

$$\Delta p^{n+1}=-\rho \nabla {\mathbf{u}}^{n+1}:{\left(\nabla {\mathbf{u}}^{n+1}\right)}^T+\nabla \cdot {\mathbf{F}}^{n+1}+\alpha (\mathbf{x})\nabla \cdot {\mathbf{u}}^{n+1}.$$

$${\mathbf{u}}_h^n=\sum _{k=1}^d\sum _{j=1}^N{u}_{k_j}^n{{\mathbf{\Phi }}_k}_j,p_h^n=\sum _{j=1}^N{p}_j^n{\phi }_j.$$

$$\forall f,g\in L_2(\Omega ):$$

$$\forall \mathbf{f},\mathbf{g}\in L_2(\Omega {)}^d:$$

$$\{\begin{array}{ll}\text{vskip}6.0ptplus2.0ptminus2.0pt{\mathbf{u}}_h^p({\mathbf{x}}_b)={\mathbf{g}}^{n+1}({\mathbf{x}}_b),& \forall {\mathbf{x}}_b\in \partial \Omega ,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\frac{\rho }{\Delta t}\left({\mathbf{u}}_h^p-{\mathbf{u}}_h^n,{\mathbf{v}}_h\right)=\frac{3}{2}\left[\left(L{\mathbf{u}}_h^n,{\mathbf{v}}_h\right)+({\mathbf{F}}^n,{\mathbf{v}}_h\right)]-\frac{1}{2}\left[(L{\mathbf{u}}_h^{n-1},{\mathbf{v}}_h)+({\mathbf{F}}^{n-1},{\mathbf{v}}_h)\right],& \forall {\mathbf{v}}_h\in V_{h0}^d.\end{array}$$

$$\{\begin{array}{l}\text{vskip}6.0ptplus2.0ptminus2.0pt-\left(\nabla p_h^p,\nabla q_h\right)+\lambda \left(1,q_h\right)=\left(-\rho \nabla {\mathbf{u}}_h^p:{\left(\nabla {\mathbf{u}}_h^p\right)}^T+\nabla \cdot {\mathbf{F}}^{n+1}+\alpha \nabla \cdot {\mathbf{u}}_h^p,q_h\right)-⟨\frac{\partial p_h^p}{\partial n},q_h⟩,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt(p_h^p,1)=0.\end{array}$$

$$⟨\frac{\partial p_h^p}{\partial n},q_h⟩=⟨\mathbf{n}\cdot \left(-\rho \frac{\partial {\mathbf{g}}^{n+1}}{\partial t}-\rho {\mathbf{g}}^{n+1}\cdot \nabla {\mathbf{u}}_h^p+{\mathbf{F}}^{n+1}\right),q_h⟩+\mu ⟨\nabla \times {\mathbf{u}}_h^p,\mathbf{n}\times \nabla q_h⟩,$$

$$⟨\mathbf{n}\cdot \nabla \times \nabla \times \mathbf{u},q⟩=\left(\nabla \times \nabla \times \mathbf{u},\nabla q\right)=-⟨\nabla \times \mathbf{u},\mathbf{n}\times \nabla q⟩.$$

$$\{\begin{array}{ll}\text{vskip}6.0ptplus2.0ptminus2.0pt{\mathbf{u}}_h^{n+1}({\mathbf{x}}_b)={\mathbf{g}}^{n+1}({\mathbf{x}}_b),& \forall {\mathbf{x}}_b\in \partial \Omega ,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\frac{\rho }{\Delta t}\left({\mathbf{u}}_h^{n+1}-{\mathbf{u}}_h^n,{\mathbf{v}}_h\right)=\frac{1}{2}\left[\left(L{\mathbf{u}}_h^n,{\mathbf{v}}_h\right)+({\mathbf{F}}^n,{\mathbf{v}}_h\right)]+\frac{1}{2}\left[(L{\mathbf{u}}_h^p,{\mathbf{v}}_h)+({\mathbf{F}}^{n+1},{\mathbf{v}}_h)\right],& \forall {\mathbf{v}}_h\in V_{h0}^d.\end{array}$$

$$\rho \left(\frac{\partial u_1}{\partial t}+\mathbf{u}\cdot \nabla u_1\right)=-\frac{\partial p}{\partial x_1}+\mu \left(\frac{{\partial }^2{u}_1}{\partial x_2^2}-\frac{{\partial }^2{u}_2}{\partial x_1\partial x_2}\right),$$

$$\rho \left(\frac{\partial u_2}{\partial t}+\mathbf{u}\cdot \nabla u_2\right)=-\frac{\partial p}{\partial x_2}+\mu \left(\frac{{\partial }^2{u}_2}{\partial x_1^2}-\frac{{\partial }^2{u}_1}{\partial x_1\partial x_2}\right).$$

$$\left(\frac{\partial p_h}{\partial x_1},{\phi }_{i_b}\right)$$

$$\left(\frac{\partial p_h}{\partial x_2},{\phi }_{i_b}\right)$$

$$\left({\mathbf{n}}^{i_b}\cdot \nabla p_h,{\phi }_{i_b}\right)=-\rho \left({\mathbf{n}}^{i_b}\cdot \left(\frac{\partial {\mathbf{u}}_h}{\partial t}+{\mathbf{u}}_h\cdot \nabla {\mathbf{u}}_h\right),{\phi }_{i_b}\right)+\mu \left(\nabla \times {\mathbf{u}}_h,{\mathbf{n}}^{i_b}\times \nabla {\phi }_{i_b}\right)+I_b,$$

$$\{\begin{array}{l}\text{vskip}6.0ptplus2.0ptminus2.0pt\text{For ∀i such that Pi∈Ω∂Ω and ∀ib such that Pib∈∂Ω},\\ \text{vskip}6.0ptplus2.0ptminus2.0pt-\left(\nabla p_h^{n+1},\nabla {\phi }_i\right)+\lambda \left(1,{\phi }_i\right)=\left(-\rho \nabla {\mathbf{u}}_h^{n+1}:{\left(\nabla {\mathbf{u}}_h^{n+1}\right)}^T+\alpha (\mathbf{x})\nabla \cdot {\mathbf{u}}_h^{n+1},{\phi }_i\right),\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\left({\mathbf{n}}^{i_b}\cdot \nabla p_h^{n+1},{\phi }_{i_b}\right)=-\rho \left({\mathbf{n}}^{i_b}\cdot \left(\frac{\partial {\mathbf{u}}_h^{n+1}}{\partial t}+{\mathbf{u}}_h^{n+1}\cdot \nabla {\mathbf{u}}_h^{n+1}\right),{\phi }_{i_b}\right)+\mu \left(\nabla \times {\mathbf{u}}_h^{n+1},{\mathbf{n}}^{i_b}\times \nabla {\phi }_{i_b}\right),\\ \text{vskip}6.0ptplus2.0ptminus2.0pt(p_h^{n+1},1)=0.\end{array}$$

$$\frac{\partial u}{\partial t}=-\frac{\partial p}{\partial x}+\nu \Delta u+f_u,$$

$$\frac{\partial v}{\partial t}=-\frac{\partial p}{\partial y}+\nu \Delta v+f_v,$$

$$\Delta p=\nabla \mathbf{f}+\alpha \nabla \cdot \mathbf{u},$$

$$u(x,0,t)=v(x,0,t)=0\text{and}$$

$$u(x,y,0)=v(x,y,0)=p(x,y,0)=0.$$

$$\{\begin{array}{l}\text{vskip}6.0ptplus2.0ptminus2.0pt\frac{\partial \hat{u}}{\partial t}=-ik\hat{p}-\nu k^2\hat{u}+\nu \frac{{\partial }^2\hat{u}}{\partial y^2}+{\hat{f}}_u,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\frac{\partial \hat{v}}{\partial t}=-\frac{\partial \hat{p}}{\partial y}-\nu k^2\hat{v}+\nu \frac{{\partial }^2\hat{v}}{\partial y^2}+{\hat{f}}_v,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt-k^2\hat{p}+\frac{{\partial }^2\hat{p}}{\partial y^2}=ik{\hat{f}}_u+\frac{\partial {\hat{f}}_v}{\partial y}+\alpha ik\hat{u}+\alpha \frac{\partial \hat{v}}{\partial y},\end{array}$$

$$\{\begin{array}{l}\text{vskip}6.0ptplus2.0ptminus2.0pt\hat{u}(k,0,t)=\hat{v}(k,0,t)=0,\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\frac{\partial \hat{p}}{\partial y}(k,0,t)=-\nu ik\frac{\partial \hat{u}}{\partial y}(k,0,t)+{\hat{f}}_v(k,0,t)\\ \text{vskip}6.0ptplus2.0ptminus2.0pt\hat{u}(k,y,0)=\hat{v}(k,y,0)=\hat{p}(k,y,0)=0.\end{array}$$

Full text

1. 1 Introduction

2. 2 Navier-Stokes equations in velocity-pressure form

3. 3 Temporal discretization

4. 4 Spatial discretization

5. 5 The complete numerical algorithm

6. 6 Weighted average over boundary elements (WABE) boundary condition

7. 7 A model problem for analysis

8. 7.1 Consistency and Error Equations

9. 7.2 Stability analysis

10. 7.2.1 Without divergence damping ($\alpha=0$)

11. 7.2.2 With divergence damping ($\alpha\neq 0$)

12. 7.3 Local Stability and Accuracy

13. 8 Numerical results

14. 8.1 Manufactured Solutions

15. 8.1.1 Periodic in $x$ direction

16. 8.1.2 No-slip boundary conditions on all boundaries

17. 8.2 Modified Lid-Driven Cavity

18. 8.3 Flow Past a Cylinder

19. 9 Conclusion

20. 10 Acknowledgment

21. A

1 Introduction

The development and analysis of numerical schemes for the fast and accurate solution of incompressible Navier-Stokes (INS) equations have long been a very active area of research, see for example [1, 2, 3, 4, 5, 6, 7] and the references therein. Popular existing numerical algorithms for solving INS equations include but are not limited to the following pioneering methods and their follow-up variants: (i) the MAC method that uses staggered grids for discretization [8]; (ii) projection methods [1] and their extension to an implicit fractional-step method [9]; (iii) the method of artificial compressibility [10]; (iv) split-step methods that solve an equivalent velocity-pressure reformulation of the INS equations [4, 6, 7]. There are also numerous other approaches based on common discretizations such as finite difference, finite element, finite volume, spectral element, and discontinuous Galerkin method; to name just a few, see [11, 12, 13, 14, 15, 16, 17].

The focus here is on solving the INS equations using finite element methods (FEM). The popular ones include methods based on the weak formulation of the INS equations in the primitive variables (also referred to as the velocity-divergence formulation in the literature), which often employ the $H^{1}(\Omega)^{d}\times L^{2}(\Omega)$ conforming elements for spatial discretization. Here $d=2$ or $3$ denotes the spatial dimension and $\Omega$ represents the fluid domain. Methods using the $H^{1}(\text{div};\Omega)\times L^{2}(\Omega)$ conforming finite elements of the Raviart-Thomas type are also proposed in [18] that better satisfy the divergence free condition of the fluid velocity. However, solving the INS equations in the primitive variables leads to a discrete saddle point problem; and, in order for this problem to have a solution, the choice of finite-element spaces for velocity and pressure must satisfy the well-known $\inf$-$\sup$ or Ladyzenskaja-Babuška-Brezzi (LBB) condition; these methods are often called mixed finite-element methods. For example, $\mathbb{Pe}_{2}$/$\mathbb{Pe}_{1}$ is a mixed velocity/pressure finite-element pair satisfying the LBB condition. The additional complexity posed by the LBB condition makes it hard to utilize these schemes for more complicated multi-physics applications such as Fluid-Structure interaction (FSI) problems. More discussion about popular finite element methods for solving the INS equations can be found in the classical book by Girault and Raviart [2].

Some projection methods based on the pioneering work of Chorin [1] avoid solving a saddle point problem by performing separate pressure updates through a splitting strategy. Therefore, the solution of these discrete problems is not subject to the LBB condition and this suggests more flexibility in choosing finite-element spaces for spatial discretization. However, the projection methods have their own drawback; that is, the pressure solution near boundaries is found to be less accurate with the presence of numerical boundary-layer errors [3, 5, 19]. Numerous techniques were developed attempting to reduce these numerical boundary layers in modern second order projection methods [9, 3, 20, 21]. Especially, Brown et al. [22] proposed a second order implicit projection method that is able to achieve full accuracy for both the velocity and pressure in rectangular domains; however, within a finite-element setting, it is not straightforward that their method is applicable for more general geometries due to the high order spatial derivatives required in their pressure update. More recently, works of Liu et al. have led to the development of a series of INS algorithms using either a velocity-pressure reformulation of the original INS equations or a variation of projection methods based on the Laplacian and Leray projection operators; these methods have been demonstrated to achieve promising results using standard finite-element discretization [7, 23, 24].

This paper concerns the development of a new FEM based INS algorithm that is suitable for FSI applications. A long-term motivator for the work is to extend the recently developed Added-Mass Partitioned (AMP) schemes for FSI simulations to the finite element framework. The AMP schemes were developed based on an interface condition derived at the continuous level by matching the time derivative of the kinematic interface condition. The AMP condition, which is a non-standard Robin-type boundary condition involving the fluid stress tensor, requires no adjustable parameters and in principle is applicable at the discrete level to couple the fluid and structure solvers of any accuracy and of any approximation methods (finite difference, finite element, finite volume, spectral element methods, etc). Within the finite-difference framework, the AMP algorithms have been developed and implemented to solve FSI problems involving the interaction of incompressible flows with a wide range of structures, such as elastic beams/shells [25, 26], bulk solids [27, 28, 29] and rigid bodies [30, 31, 32]. It has been shown in these works that the AMP schemes are second-order accurate and stable without sub-time-step iterations, even for very light structures when added-mass effects are strong. In contrast, the state-of-the-art finite-element based loosely-coupled partitioned FSI algorithms have yet achieved second-order accuracy for all the solution components of the fluid-structure system [33, 34, 35, 36, 37, 38, 39, 40]. Given the promising results achieved using the AMP schemes within the finite-difference framework, the goal of developing FEM based AMP schemes has the potential to improve the accuracy and efficiency of the state-of-the-art FSI simulations within the finite-element framework. However, the extension to FEM is nontrivial, and it poses specific requirements for the underlying INS solver. With this long-term goal in mind, the algorithm presented in this paper is designed to deal with a number of fundamental issues:

the scheme should be able to address the checker-board instability (or correspondingly the LBB stability condition in finite elements) so that the pressure solution is free of spurious oscillations;

• 2.

appropriate boundary conditions should be prescribed for an intermediate velocity field for projection type methods or for the pressure if using split-step type methods based on a velocity-pressure reformulation of the INS equations. These boundary conditions are essential to keep the boundary-layer errors small, and their choices are non-trivial;

• 3.

the scheme should be efficient and accurate, and it is often useful to decouple the solution of the velocity from the solution of the pressure;

• 4.

the scheme should be able to keep the discrete divergence small if it is not strictly enforced to be zero;

• 5.

the scheme should have the flexibility to work with non-standard boundary conditions.

Motivated by the success with the split-step finite difference algorithm [6], an FEM based split-step algorithm is developed for solving a pressure-velocity reformulation of the INS equations. The split-step strategy that separates the solution of pressure from that of the velocity variables enables more flexibility of choosing finite-element space for spatial discretization in the same manner as projection methods by avoiding a saddle point problem. Therefore, Lagrange (piecewise-polynomial) finite-elements of equal order for both velocity and pressure can be used by our method for efficiency. It is important to note that these standard elements can not be used in a straightforward way when solving the INS equations in the velocity-divergence formulation because they fail to satisfy the LBB condition. Therefore, our algorithm has the potential to be more efficient than many existing FEM based fluid solvers.

Special attention has been paid to investigate accurate numerical boundary conditions for the pressure equation. The curl-curl boundary condition appeared in the velocity-pressure reformulation is a compatibility boundary condition, which is derived by applying the normal component of the momentum equations on the boundary. The correct implementation of this compatibility boundary condition as a numerical boundary condition for the pressure is found to be crucial for the stability and accuracy of similar finite difference schemes [4, 41, 6]. To this end, two approaches of incorporating the curl-curl boundary condition within the finite-element context are considered that lead to two numerical boundary conditions of different accuracies. The most straightforward way to incorporate the curl-curl boundary condition, a Neumann boundary condition for the pressure, is to implement it as a natural boundary condition within the weak formulation of the pressure-Poisson equation. This approach, though simple and straightforward, is found to be less accurate since a slight degradation of the pressure accuracy near the boundary is observed in numerical experiments. Motivated by a post-processing technique that produces a super-convergent flux from finite element solutions [42], we propose an alternative way to implement the compatibility boundary condition as a more accurate numerical boundary condition for the pressure that alleviates the boundary-layer errors in the pressure solution.

The remainder of the paper is organized as follows. In Section 2, we discuss a pressure-velocity reformulation of the INS equations. A split-step scheme consisting of a second-order accurate predictor-corrector algorithm is described in Section 3 for the temporal discretization of the problem, and the discussion of the spatial discretization using the standard finite-element method follows in Section 4. The complete discrete algorithm is summarized in Section 5. In Section 6, novel numerical boundary conditions that help the pressure solution achieve better accuracy near the boundary are presented. The numerical properties of the scheme and the two pressure boundary conditions are analyzed for a simplified model problem discretized on a uniform mesh in Section 7. Careful numerical validations are conducted in Section 8. Finally, concluding remarks are made in Section 9.

2 Navier-Stokes equations in velocity-pressure form

An equivalent form, referred to as the “velocity-pressure” reformulation of the INS equations, is considered. Let $\Omega\subset\mathbb{R}^{d}$ ($d=2$ or $3$) denote a bounded open domain and $\partial\Omega$ be the boundary of $\Omega$. The velocity-pressure form are given by

[TABLE]

where

[TABLE]

Here $\mathbf{u}=(u_{1},\dots,u_{d})$ and $\mathbf{x}=(x_{1},\dots,x_{d})$ are the velocity and position in $d$-dimensional space; $p$ is the pressure; $\rho$ is the fluid density; $\mu$ is the coefficient of viscosity; and $\mathbf{F}(\mathbf{x},t)=\left(F_{1}(\mathbf{x},t),\dots,F_{d}(\mathbf{x},t)\right)$ is the external force. $\mathcal{B}(\mathbf{u},p)=\mathbf{g}(\mathbf{x},t)$ represents appropriate boundary conditions with $\mathbf{g}(\mathbf{x},t)=\left(g_{1}(\mathbf{x},t),\dotsm,g_{d}(\mathbf{x},t)\right)$ a given vector valued function. The given initial conditions, which are assumed to be divergence free, are represented by $\mathbf{f}(\mathbf{x})=\left(f_{1}(\mathbf{x}),\dots,f_{d}(\mathbf{x})\right)$.

Equation (2) is the pressure-Poisson equation (PPE), which is obtained by taking the divergence of the momentum equation together with the divergence free condition ($\nabla\cdot\mathbf{u}=0$). Note that a linear damping term $\alpha(\mathbf{x})\nabla\cdot\mathbf{u}$ is included in the PPE for numerical purposes. The damping term, referred to as the divergence damping, has no effect at the continuous level since $\nabla\cdot\mathbf{u}=0$; however, it helps to suppress the divergence at the discrete level, since the numerical solution is not exactly divergence-free due to properties of the numerical approximation. As is seen later in Section 5, our numerical algorithm does not enforce the divergence-free condition for the interior of the domain, instead this condition is only implicitly guaranteed by the curl-curl pressure boundary condition given below in equation (6). Discretization errors could result in the growth of the divergence of the velocity in the numerical solution. Therefore, it is important to include the divergence damping in the PPE to keep $\nabla\cdot\mathbf{u}$ small for the whole computation. Alternatively, one could perform an extra projection after every time step to map the velocity solutions into a divergence free space at the expense of a significant amount of additional computations.

The velocity-pressure formulation of the INS equations requires an extra boundary condition, and an appropriate choice is to set the divergence of the velocity to be zero on the boundary, or its normal derivative [6]. The split-step method considered in this paper separates the updates for velocity and pressure; as a consequence, a Poisson problem for the pressure needs to be solved explicitly at each time step. The divergence boundary condition given in (4) is, however, not convenient to be used as a boundary condition for the PPE (2) and hence an alternative condition is utilized. Following the studies in [12, 41, 6], the curl-curl boundary condition,

[TABLE]

has been used in place of (4) during the pressure update stage of our algorithm; noting that this curl-curl boundary condition is obtained by using the normal component of the momentum equations (1) as a compatibility boundary condition,

[TABLE]

and then replacing the diffusion term with the following vector identity,

[TABLE]

It is important to point out that the curl-curl boundary condition (6) has the divergence free condition implicitly implemented.

We remark that the INS equations in the velocity-pressure form were considered in [4, 43, 6], in which second- and forth-order accurate finite-difference based schemes have been developed and analyzed. In addition, the stability analysis of the curl-curl boundary condition in the context of a centered finite difference discretization is available in [41]. A similar velocity-pressure reformulation, without the divergence damping, was also investigated by Johnston and Liu [7]. It is interesting to note that, without the divergence damping, the INS equations in velocity-pressure form (5) are not equivalent to those in the primitive variables (also known as the velocity-divergence form) in the case of steady-state flows because solutions with constant $\nabla\cdot\mathbf{u}$ are possible for the velocity-pressure form; however, $\nabla\cdot\mathbf{u}$ is forced to be identically zero if the divergence damping is present. And we will see in the analysis and results sections that including the divergence damping is essential for our scheme to achieve optimal order of accuracy.

For simplicity, no-slip walls are assumed throughout this paper, in which case the velocity boundary conditions (3) are specifically given by

[TABLE]

3 Temporal discretization

Following [6], we use a split-step strategy to separate the updates of pressure and velocity components in the velocity-pressure reformulation (5). The split-step method is an explicit predictor-corrector method that consists of a second-order Adam-Bashforth (AB2) predictor and a modified second-order Adam-Moulton (AM2) corrector. We note that this AB2-AM2 time-stepping method has been successfully employed to solve (5) within a finite-difference framework in [6]. In this paper, we are interested in extending it to the finite-element framework. The spatial discretization using a standard finite element method will be discussed in section 4, and the summary of the full discrete algorithm follows in section 5.

To be specific, the velocity-pressure formulation (5) of the INS equations are advanced in time using the following time-stepping scheme. For simplicity, the algorithm is written for a fixed time-step, $\Delta t$, so that $t_{n}=n\Delta t$. We note that the algorithm can be extended to a variable $\Delta t$ with some reasonable strategies to dynamically determine the step size and the update frequency, but this case is not considered for the scope of this paper. Given solutions $(\mathbf{u}^{n-1},p^{n-1})$ and $(\mathbf{u}^{n},p^{n})$ at time levels $t_{n-1}$ and $t_{n}$, we first predict the velocity using the AB2 method,

[TABLE]

The pressure prediction is followed by solving the PPE with the predicted velocity solutions,

[TABLE]

The velocity is then corrected using the following modified AM2 method

[TABLE]

Note that the modified AM2 method is explicit since the predicted velocity $\mathbf{u}^{p}$ is used on the right hand side of the above equation. The pressure correction follows,

[TABLE]

We emphasize that this algorithm is stable without the use of the corrector step. Typically, the corrector step is included since the scheme has a larger stability region than the predictor alone, and the stability region includes the imaginary axis so that the scheme can be used for inviscid problems ($\mu=0$). The time step $\Delta t$ is determined by a diffusive stability constraint ($\Delta t\sim h^{2}$) for the explicit AB2-AM2 method, where $h$ is the grid spacing. Nevertheless, the time-step restriction can be alleviated if we treat the viscous term of the momentum equation implicitly using a Crank-Nicholson method; the time step for the semi-implicit scheme is determined by a convective stability constraint ($\Delta t\sim h$). One can refer to [41] for a proof of the time-step constraints for finite-difference based schemes.

4 Spatial discretization

It is important to note that the temporal discretization introduced in section 3 fully decouples the update of pressure from that of the velocity. When the pressure equation is formed explicitly, Lagrange finite elements of equal order can be used to discretize the velocity and pressure equations in a stable manner. This is in contrast to discretizing the velocity-divergence formulation in which case the standard Lagrange basis leads to an unstable scheme that does not satisfy the LBB condition.

In this paper, the standard Lagrange finite elements of equal order for both velocity and pressure are employed to discretize the above time-difference scheme in space; that is, we look for finite element solutions in the finite dimensional space $V_{h}^{d}\times Q_{h}=\mathbb{P}_{n}^{d}\times\mathbb{P}_{n}$. Here $\mathbb{P}_{n}$ is the piecewise polynomial finite-element space of degree $n$ that is defined on a triangulation of the domain $\mathcal{T}_{h}(\Omega)$. If $\{\varphi_{j}\}_{j=1}^{N}$ denote the basis functions of $\mathbb{P}_{n}$ with $N$ being the number of degrees of freedom, then a set of basis functions for $\mathbb{P}_{n}^{d}$ can be conveniently formed as ${\mathbf{\Phi}_{k}}_{j}=\varphi_{j}\mathbf{e}_{k}$ ( $k=1,\dots,d$ and $j=1,\dots,N$), where $\mathbf{e}_{i}$’s are the standard bases of $\mathbb{R}^{d}$. Thus, the finite-element approximations to velocity and pressure solutions at time $t_{n}$ can be represented as

[TABLE]

5 The complete numerical algorithm

To simplify discussion, the following notations are introduced for inner products defined over the domain $\Omega$ and its boundary $\partial\Omega$,

[TABLE]

And we denote $V_{h\mathbf{0}}^{d}$ the subspace of $V_{h}^{d}$ that vanishes on the boundary.

Given the finite element solutions $(\mathbf{u}^{n}_{h},p_{h}^{n})\in V_{h}^{d}\times Q_{h}$ at the current time $t_{n}$, and $(\mathbf{u}^{n-1}_{h},p_{h}^{n-1})\in V_{h}^{d}\times Q_{h}$ at one previous time level $t_{n-1}$, the goal of the algorithm is to determine the solution at time $t_{n+1}$. The complete discrete scheme using the above predictor-corrector time stepping method and the FEM spatial discretization is as follows.

Begin predictor.

Stage I - velocity prediction: We predict the velocity solution by solving for $\mathbf{u}^{p}_{h}\in V_{h}^{d}$ such that

[TABLE]

Here $(L\mathbf{u}_{h},\mathbf{v}_{h})=-\rho(\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{h},\mathbf{v}_{h})-(\nabla p_{h},\mathbf{v}_{h})-\mu(\nabla\mathbf{u}_{h},\nabla\mathbf{v}_{h})$. The superscripts over the given functions $\mathbf{F}$ and $\mathbf{g}$ indicate evaluating the functions at the corresponding time level.

Stage II - pressure update: It is important to point out that, with the no-slip boundary condition, the pressure is only determined up-to an additive constant in the PPE (other boundary conditions may remove this singularity). For those boundary conditions that imply a singular Poisson problem, we add an additional constraint $\left(p_{h},1\right)=0$ that sets the mean value of pressure to zero to the PPE as a Lagrange multiplier [6]. Specifically, we solve for $p^{p}_{h}\in Q_{h}$ and $\lambda\in\mathbb{R}$ such that, for $\forall\,q_{h}\in Q_{h}$,

[TABLE]

In practice, we chose $\alpha$ to be inversely proportional to the square of the mesh spacing. For a nonuniform mesh, the minimum of the spacings, $h_{\min}$, is used; that is, $\alpha=C_{d}h_{\min}^{-2}.$ In addition, the boundary integral is given by

[TABLE]

which is derived by utilizing the curl-curl boundary condition (6) and the following vector identity,

[TABLE]

We note that the vector identity (11) plays an important role in the algorithm since it reduces the regularity requirement for the admissible finite-element space that makes it possible to use $\mathbb{P}_{1}$ finite elements.

Remark: we refer to the pressure boundary condition (10) as the traditional Neumann (TN) boundary condition, which arises naturally from testing the PPE and integration by parts. However, a slight degradation of the pressure accuracy near the boundary is observed in numerical experiments. We think this degradation stems from the fact that we need to evaluate $\left<\nabla\times\mathbf{u}_{h},\mathbf{n}\times\nabla q_{h}\right>$ for the TN pressure boundary condition, and direct evaluation of $\nabla\times\mathbf{u}_{h}$ from the finite element solution for the velocity field $\mathbf{u}_{h}$ is of sub-optimal order of accuracy. Motivated by Carey’s post-processing technique that produces a super-convergent flux from finite element solutions [42], we propose an alternative compatibility boundary condition for the pressure at the discrete level. This new condition referred to as WABE boundary condition improves the pressure accuracy near the boundary, and its detail will be discussed in Section 6.

Begin corrector.

Stage III - velocity correction: To correct the velocity, we solve for $\mathbf{u}^{n+1}_{h}\in V_{h}^{d}$ such that

[TABLE]

Stage IV - pressure update: Finally, we obtain the updated pressure solution $p^{n+1}_{h}\in Q_{h}$ by solving the same equation (9) and (10) but with the corrected velocity solution $\mathbf{u}_{h}^{n+1}$ supplied on the right hand side.

6 Weighted average over boundary elements (WABE) boundary condition

In this section, an alternative numerical boundary condition for the PPE is introduced, which is essential for our algorithm to achieve high-order accuracy up-to the boundary. We derive this numerical boundary condition from a weighted average of the momentum equations over boundary elements; and therefore, we refer to it as WABE boundary condition for short. The WABE boundary condition avoids direct evaluation of $\nabla\times\mathbf{u}_{h}$ in leading order terms, and thus prevents degradation of the pressure accuracy from the boundary. Unless otherwise noted, the discussion is restricted to 2D in this paper.

Using the continuity equation in the Laplacian term of the momentum equations implies

[TABLE]

Not that the external forcing $\mathbf{F}$ is omitted here to save space, which can be easily included in the WABE boundary condition derived below. On the finite dimensional space $V_{h}^{d}\times Q_{h}$, testing the momentum equations (13) and (14) with the basis function $\varphi_{i_{b}}$ and integrating by parts, we have

[TABLE]

Here $i_{b}$ denotes the index of any boundary degree of freedom, i.e., $\forall P_{i_{b}}\in\partial\Omega$. Since $\varphi_{i_{b}}$ is nonzero only on the elements that contain the node $P_{i_{b}}$, the equations (15) and (16) are essentially a weighted average of the momentum equations over those boundary elements. Let $\mathbf{n}^{i_{b}}=\left(n_{1}^{i_{b}},n_{2}^{i_{b}}\right)$ be the unit outward normal vector at the boundary node $P_{i_{b}}$. To mimic the curl-curl boundary condition (6), the averaged momentum equations (15) and (16) are combined in $\mathbf{n}^{i_{b}}$ direction, and the WABE boundary condition is obtained:

[TABLE]

where $I_{b}=\mu\left<\left(\mathbf{n}\times\mathbf{n}^{i_{b}}\right)\cdot\left(\nabla\times\mathbf{u}_{h}\right),\varphi_{i_{b}}\right>$ accounts for the boundary integral resulted from the integration by parts. Here $\mathbf{n}$ is the normal on the edge and $\mathbf{n}^{i_{b}}$ is the normal on the boundary node. The less accurate boundary integral $I_{b}$ vanishes if the boundary of the domain is a straight line since $\mathbf{n}\times\mathbf{n}^{i_{b}}=0$; however, if the boundary is not a straight line but a smooth curve, we have $\mathbf{n}\times\mathbf{n}^{i_{b}}=\mathcal{O}\left(h^{2}\right)$ and thus the error introduced by $\nabla\times\mathbf{u}_{h}$ will be scaled down by an order of $h^{2}$. In either case, the WABE boundary condition (17) should be more accurate than the discrete TN boundary condition (10), and thus improves the boundary-layer errors that appear in the pressure solution. We note that $I_{b}$ can be ignored for second order accurate methods in practice.

To solve for pressure with the WABE boundary condition, we only need to replace the discrete pressure equations on the boundary nodes with the equations given by the WABE boundary condition (17). To be specific, we solve the following modified discrete PPE at the stages of pressure update as described in the previous section,

[TABLE]

It is important to remark that the implementation of the WABE boundary condition is not computationally more expensive than the TN boundary condition since there is no additional calculation needed to facilitate the WABE condition, and all the data needed, except for ${\partial\mathbf{u}^{n+1}_{h}}/{\partial t}$, are directly available from the previous velocity updates. However, a sufficiently accurate value for the purpose of the boundary condition can be obtained using a forward finite difference formula in time, i.e., ${\partial\mathbf{u}^{n+1}_{h}}/{\partial t}=(\mathbf{u}_{h}^{n+1}-\mathbf{u}_{h}^{n})/\Delta t$.

7 A model problem for analysis

In this section, we perform a normal-mode analysis on a model problem to reveal the numerical properties of the algorithm. The motivation that we analyze the model problem for a particular geometric domain using finite difference theory as opposed to using a more traditional approximation theory and energy estimation is as follows. As we know, the standard energy estimate method concerns the accuracy results in $L_{2}$ norm, which averages the numerical errors over the entire domain. However, one of the main focuses of this paper is to address the boundary-layer errors for the pressure solution; therefore, it is more appropriate to analyze numerical errors in $L_{\infty}$ norm. Using a normal-mode analysis for a particular geometric domain, we essentially solve the model problem analytically. Even though this analysis does not apply to general geometries and unstructured meshes, with the analytical solution of a model problem, we are able to identify some subtleties that are buried inside the $L_{2}$ averaging process, i.e., the boundary-layer errors in the pressure solution. The idea of rewriting finite element schemes on a uniform mesh as finite difference ones for analytical purposes was also employed by other researchers; for instance, in [47], the authors rewrote their discontinuous Galerkin schemes as finite difference ones, and then performed a Fourier type analysis since it was not easy for them to use standard finite element techniques to prove the inconsistency and weak instability of their schemes. Utilizing a mode analysis, they were able to analytically identify the weak instability that was observed in their numerical experiments.

Here Stokes equations are used as a model problem for the INS equations since the nonlinear convection terms can be regarded as lower order terms that do not affect the stability of the scheme. To further simplify the discussion, the model problem is assumed to be $2\pi$-periodic in $x$ direction on a semi-infinite domain $[0,2\pi]\times[0,\infty]$. Although the model problem drastically simplifies the original INS equations, the analysis performed here can shed some light on the stability and accuracy of our proposed method. Specifically, we consider the following initial-boundary value problem (IBVP) with $\rho=1$ and $\nu=\mu/\rho$,

[TABLE]

We impose the no-slip and curl-curl pressure boundary conditions at the boundary $y=0$,

[TABLE]

The homogeneous initial conditions are imposed to complete the statement of the model problem,

[TABLE]

Utilizing the assumption of periodicity, the IBVP can be Fourier transformed in $x$ direction; thus, in the Fourier space, we have

[TABLE]

subject to the transformed boundary and initial conditions,

[TABLE]

Here $\hat{u}(k,y,t)$, $\hat{v}(k,y,t)$, and $\hat{p}(k,y,t)$ are Fourier transformations of $u(x,y,t)$, $v(x,y,t)$ and $p(x,y,t)$ with wave number $k\in\mathbb{Z}$. For regularity, we also require

[TABLE]

where the norm is the standard $L_{2}$ function norm.

To analyze our numerical scheme, we discretize the above transformed equations using $\mathbb{Pe}_{1}$ finite elements on a uniform Cartesian grid, $\mathcal{G}=\{y_{j}=jh~{}|~{}\forall j\in\mathbb{N}\}$. Here $j=0$ corresponds to the boundary node. Thus, the finite element approximations to $\hat{u}(k,y,t)$, $\hat{v}(k,y,t)$, and $\hat{p}(k,y,t)$ can be represented as

[TABLE]

with $\varphi_{j}(y)$ denoting the $j$th basis function for $\mathbb{Pe}_{1}$ that is defined by

[TABLE]

For any grid function $f_{j}$, we introduce the following difference operators

[TABLE]

where $M$ is an average operator, and $D_{+}$, $D_{-}$ and $D_{0}$ are the forward, backward and centered divided difference operators, respectively. With these difference operators, it is easily seen that

[TABLE]

Similar relations can be found for $\hat{v}_{h}$ and $\hat{p}_{h}$ as well. Therefore, the finite element discretization of the transformed equations can be regarded as a finite difference scheme, and thus can be analyzed using finite difference techniques.

Specifically, we rewrite our finite element scheme as the following finite difference scheme,

[TABLE]

subject to the no-slip boundary conditions $u_{0}=v_{0}=0$, and either one of the discrete TN (23) and WABE (24) boundary conditions that are described below. Here $\alpha=C_{d}/h^{2}$ is assumed where $C_{d}$ is a constant independent of $h$.

To derive the discrete TN boundary condition for the model problem, we test the pressure equation in (19) with $\varphi_{0}$; that is

[TABLE]

With the regularity condition (21) and the curl-curl condition given in (20), we have

[TABLE]

thus the TN boundary condition can be written as

[TABLE]

On the other hand, the WABE condition for the transformed model problem is given by

[TABLE]

which can be written as the following finite difference form,

[TABLE]

7.1 Consistency and Error Equations

Assuming that the continuous problem given in (19) and (20) has a smooth solution, $(\hat{u},\hat{v},\hat{p})$, we introduce it into the difference equations (22) and obtain that the order of the truncation error for the discretization is $\mathcal{O}\left(h^{2}\right)$. In addition, introducing the smooth solution into the boundary conditions (23) and (24), we obtain the orders of the truncation errors for the TN and WABE boundary conditions are $\mathcal{O}\left(h\right)$ and $\mathcal{O}\left(h^{2}\right)$, respectively. Furthermore, we can readily write down the equations for the errors $U_{i}=u_{i}-\hat{u}_{i}$, $V_{i}=v_{i}-\hat{v}_{i}$ and $P_{i}=p_{i}-\hat{p}_{i}$, where $\hat{u}_{j}=\hat{u}(k,y_{i},t)$, $\hat{v}_{j}=\hat{v}(k,y_{i},t)$ and $\hat{p}_{j}=\hat{p}(k,y_{i},t)$ are the exact solutions evaluated at $y_{j}$ for wavenumber $k$. The consistency and error equations of the discretized problem are summarized in the following proposition; more detail of the analysis can be found in A

Proposition 1**.**

The numerical scheme (22) subject to either the TN (23) or the WABE (24) boundary conditions is consistent. The discretization is formally $\mathcal{O}\left(h^{2}\right)$ for the interior equations, $\mathcal{O}\left(h\right)$ for the TN boundary condition and $\mathcal{O}\left(h^{2}\right)$ for the WABE boundary condition. The error equations for the discretized problem are given by

[TABLE]

where $F_{u}$, $F_{v}$ and $F_{p}$ are some functions of $\mathcal{O}\left(1\right)$, and the boundary errors are given by

[TABLE]

where $r=1$ for TN condition and $r=2$ for WABE condition with some function $g_{0}=\mathcal{O}\left(1\right)$.

As a remark, we look at the order of the truncation error of the divergence damping term to see why the choice of a large coefficient, $\alpha=C_{d}/h^{2}$, does not affect the order of the truncation error of the whole scheme. The damping term appears in the pressure equation as well as in the TN boundary condition. In the pressure equation, the expansion of the damping term leads up to

[TABLE]

The continuity condition implies

[TABLE]

So the divergence damping term with the choice of $\alpha=C_{d}/h^{2}$ contributes an $\mathcal{O}\left(h^{2}\right)$ error that is in line with the accuracy of the other difference operators in the scheme. Similarly, for the divergence damping term in the TN boundary condition, we have

[TABLE]

Thus, with the continuity condition, we see that the error contributed by the divergence damping is $\mathcal{O}\left(h\right)$ that is also consistent with the accuracy of the discrete TN boundary condition.

7.2 Stability analysis

It suffices to consider the homogeneous version of the problem (25) & (26) when analyzing the stability of the scheme. For simplicity, we do not discretize in time and analyze the stability properties of the semi-discrete problem directly using Laplace transformation method and normal-mode analysis. It is mentioned in [43, 44, 45] that any dissipative time discretization can be used and the resulting fully discrete problem will be stable provided the semi-discrete problem is stable.

As is pointed out in [45], there are several possible stability definitions. Here we show the semi-discrete problem is stable in the sense of Godunov-Ryabenkii condition; that is, we demonstrate its stability by showing that there is no eigenvalue $s$ with $\Re(s)>0$ for a related eigenvalue problem, which is obtained by Laplace transforming the semi-discrete problem (25) in time with $s$ denoting the dual variable.

To be specific, after Laplace transforming the homogeneous version of (25) & (26), we obtain the eigenvalue problem,

[TABLE]

with the boundary conditions,

[TABLE]

and the regularity condition,

[TABLE]

Here ($\tilde{U}_{j}$, $\tilde{V}_{j}$, $\tilde{P}_{j}$) denote the transformed solutions and $||\cdot||_{h}$ represents the discrete $L_{2}$ norm defined on the grid $\mathcal{G}$.

Note that if $k=0$, we have a non-trivial solution to the eigenvalue problem (27) & (28) , i.e.,

[TABLE]

This solution should be excluded due to the regularity condition $||\tilde{P}||_{h}<\infty$. This case corresponds to the undetermined constant in the pressure which is regularized in our algorithm by enforcing the mean of the pressure to be zero as is described in (9). So we proceed the stability analysis assuming $k\neq 0$.

7.2.1 Without divergence damping ($\alpha=0$)

First, let us consider $\alpha=0$. In this case, the pressure is decoupled from the velocity equations, and the solution of the eigenvalue problem can be found explicitly, which is described in the following lemma.

Lemma 1**.**

If $\alpha=0$, the solution of the eigenvalue problem (27) & (28) is found to be

[TABLE]

where $\xi$ and $\gamma$ satisfy

[TABLE]

The remaining coefficient $C_{p}$ will be determined by the pressure boundary condition in (28), which implies

[TABLE]

If we let

[TABLE]

then (30) can be written as $C_{p}q(s)=0$. For $\Re(s)>0,$ we have the following lemmas that imply $q(s)\neq 0$; therefore, we conclude that $C_{p}=0$ and the solution given in (29) is trivial.

Lemma 2**.**

If $q_{1}(s)$ is real, then $s$ is real.

Lemma 3**.**

$q(s)\neq 0$* for $\Re(s)>0$*

For a concise presentation of the main results, the more technical proofs for these Lemmas are shown later in A, and here some examples of $q_{1}(s)$ and $q(s)$ for $s>0$ and various wavenumbers are plotted in Figure 1 to help demonstrate the analytical results. From the plots, we see that both $q_{1}(s)$ and $q(s)$ are decreasing and no root of $q(s)$ is observed for $s>0$, which are consistent with the analysis.

In conclusion, the eigenvalue problem given in (27) & (28) has no eigenvalue $s$ with $\Re(s)>0$, and hence the following proposition concerning the stability of the scheme holds.

Proposition 2**.**

If the divergence damping coefficient $\alpha=0$, the semi-discrete problem (25) with boundary condition (26) is stable in the sense of Godunov-Ryabenkii condition.

7.2.2 With divergence damping ($\alpha\neq 0$)

Now, we consider $\alpha\neq 0$. In this case, the pressure and velocity components are coupled. To solve the difference equations, we make the ansatz

[TABLE]

Since we want the solutions to be bounded at $y=\infty$, we look for solutions with $|\lambda|<1$.

Inserting the ansatz into the eigenvalue problem (27), we have

[TABLE]

where the following observations have been utilized,

[TABLE]

The characteristic equation of (31) is

[TABLE]

Noticing that $d_{1}^{2}=d_{2}(d_{2}+4)/4$, the characteristic equation is in fact a cubic equation for $d_{2}$; namely,

[TABLE]

So there are three roots for $d_{2}$:

[TABLE]

Note that the three roots are distinct if $\Re(s)>0$. For each $d_{2}^{(n)}$, we have an equation for $\lambda$,

[TABLE]

The $\lambda$ equation has two reciprocal roots, and the root with magnitude less than one is denoted as $\lambda_{(n)}$. The corresponding solutions for $(\tilde{U}_{0},\tilde{V}_{0},\tilde{P}_{0})^{T}$ are

[TABLE]

where $d_{1}^{(n)}=\left(\lambda_{(n)}-\lambda^{-1}_{(n)}\right)/2$.

Therefore, the general solution to (27) is

[TABLE]

After applying the boundary conditions (28), we have a system of equations for $\mathbf{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})^{T}$,

[TABLE]

where

[TABLE]

The numerical scheme given in (25) & (26) is stable in the Godunov-Ryabenkii sense if we can show the determinant condition $\det(Z)\neq 0$ holds for $\Re(s)>0$. Note that the determinant condition implies $\mathbf{\sigma}=\mathbf{0}$ and the solution (35) is trivial; thus, there is no eigenvalue $s$ with $\Re(s)>0$.

Here we provide some evidence that the determinant condition is not violated by plotting zero contours of the real and imaginary parts of $\det(Z(s)),\forall s\in\mathbb{C}$ for some wavenumbers ($k=1,5,10,100$) in Figure 2. Note that the intersections of $\Re(\det(Z))=0$ and $\Im(\det(Z))=0$ indicate the values of $s$ that makes $\det(Z)=0$. As is shown in Figure 2, no intersections are found on the right half of the complex plane (i.e.,$\Re(s)>0$) for all the examples considered.

Remark: it is hard to prove the the determinant condition for the general case. But following [44, 4], we are able to establish the so-called “local stability” for the scheme. Further, by solving the leading order terms of the error equations explicitly, we can identify that with $\alpha\sim 1/h^{2}$ the scheme is second-order accurate because the error introduced on the boundary will be damped out by the divergence damping term. However, since TN boundary condition corresponds to a first-order discretization, the boundary-layer errors in the pressure solution is expected for this case. Details of the local stability and accuracy are discussed in the following subsection.

7.3 Local Stability and Accuracy

In this section, we show stability and accuracy of the scheme assuming that $h\sqrt{s/\nu+k^{2}}\ll 1$. Stability results established under this assumption is referred to as local stability in the literature [44, 4]. A scheme that is locally stable but not stable in the global sense can be quickly identified in computations since the unstable modes occur at high frequencies.

Proposition 3**.**

Assuming that $h\sqrt{s/\nu+k^{2}}\ll 1$ and $\alpha=C_{d}/h_{0}^{2}$ for some fixed grid spacing $h_{0}$, there exists a constant $h_{c}$ such that the determinant condition $\det(Z)\neq 0$ holds for $\forall h<h_{c}$ for $\Re(s)>0$ and a fixed wavenumber $k$.

Proof.

By solving equation (34) with the assumption $\Re(s)>0$, we know the root with magnitude less one is of the following form:

[TABLE]

From (33), we have

[TABLE]

Further, we can show that

[TABLE]

Therefore, we have

[TABLE]

and hence

[TABLE]

Obviously, if $\Re(s)>0$, we have

[TABLE]

Since $\det(Z)$ continuously depends on $h$, there exists a constant $h_{c}$ such that $\det(Z)\neq 0,~{}~{}\forall h<h_{c}.$ This completes the proof. ∎

As an example, we plot the zero contours of the real and imaginary parts of $\lim_{h\rightarrow 0}\det(Z)$ with wavenumber $k=1$ and $10$ in Figure 3. As expected, no intersection is observed when $\Re(s)>0$ for both cases.

With the local stability assumption, we can further show how the TN and WABE boundary conditions affect the accuracy of the numerical results by solving the leading-order terms of the error equation (27) with the non-homogeneous boundary condition that accounts for the truncation error of the corresponding boundary conditons, i.e.,

[TABLE]

where $r=1$ for TN condition and $r=2$ for WABE condition; here $\tilde{g}_{0}=\mathcal{O}\left(1\right)$. Notice that this boundary condition is the Laplace transformation of (26). It has already been shown that the general solution to (27) is (35). After applying the no-slip boundary conditions $\tilde{U}=\tilde{V}=0$ and the boundary condition (38), we have

[TABLE]

where $Z$ is the same as (36) and its leading order is given by equation (37). The solutions to the leading order equations are

[TABLE]

Thus, we see that if $g_{0}\neq 0$, the error introduced by the boundary condition can affect the interior. However, with large $\alpha$, the boundary error is rapidly damped out producing a boundary layer of order $r$ ($r=1$ for TN and $r=2$ for WABE).

Remark: the error estimate of the IBVP (25) & (26) can be obtained in two stages. First, we obtain estimates for a pure initial-value problem on a periodic domain satisfying the forcing terms (i.e., $h^{2}F_{u}$, $h^{2}F_{u}$, $h^{2}F_{p}$); this problem is the same as the 2nd-order finite-difference scheme, and it is shown in [44] that the error estimate is $\mathcal{O}\left(h^{2}\right)$. And then, after subtracting the solutions of the pure initial-value problem from the IBVP, we have a new IBVP with zero forcing on the interior equations and inhomogeneous boundary conditions. We have already shown that the boundary errors are quickly damped out producing a boundary layer of order $r$. Therefore, after using the Parseval’s relation, we know that the scheme is $\mathcal{O}\left(h^{2}\right)$ in $L_{2}$ norm, but with the existence of numerical boundary-layer errors that is $\mathcal{O}\left(h^{r}\right)$ ($r=1$ for TN and $r=2$ for WABE). The results are supported by careful numerical mesh refinement studies shown below.

8 Numerical results

We now present the results for a series of simulations chosen to demonstrate the properties of our numerical approach. We first consider the INS equations on a sequence of refined unit square meshes to study the accuracy of the scheme; cases with and without the divergence damping are considered subject to both the TN and WABE boundary conditions. Some benchmark problems are also considered to further illustrate the numerical properties of our scheme and to compare with existing results. Finally, as a demonstration that our scheme can be easily extended to work with higher-order elements, we solve the classical flow-past-a-cylinder problem using $\mathbb{Pe}_{n}$ finite elements with $n\geq 1$.

Remark: for all the test problems with known exact solutions, errors of the numerical solutions are measured using both $L_{\infty}$ and $L_{2}$ norms. To be specific, given an exact solution $v$ and its FEM approximation $v_{h}$ in the finite space $V_{h}=\mathrm{span}\{\varphi_{1},\varphi_{2},\dots,\varphi_{n}\}$, we define the error function as

[TABLE]

where $v_{e}$ is the projection of the exact solution $v$ onto the finite space $V_{h}$, i.e., $v_{e}=\sum_{i=1}^{n}v(\mathbf{x}_{i})\varphi_{i}$. Here $\mathbf{x}_{i}$ is the coordinates of the corresponding degree of freedom. In $V_{h}$, the $L_{\infty}$ and $L_{2}$ norms of the error function is given by

[TABLE]

A numerical quadrature with sufficient order of accuracy is used to compute the integral; for example, for $\mathbb{Pe}_{1}$ elements, a third order accurate quadrature rule is used.

8.1 Manufactured Solutions

To numerically investigate the accuracy and stability of the finite element scheme, we perform careful mesh refinement study using the method of manufactured solutions [46]. Exact solutions of the INS equations can be constructed by adding forcing functions to the governing equations. The forcing is specified so that a chosen function becomes an exact solution to the forced equation. Here we use the following trigonometric functions as the exact solutions for our convergence tests,

[TABLE]

Note that the exact solutions are chosen to be divergence free. Parameters for the exact solutions are specified as $a=0.5,f_{x}=2,f_{y}=2,$ and $f_{t}=2.$

For simplicity, the test problems are solved on a unit square domain ($\Omega=[0,1]\times[0,1]$) with uniform triangular meshes using grid spacings $h=1/(10j),j=1,2,...$. The mesh with $h=1/10$ is shown in Figure 4. The INS equations are discretized using $\mathbb{Pe}_{1}$ finite element and numerically solved using the algorithm described in Section 5. We note that the discrete system obtained using $\mathbb{Pe}_{1}$ element on this particular mesh is equivalent to the 2nd-order centered finite difference scheme after lumping the mass matrix. Therefore, the numerical tests considered here are directly related to the normal-mode analysis conducted in Section 7. In order to respectively study the effects of the divergence damping and the pressure boundary conditions, we consider the following four cases: (i) $\alpha=0$ with TN boundary condition; (ii) $\alpha=1/h^{2}$ with TN boundary condition; (iii) $\alpha=0$ with WABE boundary condition; and (iv) $\alpha=1/h^{2}$ with WABE boundary condition.

8.1.1 Periodic in $x$ direction

We begin the convergence study by assuming periodicity in $x$ direction and no-slip boundary conditions on the other boundaries (i.e., $y=0$ and $y=1$). This test is designed to match the assumptions of the analysis discussed in Section 7. Results for the cases (i) TN boundary condition without divergence damping, (ii) TN boundary condition with divergence damping and (iv) WABE boundary condition with divergence damping are collected in in Figure 5. The errors for the velocity component $v$ and the pressure $p$ are plotted with the first row of images for case (i), the second row for case (ii), and the third row for case (iv). The solutions plotted here are obtained using the uniform square mesh with grid spacing $h=1/160$. Here we observe that the errors in all cases are well behaved in that the magnitudes are small and they are smooth throughout the domain except near the boundary layers for cases (ii) and (iv). We also observe that the interior accuracy is improved by including divergence damping, and the boundary accuracy is further improved by using the WABE condition, which is consistent with our analytical results.

A convergence study for all the four cases is shown in Figure 6. We see that, without divergence damping (i.e., $\alpha=0$), the errors for all the components are about first order accurate regardless of the pressure boundary conditions used. For the cases with the divergence damping turned on (i.e., $\alpha=1/h^{2}$), we observe that the accuracies for the velocity components are second order for both TN and WABE boundary conditions. However, the pressure accuracy is first order in $L_{\infty}$ norm and a little bit inferior to second order in $L_{2}$ norm if TN boundary condition is implemented for pressure; in contrast, the pressure accuracy is second order in both norms if WABE boundary condition is used. The performance of the TN boundary condition can be explained by looking at the second row of the error plots in Figure 5. We see that it improves the accuracy in the interior of the domain by adding divergence damping, so boundary layers are observed in both $E(v)$ and $E(p)$. However, since the velocity error is still dominated by the interior, we observe second order accuracy in both norms, while the error for the pressure is dominated by the boundary-layer errors so we observe first order in $L_{\infty}$ norm and almost second order in $L_{2}$ norm.

We note that the divergence $\nabla\cdot\mathbf{u}$ is always first order. This is because we simply evaluate $\nabla\cdot\mathbf{u}$ from the finite element solution for the velocity, and it is well-known that the derivative of a function represented in $\mathbb{Pe}_{1}$ elements is of first order accuracy. This can be improved by using some other post-processing techniques that compute $\nabla\cdot\mathbf{u}$ more accurately. However, since we only want to keep track of the magnitude of the divergence to make sure it remains small throughout the computation and it does not affect the accuracy of our scheme, it suffices for us to stick with this simple approach.

8.1.2 No-slip boundary conditions on all boundaries

We then consider the convergence study with no-slip boundary conditions enforced on all boundaries with the rate of convergence for all four cases shown in Figure 7. Similar convergence properties are observed for the cases (i), (ii) and (iii). For case (iv) (the lower right plot in Figure 7), we still observe second order accuracy for the solutions $u,v$ and $p$ in $L_{2}$ norm; however, we see only first order accuracy for $p$ in $L_{\infty}$ norm. By looking at the error plots in Figure 8 for this case, we see steep gradients near the corners for both $E(v)$ and $E(p)$. For $E(v)$, the error is still dominated by the interior; however, for $E(p)$, the error is dominated by the corner spikes. Thus, the max-norm error for the pressure is strongly affected by its behavior in the corners. We note that similar corner behavior for the pressure solution is also reported in [24] for numerical examples computed on domains with sharp corners. In fact, there are two independent normals at corners where two edges meet; these independent normals cause trouble when implementing Neumann boundary conditions there. Thus, we can also see the corner spikes in the error plot when solving Poisson equation with pure Neumann boundary condition on a square domain using finite element method. Fortunately, the corner issue only affects the accuracy of the scheme locally. Globally, the scheme is still well-behaved and second-order accurate, which is confirmed by the $L_{2}$-norm errors.

8.2 Modified Lid-Driven Cavity

To further verify the accuracy and efficiency of the proposed scheme, we solve a modified lid-driven cavity problem. Specifically, we consider the flow in the square domain $[0,1]\times[0,1]$. The associated boundary conditions are $\mathbf{u}=(u_{0},0)$ on the top of the domain (i.e., $y=1$) and $\mathbf{u}=(0,0)$ on the other three sides. It is well known that the classical lid-driven cavity problem that specifies $u_{0}=1$ introduces singularities at the corners since the horizontal velocity at the top corners suddenly changes from [math] to $1$. In spite of the singularities, this problem is popular for testing and evaluating numerical methods. We will compare our results with that from [47, 48, 24]. For results reported in [47, 24], the authors did nothing to suppress the corner singularities. However, a spectral method has been used in [48] so the singularities at corners cannot be ignored because the “spectral” accuracy is generally associated with the smoothness of the solution. To minimize the effect of the singularities, the authors extracted analytically the corner singularities from the dependent variables of the problem. Here we take a different approach to remove the singularities by modifying the boundary condition. Specifically, we define $u_{0}(x)$ such that its value smoothly transitions from [math] to $1$ when $x$ is away from the ends, i.e.,

[TABLE]

In the right image of Figure 9, we show the horizontal velocity (u) on part of the top boundary grids near the top-left corner (0,1) to illustrate how the boundary condition is smoothed over a couple of grid points.

The problem for $\nu=1/1000$ is solved using our numerical methods with both TN and WABE boundary conditions. The mesh used for computation consists of 4225 degrees of freedom (dof), and the maximum and minimum values of the grid spacings are $\max(h)=0.027616$ and $\min(h)=0.011561$, respectively. In Figure 9, a coarsened version of the computational mesh is shown. Note that the grids are stretched to cluster towards the boundaries for numerical purposes.

In Figure 10, we show the streamlines of the lid-driven cavity flow at $t=50$. And in Figure 11, we plot the velocity components $u$ and $v$ along the vertical and horizontal lines through the geometric center; i.e., $u(0.5,y)$ and $v(x,0.5)$. Reference data from [47] are also plotted on top of our results for comparison. We see that our results using the proposed scheme with the divergence damping coefficient $\alpha=1/h_{\min}^{2}$ match very well with existing computations reported in [47, 48, 24].

Following [48, 24], we also show the vorticity contour at levels $[-5,-4,-3,-2,-1,-0.5,0,0.5,1,2,3]$ and the pressure contour at levels $[0.3,0.17,0.12,0.11,0.09,0.07,0.05,0.02,0,-0.002]$. Plots from the solutions of both the TN and WABE boundary conditions are collected in Figure 12. In referring to the results presented in [48, 24], we find that streamlines shown in Figure 12 are in excellent qualitative agreement with the ones presented there.

8.3 Flow Past a Cylinder

As an example to illustrate the efficiency and robustness of our approach when used with higher-order finite elements, we solve a classical flow-past-cylinder problem using $\mathbb{Pe}_{n}$ finite elements with $n=1,2,4$. The settings of the test problem follow the example in [49, 50, 24]. To be specific, the domain of the problem is $\Omega=[0,2.2]\times[0,0.41]\backslash\left\{(x,y)|(x-0.2)^{2}+(y-0.2)^{2}<0.05^{2}\right\}$. The inflow and outflow velocity profiles are prescribed as a time-dependent function, $\mathbf{u}(0,y,t)=\mathbf{u}(2.2,y,t)=[0.41^{-2}\sin(\pi t/8)(6y(0.41-y)),0]^{T}.$ The top and bottom boundaries are enforced as no-slip walls. We numerically solve the INS equations using our algorithm and $\mathbb{Pe}_{n}$ finite elements in the most straightforward way; that is, the numerical solutions are represented in (7) with the basis function $\varphi_{i}\in\mathbb{Pe}_{n}$. Both the TN and WABE conditions are considered.

The coarsest computational mesh ($\mathcal{G}_{1}$) used to solve this problem is shown in the top-left image of Figure 13, which consists of 814 triangles and 486 vertices with the largest and smallest grid spacings being $\max(h)=0.1055$ and $\min(h)=0.0078$, respectively. For convergence study, we refine $\mathcal{G}_{1}$ by splitting each side of the triangles into $n$ equal segments, and denote the refined mesh as $\mathcal{G}_{n}$. We study the convergence with mesh refinement by solving the problem using ${\mathbb{Pe}}_{1}$ elements on $\mathcal{G}_{1}$, $\mathcal{G}_{2}$ and $\mathcal{G}_{4}$. We also study the convergence related to the order of the elements by solving the problem on the same mesh ($\mathcal{G}_{1}$) using finite elements with increasing orders ($\mathbb{Pe}_{1}$, $\mathbb{Pe}_{2}$ and $\mathbb{Pe}_{4}$). Note that $(\mathbb{Pe}_{n},\mathcal{G}_{k})$ indicates the element and mesh used for a particular simulation.

We show the streamlines of the flow from the simulations using $\mathbb{Pe}_{n}(n=1,2,4)$ elements in Figure 13; for comparison purposes, the results shown here maintain the same number (6828) of dofs and the same damping coefficient $\alpha=5521.08$. We can see that these solutions are consistent with each other, and our solutions are comparable with reference results given in [24].

To further validate our algorithm, we compute the drag and lift coefficients at the cylinder, denoted by $C_{d}(t)$ and $C_{l}(t)$, and the pressure difference between the front and the back of the cylinder, $\Delta p(t)$. The evolution of these variables are shown in Figure 14 and Figure 15. Here Figure 14 collects the results obtained using $\mathbb{P}_{1}$ finite elements with TN boundary condition on a sequence of refined meshes ( $\mathcal{G}_{1},\mathcal{G}_{2}$, and $\mathcal{G}_{4}$), while Figure 15 collects the results obtained using $\mathbb{P}_{1}$, $\mathbb{P}_{2}$ and $\mathbb{P}_{4}$ finite elements with WABE boundary condition on the same mesh $\mathcal{G}_{1}$. We also calculate the maximum values of $C_{d}(t)$ and $C_{l}(t)$ and the times when they occur. All the results obtained with 6828 dofs are tabulated together with reference values in the literature in Table 1. The sense of convergence for these quantities can be clearly seen in Figure 14 and Figure 15. Although we use a rather coarse triangular mesh, our results agree quite well with the range of the reference values.

We have to note that we can not claim optimal order of accuracy is achieved for $\mathbb{Pe}_{n}$ elements with $n>1$ at this point even though we have obtained accurate and comparable results in the flow-past-a-cylinder example. This is because, for all the shown computations, triangular meshes are used, which means that the cylinder is not a real cylinder but a polygon instead. Therefore, the computational mesh contributes an $\mathcal{O}\left(h^{2}\right)$ error. To achieve higher than second-order accuracy, we need to adapt our algorithms for isoparametric finite elements. And this will be left for future work.

9 Conclusion

In this paper, we propose an algorithm that solves the incompressible Navier-Stokes equations in the velocity-pressure formulation using a split-step method that separates the updates for velocity and pressure at each time step. The separation of the pressure solution is the key to avoid solving a saddle point problem whose solution depends on the choice of finite-element spaces for velocity and pressure that is subject to the LBB condition. Therefore, our algorithm has more flexibility of choosing finite-element spaces; for efficiency and robustness, Lagrange (piecewise-polynomial) finite elements of equal order for both velocity and pressure are used in the algorithm.

We also include a divergence damping term into our formulation, this linear damping term plays no role at the PDE level, but helps suppress the numerical divergence in the discretized equations, and more importantly it improves the accuracy of the scheme. Motivated by a post-processing technique that produces super-convergent derivatives from finite-elememt solutions, we formulate an alternative compatibility boundary condition at the discrete level for the pressure equation. The new pressure boundary condition, referred to as WABE boundary condition, are shown to help the pressure solution achieve better accuracy near the boundary.

An important feature of the paper is that we use the normal-mode analysis, a technique that is often used for the analysis of finite difference schemes, to reveal the stability and accuracy properties of our finite element scheme via a simplified model problem. The model problem is discretized on a uniform mesh using $\mathbb{Pe}_{1}$ finite elements, so that we can rewrite the scheme as a finite difference method and then perform the normal-mode analysis to the resulted discrete system. The analysis shows that the scheme for the model problem is locally stable with the presence of a large divergence damping term for both TN and WABE pressure boundary conditions. Further, by obtaining the leading order solutions of the error equations, we find that the error introduced by the boundary condition is rapidly damped out by the divergence damping term producing boundary-layer errors. Since the WABE boundary condition is more accurate than the TN boundary condition in terms of truncation error, it is expected to help alleviate the boundary-layer errors for the pressure solution.

Moreover, we conduct careful numerical tests to verify the stability and accuracy of our scheme. Through convergence studies using the method of manufactured solutions, we find that the interior accuracy is improved by including divergence damping, and the boundary accuracy is further improved by using the WABE pressure boundary condition. Mesh refinement study using $\mathbb{Pe}_{1}$ finite elements confirms that, with both divergence damping and WABE condition, our scheme is 2nd order accurate up-to the boundary, which is optimal for the elements used. The numerical results agree with the analysis. To further validate our scheme, benchmark problems such as lid-drive cavity and flow-past-a-cylinder are also considered; we have shown that solutions obtained using our algorithm are in excellent agreement with those reported in the literature.

The split-step finite-element method developed in this paper exhibits good numerical properties. In particular, the flexibility of using the standard Lagrange (piecewise-polynomial) finite elements of equal order for both velocity and pressure, which violates the classic LBB stability condition, makes it much easier to couple our fluid solver with a structure solver for solving FSI problems. Furthermore, the ability to suppress the boundary-layer errors in the pressure solution using WABE boundary condition ensures that information can be accurately transferred across the fluid-structure interface, which could be crucial for an FSI algorithm to maintain high-order accuracy. In the future, we will investigate extending our AMP FSI schemes using this finite element INS algorithm to develop an accurate and efficient partitioned FSI algorithm within the finite element framework.

In terms of improving the INS algorithm itself, we will use isoparametric elements to achieve higher order accuracy. In addition to the WABE boundary condition, we will also investigate the possibility of addressing the pressure boundary-layer issue using $p$-refinement by increasing the polynomial degree for the basis functions on the boundary nodes, thus we could obtain more accurate approximation of $\nabla\times\mathbf{u}_{h}$ on the boundary. This $p$-refinement strategy is also computationally efficient since the degrees of freedom are only increased on the boundary with all the other basis functions in the interior being unchanged. Analysis for this new finite element framework is also under investigation; we are interested in deriving some energy estimates for our scheme under more general assumptions.

10 Acknowledgment

The author would like to acknowledge Professors W. D. Henshaw, Fengyan Li, J. W. Banks and D. W. Schewendeman of Rensselaer Polytechnic Institute (RPI) for helpful conversations. Portions of this research were conducted with high performance computational resources provided by the Louisiana Optical Network Infrastructure (http://www.loni.org).

Appendix A

In this appendix, we show technical details for some of the results presented in Section 7.

Proposition 1

Proof.

Inserting the errors $U_{j}=u_{j}-\hat{u}_{j}$, $V_{j}=v_{j}-\hat{v}_{j}$ and $P_{j}=p_{j}-\hat{p}_{j}$ into the numerical scheme (22) and expanding the exact solutions at grid $j$, we arrive at the following error equations

[TABLE]

Similarly, for the the TN condition, we have

[TABLE]

and, for the WABE condition, we have

[TABLE]

Realizing that $Mf_{j}=f_{j}+\mathcal{O}\left(h^{2}\right)$, and representing $\mathcal{O}\left(h^{r}\right)$ as $h^{r}F$ with some $\mathcal{O}\left(1\right)$ function $F$, we derive the error equations (25) and their boundary conditions (26), which completes the proof. ∎

Note that replacing $Mf_{j}$ with $f_{j}$ is in the same spirit of mass lumping, a technique frequently used in finite-element methods. In addition, to see the second order accuracy of the WABE condition, it is most convenient to expand the exact solutions about the point $y=1/3$.

• 2.

Lemma 1

Proof.

If $\alpha=0$, the general solution to the pressure equation in the eigenvalue problem (27) is found to be

[TABLE]

where $\xi$ satisfies

[TABLE]

Substituting the pressure solution into the velocity equations, we have

[TABLE]

We note that both velocity equations have the same homogeneous part, which is solved by

[TABLE]

Here $\gamma$ satisfies

[TABLE]

Note that $\gamma=\gamma(s)$ depends on $s$. The particular solutions of the velocity equations have the forms

[TABLE]

Substituting them into (40) and (41), respectively, we get

[TABLE]

Using (39), we have

[TABLE]

The general solutions of (40) and (41) are then given by

[TABLE]

Implementing the no-slip boundary condition, $\tilde{U}_{0}=\tilde{V}_{0}=0$, we have

[TABLE]

Therefore, we have found the solution given in (29)

∎

• 3.

Lemma 2

Proof.

If $q_{1}(s)$ is real, there is a real number $c$ such that

[TABLE]

Then we have

[TABLE]

and squaring (43) implies

[TABLE]

From (39) and (42), we have

[TABLE]

After inserting (43) and (46) into (44) to eliminate $e^{-\gamma h}$ and $e^{-2\gamma h}$, we have

[TABLE]

Simplifying the above equation using (45), we arrive at

[TABLE]

This is a quadratic equation for $s$ with real coefficients. Obviously, $s=0$ is a root; it follows that the other root must also be real, which proves the lemma. ∎

• 4.

Lemma 3

Proof.

If $s$ is a root for $q(s)$, then $0=q(s)=\left(e^{-\xi h}-1\right)+\nu k^{2}q_{1}(s)$ implies $q_{1}(s)$ is real. According to Lemma 2, $s$ must be real. So it suffices to consider $s>0$ to prove this lemma.

For $s>0$, we solve (45) and (46) and obtain

[TABLE]

We note that both (45) and (46) have two roots that are reciprocal of each other. Because of the regularity conditions at infinity, the roots with magnitude less than one are kept in the above expressions. Therefore, we have

[TABLE]

Then the derivative of $q_{1}(s)$ can be obtained explicitly as following:

[TABLE]

where

[TABLE]

It is easily seen that

[TABLE]

Notice that $N_{1}>0$ and $N_{2}>0$, so we have $N_{1}-N_{2}>0$. Since $N_{3}>0$, we conclude that $q_{1}^{\prime}(s)<0$. So $q_{1}(s)$ is decreasing, and we have $q_{1}(s)<q_{1}(0),$ where

[TABLE]

The limit is evaluated using L’Hospital’s rule. Furthermore, $q^{\prime}(s)=\nu k^{2}q^{\prime}_{1}(s)<0$ implies that $q(s)<q(0).$ Here

[TABLE]

Therefore, we have shown that $q(s)<0$ for $s>0$, and this proves the lemma. ∎