A split-step finite-element method for incompressible Navier-Stokes equations with high-order accuracy up-to the boundary
AA split-step finite-element method for incompressible Navier-Stokesequations with high-order accuracy up-to the boundary
Longfei Li a,1,2, ∗ a Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504, USA.
Abstract
An efficient and accurate finite-element algorithm is described for the numerical solution of the incompress-ible Navier-Stokes (INS) equations. The new algorithm that solves the INS equations in a velocity-pressurereformulation is based on a split-step scheme in conjunction with the standard finite-element method. Thesplit-step scheme employed for the temporal discretization of our algorithm completely separates the pres-sure updates from the solution of velocity variables. When the pressure equation is formed explicitly, thealgorithm avoids solving a saddle-point problem; therefore, our algorithm has more flexibility in choosingfinite-element spaces. In contrast, popular mixed finite-element methods that solve the INS equations inthe primitive variables (or velocity-divergence formulation) lead to discrete saddle-point problems whosesolution depends on the choice of finite-element spaces for velocity and pressure that is subject to the well-known Ladyzenskaja-Babuˇska-Brezzi (LBB) (or inf-sup) condition. For efficiency and robustness, Lagrange(piecewise-polynomial) finite elements of equal order for both velocity and pressure are used. Accurate nu-merical boundary condition for the pressure equation is also investigated. Motivated by a post-processingtechnique that calculates derivatives of a finite element solution with super-convergent error estimates, analternative numerical boundary condition is proposed for the pressure equation at the discrete level. The newnumerical pressure boundary condition that can be regarded as a better implementation of the compatibilityboundary condition improves the boundary-layer errors of the pressure solution. A normal-mode analysisis performed using a simplified model problem on a uniform mesh to demonstrate the numerical propertiesof our methods. Convergence studies using P elements support the analytical results and demonstratesthat our algorithm with the new numerical boundary condition achieves the optimal second-order accuracyfor both velocity and pressure up-to the boundary. Benchmark problems are also computed and carefullycompared with existing studies. Finally, as an example to illustrate that our approach can be easily adaptedfor higher-order finite elements, we solve the classical flow-past-a-cylinder problem using P n finite elementswith n ≥ Keywords:
Navier-Stokes equations, pressure-Poisson reformulation, split-step method,finite-element method, normal-mode analysis ∗ Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504, USA
Email address: [email protected] (Longfei Li) Research supported by the Margaret A. Darrin Postdoctoral Fellowship of Rensselaer Polytechnic Institute (RPI). Research supported by the Louisiana Board of Regents Support Fund under contract No. LEQSF(2018-21)-RD-A-23.
Preprint submitted to Elsevier January 23, 2020 a r X i v : . [ m a t h . NA ] J a n ontents α = 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137.2.2 With divergence damping ( α (cid:54) = 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147.3 Local Stability and Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 x direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208.1.2 No-slip boundary conditions on all boundaries . . . . . . . . . . . . . . . . . . . . . . 208.2 Modified Lid-Driven Cavity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 248.3 Flow Past a Cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1. Introduction
The development and analysis of numerical schemes for the fast and accurate solution of incompressibleNavier-Stokes (INS) equations have long been a very active area of research, see for example [1–7] and thereferences therein. Popular existing numerical algorithms for solving INS equations include but are notlimited to the following pioneering methods and their follow-up variants: (i) the MAC method that usesstaggered 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 anequivalent velocity-pressure reformulation of the INS equations [4, 6, 7]. There are also numerous otherapproaches based on common discretizations such as finite difference, finite element, finite volume, spectralelement, and discontinuous Galerkin method; to name just a few, see [11–17].The focus here is on solving the INS equations using finite element methods (FEM). The popular onesinclude methods based on the weak formulation of the INS equations in the primitive variables (also referredto as the velocity-divergence formulation in the literature), which often employ the H (Ω) d × L (Ω) con-forming elements for spatial discretization. Here d = 2 or 3 denotes the spatial dimension and Ω representsthe fluid domain. Methods using the H (div; Ω) × L (Ω) conforming finite elements of the Raviart-Thomastype 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 orderfor this problem to have a solution, the choice of finite-element spaces for velocity and pressure must satisfythe well-known inf-sup or Ladyzenskaja-Babuˇska-Brezzi (LBB) condition; these methods are often calledmixed finite-element methods. For example, P / P is a mixed velocity/pressure finite-element pair satisfying he LBB condition. The additional complexity posed by the LBB condition makes it hard to utilize theseschemes 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 theclassical book by Girault and Raviart [2].Some projection methods based on the pioneering work of Chorin [1] avoid solving a saddle point problemby performing separate pressure updates through a splitting strategy. Therefore, the solution of these discreteproblems is not subject to the LBB condition and this suggests more flexibility in choosing finite-elementspaces for spatial discretization. However, the projection methods have their own drawback; that is, thepressure solution near boundaries is found to be less accurate with the presence of numerical boundary-layererrors [3, 5, 19]. Numerous techniques were developed attempting to reduce these numerical boundary layersin modern second order projection methods [3, 9, 20, 21]. Especially, Brown et al . [22] proposed a secondorder implicit projection method that is able to achieve full accuracy for both the velocity and pressure inrectangular domains; however, within a finite-element setting, it is not straightforward that their method isapplicable for more general geometries due to the high order spatial derivatives required in their pressureupdate. More recently, works of Liu et al. have led to the development of a series of INS algorithms usingeither a velocity-pressure reformulation of the original INS equations or a variation of projection methodsbased on the Laplacian and Leray projection operators; these methods have been demonstrated to achievepromising 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 FSIapplications. 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 basedon an interface condition derived at the continuous level by matching the time derivative of the kinematicinterface condition. The AMP condition, which is a non-standard Robin-type boundary condition involvingthe fluid stress tensor, requires no adjustable parameters and in principle is applicable at the discrete level tocouple 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, theAMP algorithms have been developed and implemented to solve FSI problems involving the interactionof incompressible flows with a wide range of structures, such as elastic beams/shells [25, 26], bulk solids[27–29] and rigid bodies [30–32]. It has been shown in these works that the AMP schemes are second-orderaccurate and stable without sub-time-step iterations, even for very light structures when added-mass effectsare strong. In contrast, the state-of-the-art finite-element based loosely-coupled partitioned FSI algorithmshave yet achieved second-order accuracy for all the solution components of the fluid-structure system [33–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 efficiencyof the state-of-the-art FSI simulations within the finite-element framework. However, the extension to FEMis nontrivial, and it poses specific requirements for the underlying INS solver. With this long-term goal inmind, 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 sta-bility condition in finite elements) so that the pressure solution is free of spurious oscillations; • appropriate boundary conditions should be prescribed for an intermediate velocity field for projectiontype methods or for the pressure if using split-step type methods based on a velocity-pressure refor-mulation of the INS equations. These boundary conditions are essential to keep the boundary-layererrors small, and their choices are non-trivial; • the scheme should be efficient and accurate, and it is often useful to decouple the solution of thevelocity from the solution of the pressure; • the scheme should be able to keep the discrete divergence small if it is not strictly enforced to be zero; • 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-stepalgorithm is developed for solving a pressure-velocity reformulation of the INS equations. The split-stepstrategy that separates the solution of pressure from that of the velocity variables enables more flexibility ofchoosing finite-element space for spatial discretization in the same manner as projection methods by avoiding saddle point problem. Therefore, Lagrange (piecewise-polynomial) finite-elements of equal order for bothvelocity and pressure can be used by our method for efficiency. It is important to note that these standardelements can not be used in a straightforward way when solving the INS equations in the velocity-divergenceformulation because they fail to satisfy the LBB condition. Therefore, our algorithm has the potential to bemore efficient than many existing FEM based fluid solvers.Special attention has been paid to investigate accurate numerical boundary conditions for the pressureequation. The curl-curl boundary condition appeared in the velocity-pressure reformulation is a compatibil-ity boundary condition, which is derived by applying the normal component of the momentum equations onthe boundary. The correct implementation of this compatibility boundary condition as a numerical bound-ary condition for the pressure is found to be crucial for the stability and accuracy of similar finite differenceschemes [4, 6, 41]. To this end, two approaches of incorporating the curl-curl boundary condition within thefinite-element context are considered that lead to two numerical boundary conditions of different accura-cies. The most straightforward way to incorporate the curl-curl boundary condition, a Neumann boundarycondition for the pressure, is to implement it as a natural boundary condition within the weak formulationof the pressure-Poisson equation. This approach, though simple and straightforward, is found to be lessaccurate since a slight degradation of the pressure accuracy near the boundary is observed in numericalexperiments. Motivated by a post-processing technique that produces a super-convergent flux from finiteelement solutions [42], we propose an alternative way to implement the compatibility boundary condition asa more accurate numerical boundary condition for the pressure that alleviates the boundary-layer errors inthe pressure solution.The remainder of the paper is organized as follows. In Section 2, we discuss a pressure-velocity refor-mulation of the INS equations. A split-step scheme consisting of a second-order accurate predictor-correctoralgorithm is described in Section 3 for the temporal discretization of the problem, and the discussion of thespatial discretization using the standard finite-element method follows in Section 4. The complete discretealgorithm is summarized in Section 5. In Section 6, novel numerical boundary conditions that help thepressure solution achieve better accuracy near the boundary are presented. The numerical properties of thescheme and the two pressure boundary conditions are analyzed for a simplified model problem discretized ona uniform mesh in Section 7. Careful numerical validations are conducted in Section 8. Finally, concludingremarks 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 con-sidered. Let Ω ⊂ R d ( d = 2 or 3) denote a bounded open domain and ∂ Ω be the boundary of Ω. Thevelocity-pressure form are given by ρ (cid:18) ∂ u ∂t + u · ∇ u (cid:19) = −∇ p + µ ∆ u + F for x ∈ Ω , t > p = − ρ ∇ u : ( ∇ u ) T + ∇ · F + α ( x ) ∇ · u for x ∈ Ω , t > B ( u , p ) = g ( x , t ) for x ∈ ∂ Ω , t > ∇ · u = 0 for x ∈ ∂ Ω , t > u ( x ,
0) = f ( x ) for x ∈ Ω , t = 0 (1e)where ∇ u : ( ∇ u ) T ≡ d (cid:88) i =1 d (cid:88) j =1 ∂u i ∂x j ∂u j ∂x i . Here u = ( u , . . . , u d ) and x = ( x , . . . , x d ) are the velocity and position in d -dimensional space; p is the pres-sure; ρ is the fluid density; µ is the coefficient of viscosity; and F ( x , t ) = ( F ( x , t ) , . . . , F d ( x , t )) is the externalforce. B ( u , p ) = g ( x , t ) represents appropriate boundary conditions with g ( x , t ) = ( g ( x , t ) , · · · , g d ( x , t )) agiven vector valued function. The given initial conditions, which are assumed to be divergence free, arerepresented by f ( x ) = ( f ( x ) , . . . , f d ( x )). quation (1b) is the pressure-Poisson equation (PPE), which is obtained by taking the divergence of themomentum equation together with the divergence free condition ( ∇ · u = 0). Note that a linear dampingterm α ( x ) ∇ · u is included in the PPE for numerical purposes. The damping term, referred to as thedivergence damping, has no effect at the continuous level since ∇ · u = 0; however, it helps to suppress thedivergence at the discrete level, since the numerical solution is not exactly divergence-free due to propertiesof the numerical approximation. As is seen later in Section 5, our numerical algorithm does not enforce thedivergence-free condition for the interior of the domain, instead this condition is only implicitly guaranteedby the curl-curl pressure boundary condition given below in equation (2). Discretization errors could resultin the growth of the divergence of the velocity in the numerical solution. Therefore, it is important to includethe divergence damping in the PPE to keep ∇ · u small for the whole computation. Alternatively, one couldperform an extra projection after every time step to map the velocity solutions into a divergence free spaceat the expense of a significant amount of additional computations.The velocity-pressure formulation of the INS equations requires an extra boundary condition, and anappropriate 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 aconsequence, a Poisson problem for the pressure needs to be solved explicitly at each time step. Thedivergence boundary condition given in (1d) is, however, not convenient to be used as a boundary conditionfor the PPE (1b) and hence an alternative condition is utilized. Following the studies in [6, 12, 41], thecurl-curl boundary condition, ∂p∂n = n · (cid:18) − ρ ∂ u ∂t − ρ u · ∇ u − µ ∇ × ∇ × u + F (cid:19) , x ∈ ∂ Ω , (2)has been used in place of (1d) 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 (1a) as acompatibility boundary condition, ∂p∂n = n · (cid:18) − ρ ∂ u ∂t − ρ u · ∇ u + µ ∆ u + F (cid:19) , and then replacing the diffusion term with the following vector identity,∆ u = ∇ ( ∇ · u ) − ∇ × ∇ × u . It is important to point out that the curl-curl boundary condition (2) has the divergence free conditionimplicitly implemented.We remark that the INS equations in the velocity-pressure form were considered in [4, 6, 43], in whichsecond- and forth-order accurate finite-difference based schemes have been developed and analyzed. Inaddition, the stability analysis of the curl-curl boundary condition in the context of a centered finite differencediscretization 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 (1) are not equivalent to those in the primitive variables (alsoknown as the velocity-divergence form) in the case of steady-state flows because solutions with constant ∇ · u are possible for the velocity-pressure form; however, ∇ · u is forced to be identically zero if the divergencedamping is present. And we will see in the analysis and results sections that including the divergencedamping 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 boundaryconditions (1c) are specifically given by u | ∂ Ω = g ( x , t ) .
3. Temporal discretization
Following [6], we use a split-step strategy to separate the updates of pressure and velocity components inthe velocity-pressure reformulation (1). The split-step method is an explicit predictor-corrector method thatconsists 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 tosolve (1) within a finite-difference framework in [6]. In this paper, we are interested in extending it tothe finite-element framework. The spatial discretization using a standard finite element method will bediscussed in section 4, and the summary of the full discrete algorithm follows in section 5.To be specific, the velocity-pressure formulation (1) of the INS equations are advanced in time using thefollowing time-stepping scheme. For simplicity, the algorithm is written for a fixed time-step, ∆ t , so that t n = n ∆ t . We note that the algorithm can be extended to a variable ∆ t with some reasonable strategiesto dynamically determine the step size and the update frequency, but this case is not considered for thescope of this paper. Given solutions ( u n − , p n − ) and ( u n , p n ) at time levels t n − and t n , we first predictthe velocity using the AB2 method, ρ u p − u n ∆ t = 32 ( L u n + F n ) − (cid:0) L u n − + F n − (cid:1) with L u = − ρ u · ∇ u − ∇ p + µ ∆ u . The pressure prediction is followed by solving the PPE with the predicted velocity solutions,∆ p p = − ρ ∇ u p : ( ∇ u p ) T + ∇ · F n +1 + α ( x ) ∇ · u p . The velocity is then corrected using the following modified AM2 method ρ u n +1 − u n ∆ t = 12 ( L u n + F n ) + 12 (cid:0) L u p + F n +1 (cid:1) . Note that the modified AM2 method is explicit since the predicted velocity u p is used on the right hand sideof the above equation. The pressure correction follows,∆ p n +1 = − ρ ∇ u n +1 : (cid:0) ∇ u n +1 (cid:1) T + ∇ · F n +1 + α ( x ) ∇ · u n +1 . We emphasize that this algorithm is stable without the use of the corrector step. Typically, the correctorstep is included since the scheme has a larger stability region than the predictor alone, and the stabilityregion includes the imaginary axis so that the scheme can be used for inviscid problems ( µ = 0). The timestep ∆ t is determined by a diffusive stability constraint (∆ t ∼ h ) 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 termof the momentum equation implicitly using a Crank-Nicholson method; the time step for the semi-implicitscheme is determined by a convective stability constraint (∆ t ∼ h ). One can refer to [41] for a proof of thetime-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 updateof pressure from that of the velocity. When the pressure equation is formed explicitly, Lagrange finite elementsof equal order can be used to discretize the velocity and pressure equations in a stable manner. This is incontrast to discretizing the velocity-divergence formulation in which case the standard Lagrange basis leadsto 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 areemployed to discretize the above time-difference scheme in space; that is, we look for finite element solutionsin the finite dimensional space V dh × Q h = P dn × P n . Here P n is the piecewise polynomial finite-element spaceof degree n that is defined on a triangulation of the domain T h (Ω). If { ϕ j } Nj =1 denote the basis functions of P n with N being the number of degrees of freedom, then a set of basis functions for P dn can be convenientlyformed as Φ kj = ϕ j e k ( k = 1 , . . . , d and j = 1 , . . . , N ), where e i ’s are the standard bases of R d . Thus, thefinite-element approximations to velocity and pressure solutions at time t n can be represented as u nh = d (cid:88) k =1 N (cid:88) j =1 u nk j Φ kj , p nh = N (cid:88) j =1 p nj ϕ j . (3) . The complete numerical algorithm To simplify discussion, the following notations are introduced for inner products defined over the domainΩ and its boundary ∂ Ω, ∀ f, g ∈ L (Ω) : ( f, g ) = (cid:90) Ω f g dX, (cid:104) f, g (cid:105) = (cid:90) ∂ Ω f g dS, ∀ f , g ∈ L (Ω) d : ( f , g ) = (cid:90) Ω f · g dX, ( ∇ f , ∇ g ) = (cid:90) Ω ∇ f : ∇ g dX. And we denote V dh the subspace of V dh that vanishes on the boundary.Given the finite element solutions ( u nh , p nh ) ∈ V dh × Q h at the current time t n , and ( u n − h , p n − h ) ∈ V dh × Q h at one previous time level t n − , the goal of the algorithm is to determine the solution at time t n +1 . Thecomplete discrete scheme using the above predictor-corrector time stepping method and the FEM spatialdiscretization is as follows. Begin predictor.Stage I - velocity prediction : We predict the velocity solution by solving for u ph ∈ V dh such that u ph ( x b ) = g n +1 ( x b ) , ∀ x b ∈ ∂ Ω ,ρ ∆ t ( u ph − u nh , v h ) = 32 [( L u nh , v h ) + ( F n , v h )] − (cid:2) ( L u n − h , v h ) + ( F n − , v h ) (cid:3) , ∀ v h ∈ V dh . (4)Here ( L u h , v h ) = − ρ ( u h ·∇ u h , v h ) − ( ∇ p h , v h ) − µ ( ∇ u h , ∇ v h ). The superscripts over the given functions F and 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 mayremove this singularity). For those boundary conditions that imply a singular Poisson problem, we addan additional constraint ( p h ,
1) = 0 that sets the mean value of pressure to zero to the PPE as a Lagrangemultiplier [6]. Specifically, we solve for p ph ∈ Q h and λ ∈ R such that, for ∀ q h ∈ Q h , − ( ∇ p ph , ∇ q h ) + λ (1 , q h ) = (cid:16) − ρ ∇ u ph : ( ∇ u ph ) T + ∇ · F n +1 + α ∇ · u ph , q h (cid:17) − (cid:28) ∂p ph ∂n , q h (cid:29) , ( p ph ,
1) = 0 . (5)In practice, we chose α to be inversely proportional to the square of the mesh spacing. For a nonuniformmesh, the minimum of the spacings, h min , is used; that is, α = C d h − . In addition, the boundary integralis given by (cid:28) ∂p ph ∂n , q h (cid:29) = (cid:28) n · (cid:18) − ρ ∂ g n +1 ∂t − ρ g n +1 · ∇ u ph + F n +1 (cid:19) , q h (cid:29) + µ (cid:104)∇ × u ph , n × ∇ q h (cid:105) , (6)which is derived by utilizing the curl-curl boundary condition (2) and the following vector identity, (cid:104) n · ∇ × ∇ × u , q (cid:105) = ( ∇ × ∇ × u , ∇ q ) = − (cid:104)∇ × u , n × ∇ q (cid:105) . (7)We note that the vector identity (7) plays an important role in the algorithm since it reduces the regularityrequirement for the admissible finite-element space that makes it possible to use P finite elements. Remark : we refer to the pressure boundary condition (6) as the traditional Neumann (TN) boundarycondition, which arises naturally from testing the PPE and integration by parts. However, a slightdegradation of the pressure accuracy near the boundary is observed in numerical experiments. Wethink this degradation stems from the fact that we need to evaluate (cid:104)∇ × u h , n × ∇ q h (cid:105) for the TNpressure boundary condition, and direct evaluation of ∇ × u h from the finite element solution for the elocity field u h is of sub-optimal order of accuracy. Motivated by Carey’s post-processing techniquethat produces a super-convergent flux from finite element solutions [42], we propose an alternativecompatibility boundary condition for the pressure at the discrete level. This new condition referred toas WABE boundary condition improves the pressure accuracy near the boundary, and its detail will bediscussed in Section 6. Begin corrector.Stage III - velocity correction : To correct the velocity, we solve for u n +1 h ∈ V dh such that u n +1 h ( x b ) = g n +1 ( x b ) , ∀ x b ∈ ∂ Ω ,ρ ∆ t (cid:0) u n +1 h − u nh , v h (cid:1) = 12 [( L u nh , v h ) + ( F n , v h )] + 12 (cid:2) ( L u ph , v h ) + ( F n +1 , v h ) (cid:3) , ∀ v h ∈ V dh . (8) Stage IV - pressure update : Finally, we obtain the updated pressure solution p n +1 h ∈ Q h by solvingthe same equation (5) and (6) but with the corrected velocity solution u n +1 h supplied on the right handside.
6. Weighted average over boundary elements (WABE) boundary condition
In this section, an alternative numerical boundary condition for the PPE is introduced, which is essentialfor our algorithm to achieve high-order accuracy up-to the boundary. We derive this numerical boundarycondition from a weighted average of the momentum equations over boundary elements; and therefore, werefer to it as WABE boundary condition for short. The WABE boundary condition avoids direct evaluationof ∇ × 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 ρ (cid:18) ∂u ∂t + u · ∇ u (cid:19) = − ∂p∂x + µ (cid:18) ∂ u ∂x − ∂ u ∂x ∂x (cid:19) , (9) ρ (cid:18) ∂u ∂t + u · ∇ u (cid:19) = − ∂p∂x + µ (cid:18) ∂ u ∂x − ∂ u ∂x ∂x (cid:19) . (10)Not that the external forcing F is omitted here to save space, which can be easily included in the WABEboundary condition derived below. On the finite dimensional space V dh × Q h , testing the momentum equations(9) and (10) with the basis function ϕ i b and integrating by parts, we have (cid:18) ∂p h ∂x , ϕ i b (cid:19) = − ρ (cid:18) ∂u h ∂t + u h · ∇ u h , ϕ i b (cid:19) − µ (cid:18) ∂u h ∂x − ∂u h ∂x , ∂ϕ i b ∂x (cid:19) + n µ (cid:28) ∂u h ∂x − ∂u h ∂x , ϕ i b (cid:29) , (11) (cid:18) ∂p h ∂x , ϕ i b (cid:19) = − ρ (cid:18) ∂u h ∂t + u h · ∇ u h , ϕ i b (cid:19) − µ (cid:18) ∂u h ∂x − ∂u h ∂x , ∂ϕ i b ∂x (cid:19) + n µ (cid:28) ∂u h ∂x − ∂u h ∂x , ϕ i b (cid:29) . (12)Here i b denotes the index of any boundary degree of freedom, i.e., ∀ P i b ∈ ∂ Ω. Since ϕ i b is nonzero only onthe elements that contain the node P i b , the equations (11) and (12) are essentially a weighted average ofthe momentum equations over those boundary elements. Let n i b = (cid:0) n i b , n i b (cid:1) be the unit outward normalvector at the boundary node P i b . To mimic the curl-curl boundary condition (2), the averaged momentumequations (11) and (12) are combined in n i b direction, and the WABE boundary condition is obtained: (cid:0) n i b · ∇ p h , ϕ i b (cid:1) = − ρ (cid:18) n i b · (cid:18) ∂ u h ∂t + u h · ∇ u h (cid:19) , ϕ i b (cid:19) + µ (cid:0) ∇ × u h , n i b × ∇ ϕ i b (cid:1) + I b , (13)where I b = µ (cid:10)(cid:0) n × n i b (cid:1) · ( ∇ × u h ) , ϕ i b (cid:11) accounts for the boundary integral resulted from the integration byparts. Here n is the normal on the edge and n i b is the normal on the boundary node. The less accurateboundary integral I b vanishes if the boundary of the domain is a straight line since n × n i b = 0; however, if the oundary is not a straight line but a smooth curve, we have n × n i b = O (cid:0) h (cid:1) and thus the error introducedby ∇ × u h will be scaled down by an order of h . In either case, the WABE boundary condition (13) shouldbe more accurate than the discrete TN boundary condition (6), and thus improves the boundary-layer errorsthat appear in the pressure solution. We note that I b can be ignored for second order accurate methods inpractice.To solve for pressure with the WABE boundary condition, we only need to replace the discrete pressureequations on the boundary nodes with the equations given by the WABE boundary condition (13). To bespecific, we solve the following modified discrete PPE at the stages of pressure update as described in theprevious section, For ∀ i such that P i ∈ Ω \ ∂ Ω and ∀ i b such that P i b ∈ ∂ Ω , − (cid:0) ∇ p n +1 h , ∇ ϕ i (cid:1) + λ (1 , ϕ i ) = (cid:16) − ρ ∇ u n +1 h : (cid:0) ∇ u n +1 h (cid:1) T + α ( x ) ∇ · u n +1 h , ϕ i (cid:17) , (cid:0) n i b · ∇ p n +1 h , ϕ i b (cid:1) = − ρ (cid:18) n i b · (cid:18) ∂ u n +1 h ∂t + u n +1 h · ∇ u n +1 h (cid:19) , ϕ i b (cid:19) + µ (cid:0) ∇ × u n +1 h , n i b × ∇ ϕ i b (cid:1) , ( p n +1 h ,
1) = 0 . (14)It is important to remark that the implementation of the WABE boundary condition is not computationallymore expensive than the TN boundary condition since there is no additional calculation needed to facilitatethe WABE condition, and all the data needed, except for ∂ u n +1 h /∂t , are directly available from the previousvelocity updates. However, a sufficiently accurate value for the purpose of the boundary condition can beobtained using a forward finite difference formula in time, i.e., ∂ u n +1 h /∂t = ( u n +1 h − u nh ) / ∆ t .
7. A model problem for analysis
In this section, we perform a normal-mode analysis on a model problem to reveal the numerical propertiesof the algorithm. The motivation that we analyze the model problem for a particular geometric domain usingfinite difference theory as opposed to using a more traditional approximation theory and energy estimationis as follows. As we know, the standard energy estimate method concerns the accuracy results in L norm,which averages the numerical errors over the entire domain. However, one of the main focuses of this paperis to address the boundary-layer errors for the pressure solution; therefore, it is more appropriate to analyzenumerical errors in L ∞ norm. Using a normal-mode analysis for a particular geometric domain, we essentiallysolve the model problem analytically. Even though this analysis does not apply to general geometries andunstructured meshes, with the analytical solution of a model problem, we are able to identify some subtletiesthat are buried inside the L averaging process, i.e., the boundary-layer errors in the pressure solution. Theidea of rewriting finite element schemes on a uniform mesh as finite difference ones for analytical purposeswas also employed by other researchers; for instance, in [47], the authors rewrote their discontinuous Galerkinschemes as finite difference ones, and then performed a Fourier type analysis since it was not easy for themto 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 intheir numerical experiments.Here Stokes equations are used as a model problem for the INS equations since the nonlinear convectionterms can be regarded as lower order terms that do not affect the stability of the scheme. To further simplifythe discussion, the model problem is assumed to be 2 π -periodic in x direction on a semi-infinite domain[0 , π ] × [0 , ∞ ]. Although the model problem drastically simplifies the original INS equations, the analysisperformed here can shed some light on the stability and accuracy of our proposed method. Specifically, weconsider the following initial-boundary value problem (IBVP) with ρ = 1 and ν = µ/ρ , ∂u∂t = − ∂p∂x + ν ∆ u + f u , for ( x, y ) ∈ [0 , π ] × [0 , ∞ ] ,∂v∂t = − ∂p∂y + ν ∆ v + f v , for ( x, y ) ∈ [0 , π ] × [0 , ∞ ] , ∆ p = ∇ f + α ∇ · u , for ( x, y ) ∈ [0 , π ] × [0 , ∞ ] . e impose the no-slip and curl-curl pressure boundary conditions at the boundary y = 0, u ( x, , t ) = v ( x, , t ) = 0 and ∂p∂y = − ν ∂ u∂x∂y + f v . The homogeneous initial conditions are imposed to complete the statement of the model problem, u ( x, y,
0) = v ( x, y,
0) = p ( x, y,
0) = 0 . Utilizing the assumption of periodicity, the IBVP can be Fourier transformed in x direction; thus, in theFourier space, we have ∂ ˆ u∂t = − ik ˆ p − νk ˆ u + ν ∂ ˆ u∂y + ˆ f u ,∂ ˆ v∂t = − ∂ ˆ p∂y − νk ˆ v + ν ∂ ˆ v∂y + ˆ f v , − k ˆ p + ∂ ˆ p∂y = ik ˆ f u + ∂ ˆ f v ∂y + αik ˆ u + α ∂ ˆ v∂y , (15)subject to the transformed boundary and initial conditions, ˆ u ( k, , t ) = ˆ v ( k, , t ) = 0 ,∂ ˆ p∂y ( k, , t ) = − νik ∂ ˆ u∂y ( k, , t ) + ˆ f v ( k, , t )ˆ u ( k, y,
0) = ˆ v ( k, y,
0) = ˆ p ( k, y,
0) = 0 . (16)Here ˆ u ( k, y, t ), ˆ v ( k, y, t ), and ˆ p ( k, y, t ) are Fourier transformations of u ( x, y, t ), v ( x, y, t ) and p ( x, y, t ) withwave number k ∈ Z . For regularity, we also require (cid:107) ˆ u (cid:107) < ∞ , (cid:107) ˆ v (cid:107) < ∞ , and (cid:107) ˆ p (cid:107) < ∞ , (17)where the norm is the standard L function norm.To analyze our numerical scheme, we discretize the above transformed equations using P finite elementson a uniform Cartesian grid, G = { y j = jh | ∀ j ∈ N } . Here j = 0 corresponds to the boundary node. Thus,the finite element approximations to ˆ u ( k, y, t ), ˆ v ( k, y, t ), and ˆ p ( k, y, t ) can be represented asˆ u h ( k, y, t ) = ∞ (cid:88) j =0 u j ( k, t ) ϕ j ( y ) , ˆ v h ( k, y, t ) = ∞ (cid:88) j =0 v j ( k, t ) ϕ j ( y ) , and ˆ p h ( k, y, t ) = ∞ (cid:88) j =0 p j ( k, t ) ϕ j ( y ) , with ϕ j ( y ) denoting the j th basis function for P that is defined by ϕ j ( y ) = h ( y − y j − ) , y ∈ [ y j − , y j ] − h ( y − y j ) , y ∈ [ y j , y j +1 ]0 , otherwise , for j > , and ϕ ( y ) = (cid:40) − h ( y − y ) , y ∈ [ y , y ]0 , otherwise . For any grid function f j , we introduce the following difference operators M f j = 16 f j − + 23 f j + 16 f j +1 , D + f j = f j +1 − f j h , D − f j = f j − f j − h , and D f j = f j +1 − f j − h , where M is an average operator, and D + , D − and D are the forward, backward and centered divideddifference operators, respectively. With these difference operators, it is easily seen that(ˆ u h , ϕ j ) = hM u j , (cid:18) ∂ ˆ u h ∂y , ϕ j (cid:19) = hD u j , and − (cid:18) ∂ ˆ u h ∂y , ∂ϕ j ∂y (cid:19) = hD + D − u j , ∀ j > . imilar relations can be found for ˆ v h and ˆ p h as well. Therefore, the finite element discretization of thetransformed equations can be regarded as a finite difference scheme, and thus can be analyzed using finitedifference techniques.Specifically, we rewrite our finite element scheme as the following finite difference scheme,for j > M ˙ u j = − ikM p j − νk M u j + νD + D − u j + ( ˆ f u , ϕ j ) /h,M ˙ v j = − D p j − νk M v j + νD + D − v j + ( ˆ f v , ϕ j ) /h, − k M p j + D + D − p j = αikM u j + αD v j + (cid:16) ik ˆ f u + ∂ ˆ f v ∂y , ϕ j (cid:17) /h, (18)subject to the no-slip boundary conditions u = v = 0, and either one of the discrete TN (19) and WABE(20) boundary conditions that are described below. Here α = C d /h is assumed where C d is a constantindependent of h .To derive the discrete TN boundary condition for the model problem, we test the pressure equation in(15) with ϕ ; that is − k (ˆ p h , ϕ ) − (cid:18) ∂ ˆ p h ∂y , ∂ϕ ∂y (cid:19) + ∂ ˆ p h ∂y ϕ (cid:12)(cid:12)(cid:12)(cid:12) ∞ = αik (ˆ u h , ϕ ) + α (cid:18) ∂ ˆ v h ∂y , ϕ (cid:19) + (cid:32) ik ˆ f u + ∂ ˆ f v ∂y , ϕ (cid:33) . With the regularity condition (17) and the curl-curl condition given in (16), we have ∂ ˆ p h ∂y ϕ (cid:12)(cid:12)(cid:12)(cid:12) ∞ = − ∂ ˆ p h ∂y ( k, , t ) ϕ (0) = νik ∂ ˆ u h ∂y ( k, , t ) − ˆ f v ( k, , t ) = νik u − u h − ˆ f v , thus the TN boundary condition can be written as D + p + νikD + u = ˆ f v + k (cid:18) p + 16 p (cid:19) h + αik (cid:18) u + 16 u (cid:19) h + α v − v (cid:32) ik ˆ f u + ∂ ˆ f v ∂y , ϕ (cid:33) . (19)On the other hand, the WABE condition for the transformed model problem is given by (cid:18) ∂ ˆ p∂y , ϕ (cid:19) = − (cid:18) ∂ ˆ v∂t , ϕ (cid:19) − νk (ˆ v, ϕ ) − νik (cid:18) ∂ ˆ u∂y , ϕ (cid:19) + ( ˆ f v , ϕ ) , which can be written as the following finite difference form, D + p + νikD + u = 2 h ( ˆ f v , ϕ ) − (cid:18)
23 ˙ v + 13 ˙ v (cid:19) − νk (cid:18) v + 13 v (cid:19) . (20) Assuming that the continuous problem given in (15) and (16) has a smooth solution, (ˆ u, ˆ v, ˆ p ), we introduceit into the difference equations (18) and obtain that the order of the truncation error for the discretizationis O (cid:0) h (cid:1) . In addition, introducing the smooth solution into the boundary conditions (19) and (20), weobtain the orders of the truncation errors for the TN and WABE boundary conditions are O ( h ) and O (cid:0) h (cid:1) ,respectively. Furthermore, we can readily write down the equations for the errors U i = u i − ˆ u i , V i = v i − ˆ v i and P i = p i − ˆ p i , where ˆ u j = ˆ u ( k, y i , t ), ˆ v j = ˆ v ( k, y i , t ) and ˆ p j = ˆ p ( k, y i , t ) are the exact solutions evaluatedat y j for wavenumber k . The consistency and error equations of the discretized problem are summarized inthe following proposition; more detail of the analysis can be found in Appendix A Proposition 1.
The numerical scheme (18) subject to either the TN (19) or the WABE (20) boundaryconditions is consistent. The discretization is formally O (cid:0) h (cid:1) for the interior equations, O ( h ) for the TNboundary condition and O (cid:0) h (cid:1) for the WABE boundary condition. The error equations for the discretizedproblem are given by for j > ˙ U j = − ikP j − νk U j + νD + D − U j + h F u , ˙ V j = − D P j − νk V j + νD + D − V j + h F v , − k P j + D + D − P j = αikU j + αD V j + h F p , (21) here F u , F v and F p are some functions of O (1) , and the boundary errors are given by (cid:40) U = V = 0 D + P + νikD + U = h r g , (22) where r = 1 for TN condition and r = 2 for WABE condition with some function g = O (1) . As a remark, we look at the order of the truncation error of the divergence damping term to see whythe choice of a large coefficient, α = C d /h , does not affect the order of the truncation error of the wholescheme. The damping term appears in the pressure equation as well as in the TN boundary condition. Inthe pressure equation, the expansion of the damping term leads up to αikM ˆ u j + αD ˆ v j = α (cid:20)(cid:18) ik ˆ u j + ∂ ˆ v j ∂y (cid:19) + 16 (cid:18) ik ∂ ˆ u∂y + ∂ ˆ v∂y (cid:19) + O (cid:0) h (cid:1)(cid:21) , ∀ j > . The continuity condition implies ik ˆ u j + ∂ ˆ v j ∂y = ik ∂ ˆ u j ∂y + ∂ ˆ v j ∂y = 0 . So the divergence damping term with the choice of α = C d /h contributes an O (cid:0) h (cid:1) error that is in linewith the accuracy of the other difference operators in the scheme. Similarly, for the divergence dampingterm in the TN boundary condition, we have αik (cid:18)
13 ˆ u + 16 ˆ u (cid:19) h + α ˆ v − ˆ v α h (cid:20) ik (cid:18)
23 ˆ u + 13 ˆ u (cid:19) + ˆ v − ˆ v h (cid:21) = α h (cid:20) ik ˆ u + ∂ ˆ v ∂y + O (cid:0) h (cid:1)(cid:21) . Thus, with the continuity condition, we see that the error contributed by the divergence damping is O ( h )that is also consistent with the accuracy of the discrete TN boundary condition. It suffices to consider the homogeneous version of the problem (21) & (22) when analyzing the stabilityof 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 mentionedin [43–45] that any dissipative time discretization can be used and the resulting fully discrete problem willbe 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-discreteproblem is stable in the sense of Godunov-Ryabenkii condition; that is, we demonstrate its stability byshowing that there is no eigenvalue s with (cid:60) ( s ) > s denoting the dual variable.To be specific, after Laplace transforming the homogeneous version of (21) & (22), we obtain the eigen-value problem, for j > s ˜ U j = − ik ˜ P j − νk ˜ U j + νD + D − ˜ U j ,s ˜ V j = − D ˜ P j − νk ˜ V j + νD + D − ˜ V j − k ˜ P j + D + D − ˜ P j = αik ˜ U j + αD ˜ V j , (23)with the boundary conditions, (cid:40) ˜ U = ˜ V = 0 ,D + ˜ P + νikD + ˜ U = 0 , (24)and the regularity condition, || ˜ U || h < ∞ , || ˜ V || h < ∞ , || ˜ P || h < ∞ . Here ( ˜ U j , ˜ V j , ˜ P j ) denote the transformed solutions and || · || h represents the discrete L norm defined on thegrid G . ote that if k = 0, we have a non-trivial solution to the eigenvalue problem (23) & (24) , i.e., ˜ U j ˜ V j ˜ P j = P for any constant ˜ P . This solution should be excluded due to the regularity condition || ˜ P || h < ∞ . This case corresponds to theundetermined constant in the pressure which is regularized in our algorithm by enforcing the mean of thepressure to be zero as is described in (5). So we proceed the stability analysis assuming k (cid:54) = 0. α = 0 ) First, let us consider α = 0. In this case, the pressure is decoupled from the velocity equations, and thesolution of the eigenvalue problem can be found explicitly, which is described in the following lemma. Lemma 1. If α = 0 , the solution of the eigenvalue problem (23) & (24) is found to be ˜ U j = − ikC p s (cid:0) e − ξy j − e − γy j (cid:1) , ˜ V j = 1 hs sinh( ξh ) C p (cid:0) e − ξy j − e − γy j (cid:1) , ˜ P j = C p e − ξy j , (25) where ξ and γ satisfy h sinh (cid:18) ξh (cid:19) = k , ξ > , h sinh (cid:18) γh (cid:19) = sν + k , (cid:60) ( γ ) > . The remaining coefficient C p will be determined by the pressure boundary condition in (24) , which implies C p h (cid:20)(cid:0) e − ξh − (cid:1) + ν k s (cid:0) e − ξh − e − γh (cid:1)(cid:21) = 0 . (26)If we let q ( s ) = 1 s (cid:0) e − ξh − e − γh (cid:1) and q ( s ) = (cid:0) e − ξh − (cid:1) + νk q ( s ) , then (26) can be written as C p q ( s ) = 0. For (cid:60) ( s ) > , we have the following lemmas that imply q ( s ) (cid:54) = 0;therefore, we conclude that C p = 0 and the solution given in (25) is trivial. Lemma 2. If q ( s ) is real, then s is real. Lemma 3. q ( s ) (cid:54) = 0 for (cid:60) ( s ) > q ( s ) and q ( s ) for s > q ( 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 (23) & (24) has no eigenvalue s with (cid:60) ( s ) >
0, and hencethe following proposition concerning the stability of the scheme holds.
Proposition 2.
If the divergence damping coefficient α = 0 , the semi-discrete problem (21) with boundarycondition (22) is stable in the sense of Godunov-Ryabenkii condition.
20 40 60 80 100 s q ( s ) k=1k=5k=10k=100 s -1-0.8-0.6-0.4-0.20 q ( s ) k=1k=5k=10k=100 Figure 1: Plots of q ( s ) (left image) and q ( s ) (right image) with h = 0 . ν = 1 for various wavenumbers. α (cid:54) = 0 ) Now, we consider α (cid:54) = 0. In this case, the pressure and velocity components are coupled. To solve thedifference equations, we make the ansatz ˜ U j ˜ V j ˜ P j = λ j ˜ U ˜ V ˜ P . Since we want the solutions to be bounded at y = ∞ , we look for solutions with | λ | < h ( s + νk ) − νd ( λ ) 0 ikh h ( s + νk ) − νd ( λ ) hd ( λ ) αikh αhd ( λ ) k h − d ( λ ) ˜ U ˜ V ˜ P = , (27)where the following observations have been utilized, hD λ j = d ( λ ) λ j , where d ( λ ) = λ − λ − ,h D + D − λ j = d ( λ ) λ j , where d ( λ ) = λ − λ − . The characteristic equation of (27) is (cid:2) h ( s + νk ) − νd (cid:3) (cid:8)(cid:2) h ( s + νk + α ) − νd (cid:3) ( k h − d ) + αh ( d − d ) (cid:9) = 0 . (28)Noticing that d = d ( d + 4) /
4, the characteristic equation is in fact a cubic equation for d ; namely, (cid:2) h ( s + νk ) − νd (cid:3) (cid:8) ( ν − αh / d − h ( s + 2 νk + α ) d + h k ( s + νk + α ) (cid:9) = 0 . So there are three roots for d : d (1)2 = h ( s + νk ) ν ,d (2)2 = h ( s + 2 νk + α ) + (cid:112) [( α + s ) + ( s + α + νk ) αh k ] h ν − αh / ,d (3)2 = h ( s + 2 νk + α ) − (cid:112) [( α + s ) + ( s + α + νk ) αh k ] h ν − αh / . (29) ote that the three roots are distinct if (cid:60) ( s ) >
0. For each d ( n )2 , we have an equation for λ , λ − (cid:16) d ( n )2 (cid:17) λ + 1 = 0 , n = 1 , , . (30)The λ equation has two reciprocal roots, and the root with magnitude less than one is denoted as λ ( n ) . Thecorresponding solutions for ( ˜ U , ˜ V , ˜ P ) T are ˜ U (1)0 ˜ V (1)0 ˜ P (1)0 = − d (1)1 hik , ˜ U (2)0 ˜ V (2)0 ˜ P (2)0 = ikd (2)1 h − ( s + νk ) + ν d (2)2 h , ˜ U (3)0 ˜ V (3)0 ˜ P (3)0 = ikαd (3)1 ( hα )[ − ( s + νk ) + ν d (3)2 h ] 1 α , where d ( n )1 = (cid:16) λ ( n ) − λ − n ) (cid:17) / ˜ U j ˜ V j ˜ P j = (cid:88) n =1 σ n λ j ( n ) ˜ U ( n )0 ˜ V ( n )0 ˜ P ( n )0 . (31)After applying the boundary conditions (24), we have a system of equations for σ = ( σ , σ , σ ) T , Zσ = , where Z = − d (1)1 h ik ikαik d (2)1 h d (3)1 hα − νik d (1)1 ( λ (1) − h (cid:34) − ( s + 2 νk ) + ν d (2)2 h (cid:35) λ (2) − h (cid:34) − ( s + 2 νk ) + ν d (3)2 h (cid:35) λ (3) − hα . (32)The numerical scheme given in (21) & (22) is stable in the Godunov-Ryabenkii sense if we can show thedeterminant condition det( Z ) (cid:54) = 0 holds for (cid:60) ( s ) >
0. Note that the determinant condition implies σ = and the solution (31) is trivial; thus, there is no eigenvalue s with (cid:60) ( s ) > Z ( s )) , ∀ s ∈ C for some wavenumbers ( k = 1 , , , (cid:60) (det( Z )) = 0 and (cid:61) (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., (cid:60) ( s ) >
0) forall the examples considered.
Remark: it is hard to prove the the determinant condition for the general case. But following [4, 44],we are able to establish the so-called “local stability” for the scheme. Further, by solving the leading orderterms of the error equations explicitly, we can identify that with α ∼ /h the scheme is second-orderaccurate 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 errorsin the pressure solution is expected for this case. Details of the local stability and accuracy are discussed inthe following subsection.
500 0 500-5000500 k=1 -500 0 500-5000500 k=5 -500 0 500-5000500 k=10 -5 0 510 -505 10 k=100 Figure 2: Zero contours of (cid:60) (det( Z )) and (cid:61) (det( Z )) with h = 0 . ν = 1 for various wavenumbers. In this section, we show stability and accuracy of the scheme assuming that h (cid:112) s/ν + k (cid:28)
1. Stabilityresults established under this assumption is referred to as local stability in the literature [4, 44]. A schemethat is locally stable but not stable in the global sense can be quickly identified in computations since theunstable modes occur at high frequencies.
Proposition 3.
Assuming that h (cid:112) s/ν + k (cid:28) and α = C d /h for some fixed grid spacing h , there existsa constant h c such that the determinant condition det( Z ) (cid:54) = 0 holds for ∀ h < h c for (cid:60) ( s ) > and a fixedwavenumber k .Proof. By solving equation (30) with the assumption (cid:60) ( s ) >
0, we know the root with magnitude less oneis of the following form: λ ( n ) = (cid:16) d ( n )2 (cid:17) − (cid:114) d ( n )2 (cid:16) d ( n )2 + 4 (cid:17) , n = 1 , , . rom (29), we have lim h → d (1)2 h = s + νk ν , lim h → d (2)2 h = s + νk + αν , lim h → d (3)2 h = k . Further, we can show thatlim h → λ ( n ) − h = − lim h → (cid:115) d ( n )2 h and lim h → d ( n )1 h = − lim h → (cid:115) d ( n )2 h . Therefore, we have lim h → Z = (cid:114) s + νk ν ik ikαik − (cid:114) s + νk + αν , − | k | α − ik ( s + νk ) ( νk − α ) (cid:114) s + νk + αν ( s + νk ) | k | α (33)and hence lim h → det( Z ) = − α ( α + s ) (cid:32) | k | (cid:114) νk + sν − k (cid:33) (cid:114) νk + α + sν . Obviously, if (cid:60) ( s ) >
0, we have lim h → det( Z ) (cid:54) = 0 . Since det( Z ) continuously depends on h , there exists a constant h c such that det( Z ) (cid:54) = 0 , ∀ h < h c . Thiscompletes the proof.As an example, we plot the zero contours of the real and imaginary parts of lim h → det( Z ) with wavenum-ber k = 1 and 10 in Figure 3. As expected, no intersection is observed when (cid:60) ( s ) > D + ˜ P + νikD + ˜ U = h r ˜ g , (34)where r = 1 for TN condition and r = 2 for WABE condition; here ˜ g = O (1). Notice that this boundarycondition is the Laplace transformation of (22). It has already been shown that the general solution to (23)is (31). After applying the no-slip boundary conditions ˜ U = ˜ V = 0 and the boundary condition (34), wehave Zσ = h r ˜ g , where Z is the same as (32) and its leading order is given by equation (33). The solutions to the leading
100 -50 0 50 100-100-50050100 k=1 -100 -50 0 50 100-100-50050100 k=10
Figure 3: Zero contours of the real and imaginary parts of lim h → (det( Z )) with α = 100 and ν = 1 for wavenumbers k = 1 (left figure) and k = 10 (right figure). order equations are σ ∼ i ˜ g h r k (cid:16) | k | − (cid:112) ( νk + α + s ) /ν (cid:17) ( α + s )( | k | (cid:112) ( νk + s ) /ν − k ) (cid:112) ( νk + α + s ) /ν ,σ ∼ − ˜ g h r ( α + s ) (cid:112) ( νk + α + s ) /ν ,σ ∼ α ˜ g h r (cid:16)(cid:112) ( νk + s ) /ν (cid:112) ( νk + α + s ) /ν − k (cid:17) ( α + s )( | k | (cid:112) ( νk + s ) /ν − k ) (cid:112) ( νk + α + s ) /ν . Thus, we see that if g (cid:54) = 0, the error introduced by the boundary condition can affect the interior. However,with large α , the boundary error is rapidly damped out producing a boundary layer of order r ( r = 1 forTN and r = 2 for WABE). Remark: the error estimate of the IBVP (21) & (22) can be obtained in two stages. First, we obtainestimates for a pure initial-value problem on a periodic domain satisfying the forcing terms (i.e., h F u , h F u , h F p ); this problem is the same as the 2nd-order finite-difference scheme, and it is shown in [44] that theerror estimate is O (cid:0) h (cid:1) . And then, after subtracting the solutions of the pure initial-value problem fromthe IBVP, we have a new IBVP with zero forcing on the interior equations and inhomogeneous boundaryconditions. We have already shown that the boundary errors are quickly damped out producing a boundarylayer of order r . Therefore, after using the Parseval’s relation, we know that the scheme is O (cid:0) h (cid:1) in L norm, but with the existence of numerical boundary-layer errors that is O ( h r ) ( r = 1 for TN and r = 2 forWABE). 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 ournumerical approach. We first consider the INS equations on a sequence of refined unit square meshes to studythe accuracy of the scheme; cases with and without the divergence damping are considered subject to boththe TN and WABE boundary conditions. Some benchmark problems are also considered to further illustratethe numerical properties of our scheme and to compare with existing results. Finally, as a demonstrationthat our scheme can be easily extended to work with higher-order elements, we solve the classical flow-past-a-cylinder problem using P n finite elements with n ≥ emark : for all the test problems with known exact solutions, errors of the numerical solutions aremeasured using both L ∞ and L norms. To be specific, given an exact solution v and its FEM approximation v h in the finite space V h = span { ϕ , ϕ , . . . , ϕ n } , we define the error function as E ( v ) = | v h − v e | , where v e is the projection of the exact solution v onto the finite space V h , i.e., v e = (cid:80) ni =1 v ( x i ) ϕ i . Here x i isthe coordinates of the corresponding degree of freedom. In V h , the L ∞ and L norms of the error functionis given by || E ( v ) || ∞ = max( E ( v )) and || E ( v ) || = (cid:18)(cid:90) Ω E ( v ) d x (cid:19) / . A numerical quadrature with sufficient order of accuracy is used to compute the integral; for example, for P elements, a third order accurate quadrature rule is used. To numerically investigate the accuracy and stability of the finite element scheme, we perform carefulmesh refinement study using the method of manufactured solutions [46]. Exact solutions of the INS equationscan be constructed by adding forcing functions to the governing equations. The forcing is specified so that achosen function becomes an exact solution to the forced equation. Here we use the following trigonometricfunctions as the exact solutions for our convergence tests, u e = a sin( f x πx ) sin( f y πy ) cos( f t πt ) ,v e = a cos( f x πx ) cos( f y πy ) cos( f t πt ) ,p e = a sin( f x πx ) cos( f y πy ) cos( f t πt ) . Note that the exact solutions are chosen to be divergence free. Parameters for the exact solutions are specifiedas a = 0 . , f x = 2 , f y = 2 , and f t = 2 . Figure 4: A unit square with uniform triangular mesh. The mesh size is h = 0 . For simplicity, the test problems are solved on a unit square domain (Ω = [0 , × [0 , h = 1 / (10 j ) , j = 1 , , ... . The mesh with h = 1 /
10 is shown inFigure 4. The INS equations are discretized using P finite element and numerically solved using thealgorithm described in Section 5. We note that the discrete system obtained using P element on thisparticular mesh is equivalent to the 2nd-order centered finite difference scheme after lumping the massmatrix. Therefore, the numerical tests considered here are directly related to the normal-mode analysisconducted in Section 7. In order to respectively study the effects of the divergence damping and the pressureboundary conditions, we consider the following four cases: (i) α = 0 with TN boundary condition; (ii) α = 1 /h with TN boundary condition; (iii) α = 0 with WABE boundary condition; and (iv) α = 1 /h withWABE boundary condition. .1.1. Periodic in x direction We begin the convergence study by assuming periodicity in x direction and no-slip boundary conditionson the other boundaries (i.e., y = 0 and y = 1). This test is designed to match the assumptions of the analysisdiscussed in Section 7. Results for the cases (i) TN boundary condition without divergence damping, (ii) TNboundary condition with divergence damping and (iv) WABE boundary condition with divergence dampingare collected in in Figure 5. The errors for the velocity component v and the pressure p are plotted with thefirst row of images for case (i), the second row for case (ii), and the third row for case (iv). The solutionsplotted here are obtained using the uniform square mesh with grid spacing h = 1 / α = 0), the errors for all the components are about first order accurate regardless of the pressureboundary conditions used. For the cases with the divergence damping turned on (i.e., α = 1 /h ), weobserve that the accuracies for the velocity components are second order for both TN and WABE boundaryconditions. However, the pressure accuracy is first order in L ∞ norm and a little bit inferior to second orderin L norm if TN boundary condition is implemented for pressure; in contrast, the pressure accuracy is secondorder in both norms if WABE boundary condition is used. The performance of the TN boundary conditioncan be explained by looking at the second row of the error plots in Figure 5. We see that it improves theaccuracy in the interior of the domain by adding divergence damping, so boundary layers are observed inboth E ( v ) and E ( p ). However, since the velocity error is still dominated by the interior, we observe secondorder accuracy in both norms, while the error for the pressure is dominated by the boundary-layer errors sowe observe first order in L ∞ norm and almost second order in L norm.We note that the divergence ∇ · u is always first order. This is because we simply evaluate ∇ · u from thefinite element solution for the velocity, and it is well-known that the derivative of a function represented in P elements is of first order accuracy. This can be improved by using some other post-processing techniquesthat compute ∇ · u more accurately. However, since we only want to keep track of the magnitude of thedivergence to make sure it remains small throughout the computation and it does not affect the accuracy ofour scheme, it suffices for us to stick with this simple approach. We then consider the convergence study with no-slip boundary conditions enforced on all boundarieswith the rate of convergence for all four cases shown in Figure 7. Similar convergence properties are observedfor the cases (i), (ii) and (iii). For case (iv) (the lower right plot in Figure 7), we still observe second orderaccuracy for the solutions u, v and p in L norm; however, we see only first order accuracy for p in L ∞ 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 dominatedby the corner spikes. Thus, the max-norm error for the pressure is strongly affected by its behavior in thecorners. We note that similar corner behavior for the pressure solution is also reported in [24] for numericalexamples computed on domains with sharp corners. In fact, there are two independent normals at cornerswhere two edges meet; these independent normals cause trouble when implementing Neumann boundaryconditions there. Thus, we can also see the corner spikes in the error plot when solving Poisson equationwith pure Neumann boundary condition on a square domain using finite element method. Fortunately, thecorner issue only affects the accuracy of the scheme locally. Globally, the scheme is still well-behaved andsecond-order accurate, which is confirmed by the L -norm errors. α = 0 with TN boundary conditioncase (ii): α = 1 /h with TN boundary conditioncase (iv): α = 1 /h with WABE boundary condition Figure 5: Absolute values of the errors for the velocity component v and the pressure p are plotted, i.e., E ( v ) = | v − v e | and E ( p ) = | p − p e | . Periodicity in x direction is assumed, while no-slip boundary conditions are enforced on theboundaries in y direction. The grid spacing is h = 1 /
160 and the time for plot is t = 0 . -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || case (i): α = 0 with TN BC case (ii): α = 1 /h with TN BCcase (iii): α = 0 with WABE BC case (iv): α = 1 /h with WABE BC Figure 6: Convergence rates of u, v, p and ∇ · u are plotted at t = 0 .
1. The errors are measured using both L ∞ and L norms with the former represented by solid lines and the latter by dashed lines. Periodicity in x direction isassumed, while no-slip boundary conditions are enforced on the boundaries in y direction. -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || -3 -2 -1 e rr o r -5 -4 -3 -2 -1 || E ( u ) || ∞ || E ( v ) || ∞ || E ( p ) || ∞ || E ( ∇ · u ) || ∞ || E ( u ) || || E ( v ) || || E ( p ) || || E ( ∇ · u ) || case (i): α = 0 with TN BC case (ii): α = 1 /h with TN BCcase (iii): α = 0 with WABE BC case (iv): α = 1 /h with WABE BC Figure 7: Convergence rates of u, v, p and ∇ · u are plotted at t = 0 .
1. The errors are measured using both L ∞ and L norms with the former represented by solid lines and the latter by dashed lines. No-slip boundary conditions areenforced on all the boundaries. case (iv): α = 1 /h with WABE BC Figure 8: Absolute values of the errors for the velocity component v and the pressure p are plotted, i.e., E ( v ) = | v − v e | and E ( p ) = | p − p e | . We enforce no-slip boundary conditions at all the boundaries. The grid spacing is h = 1 / t = 0 . .2. Modified Lid-Driven Cavity To further verify the accuracy and efficiency of the proposed scheme, we solve a modified lid-driven cavityproblem. Specifically, we consider the flow in the square domain [0 , × [0 , u = ( u ,
0) on the top of the domain (i.e., y = 1) and u = (0 ,
0) on the other three sides.It is well known that the classical lid-driven cavity problem that specifies u = 1 introduces singularities atthe corners since the horizontal velocity at the top corners suddenly changes from 0 to 1. In spite of thesingularities, this problem is popular for testing and evaluating numerical methods. We will compare ourresults with that from [24, 47, 48]. For results reported in [24, 47], the authors did nothing to suppress thecorner singularities. However, a spectral method has been used in [48] so the singularities at corners cannotbe 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 fromthe dependent variables of the problem. Here we take a different approach to remove the singularities bymodifying the boundary condition. Specifically, we define u ( x ) such that its value smoothly transitionsfrom 0 to 1 when x is away from the ends, i.e., u ( x ) = 12 (cid:20) − tanh (cid:18) | x − . | − . . (cid:19) + 1 (cid:21) . In the right image of Figure 9, we show the horizontal velocity (u) on part of the top boundary grids nearthe top-left corner (0,1) to illustrate how the boundary condition is smoothed over a couple of grid points.The problem for ν = 1 / h ) = 0 . h ) = 0 . Figure 9: Left: a coarsened version of the computational mesh with grids stretched to cluster towards the boundaries.For the actual mesh used for computation, we have dof = 4225,max( h ) = 0 . h ) = 0 . In Figure 10, we show the streamlines of the lid-driven cavity flow at t = 50. And in Figure 11, we plotthe velocity components u and v along the vertical and horizontal lines through the geometric center; i.e., u (0 . , y ) and v ( x, . α = 1 /h matchvery well with existing computations reported in [24, 47, 48].Following [24, 48], we also show the vorticity contour at levels [ − , − , − , − , − , − . , , . , , ,
3] andthe pressure contour at levels [0 . , . , . , . , . , . , . , . , , − . Figure 10: Streamline plots at t = 50. Left: TN boundary condition. Right: WABE boundary condition. u(0.5,y)u(0.5,y) (Giah et al.)v(x,0.5)v(x,0.5) (Giah et al.) u(0.5,y)u(0.5,y) (Giah et al.)v(x,0.5)v(x,0.5) (Giah et al.) Figure 11: Plots for u (0 . , y ) and v ( x, .
5) at t = 50. Left: TN boundary condition. Right: WABE boundarycondition. As an example to illustrate the efficiency and robustness of our approach when used with higher-orderfinite elements, we solve a classical flow-past-cylinder problem using P n finite elements with n = 1 , ,
4. Thesettings of the test problem follow the example in [24, 49, 50]. To be specific, the domain of the problemis Ω = [0 , . × [0 , . \ (cid:8) ( x, y ) | ( x − . + ( y − . < . (cid:9) . The inflow and outflow velocity profilesare prescribed as a time-dependent function, u (0 , y, t ) = u (2 . , y, t ) = [0 . − sin( πt/ y (0 . − y )) , T . The top and bottom boundaries are enforced as no-slip walls. We numerically solve the INS equations usingour algorithm and P n finite elements in the most straightforward way; that is, the numerical solutions arerepresented in (3) with the basis function ϕ i ∈ P n . Both the TN and WABE conditions are considered.The coarsest computational mesh ( G ) used to solve this problem is shown in the top-left image ofFigure 13, which consists of 814 triangles and 486 vertices with the largest and smallest grid spacings beingmax( h ) = 0 . h ) = 0 . G by splitting eachside of the triangles into n equal segments, and denote the refined mesh as G n . We study the convergence withmesh refinement by solving the problem using P elements on G , G and G . We also study the convergencerelated to the order of the elements by solving the problem on the same mesh ( G ) using finite elements withincreasing orders ( P , P and P ). Note that ( P n , G k ) indicates the element and mesh used for a particularsimulation. e show the streamlines of the flow from the simulations using P n ( n = 1 , ,
4) elements in Figure 13;for comparison purposes, the results shown here maintain the same number (6828) of dofs and the samedamping coefficient α = 5521 .
08. We can see that these solutions are consistent with each other, and oursolutions are comparable with reference results given in [24].To further validate our algorithm, we compute the drag and lift coefficients at the cylinder, denotedby C d ( t ) and C l ( t ), and the pressure difference between the front and the back of the cylinder, ∆ p ( t ). Theevolution of these variables are shown in Figure 14 and Figure 15. Here Figure 14 collects the results obtainedusing P finite elements with TN boundary condition on a sequence of refined meshes ( G , G , and G ), whileFigure 15 collects the results obtained using P , P and P finite elements with WABE boundary conditionon the same mesh G . We also calculate the maximum values of C d ( t ) and C l ( t ) and the times when theyoccur. All the results obtained with 6828 dofs are tabulated together with reference values in the literaturein 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 referencevalues.We have to note that we can not claim optimal order of accuracy is achieved for P n elements with n > O (cid:0) h (cid:1) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - . - . - . - . - . - . - . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - . - . - . - . - . - . - . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . WABE:TN:
Figure 12: Vorticity and pressure contours at t = 50. G : ( P , G ):( P , G ): ( P , G ): Figure 13: The computational mesh (top left) and streamlines of the solutions obtained using ( P , G ) (top right),( P , G ) (bottom left) and( P , G ) (bottom right) elements with ν = 1 × − at t = 8. C d ( t ) C l ( t ) ∆ p ( t ) Figure 14: Simulations obtained using P finite elements with TN boundary condition on a sequence of refined meshes: G , G , and G . C d ( t ) C l ( t ) ∆ p ( t ) Figure 15: Simulations obtained using P , P and P finite elements with WABE boundary condition on the mesh G . error. To achieve higher than second-order accuracy, we need to adapt our algorithms for isoparametricfinite elements. And this will be left for future work. d, max t ( C d, max ) C l, max t ( C l, max ) ∆ p (8)TN ( P , G ) 2 . . . . − . P , G ) 2 . . . . − . P , G ) 2 . . . . − . P , G ) 2 . . . . − . P , G ) 2 . . . . − . P , G ) 2 . . . . − . . . . . − . et al. [24] 2 . . . . − . et al. [49] [2 . , . . , . − . , − . Table 1: Maximum values of the drag and lift coefficients and the pressure difference at final time t = 8.
9. Conclusion
In this paper, we propose an algorithm that solves the incompressible Navier-Stokes equations in thevelocity-pressure formulation using a split-step method that separates the updates for velocity and pressureat each time step. The separation of the pressure solution is the key to avoid solving a saddle point problemwhose solution depends on the choice of finite-element spaces for velocity and pressure that is subject to theLBB condition. Therefore, our algorithm has more flexibility of choosing finite-element spaces; for efficiencyand robustness, Lagrange (piecewise-polynomial) finite elements of equal order for both velocity and pressureare used in the algorithm.We also include a divergence damping term into our formulation, this linear damping term plays no role atthe PDE level, but helps suppress the numerical divergence in the discretized equations, and more importantlyit 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 boundarycondition at the discrete level for the pressure equation. The new pressure boundary condition, referredto as WABE boundary condition, are shown to help the pressure solution achieve better accuracy near theboundary.An important feature of the paper is that we use the normal-mode analysis, a technique that is oftenused for the analysis of finite difference schemes, to reveal the stability and accuracy properties of our finiteelement scheme via a simplified model problem. The model problem is discretized on a uniform mesh using P finite elements, so that we can rewrite the scheme as a finite difference method and then perform thenormal-mode analysis to the resulted discrete system. The analysis shows that the scheme for the modelproblem is locally stable with the presence of a large divergence damping term for both TN and WABEpressure boundary conditions. Further, by obtaining the leading order solutions of the error equations, wefind that the error introduced by the boundary condition is rapidly damped out by the divergence dampingterm producing boundary-layer errors. Since the WABE boundary condition is more accurate than the TNboundary condition in terms of truncation error, it is expected to help alleviate the boundary-layer errorsfor the pressure solution.Moreover, we conduct careful numerical tests to verify the stability and accuracy of our scheme. Throughconvergence studies using the method of manufactured solutions, we find that the interior accuracy is im-proved by including divergence damping, and the boundary accuracy is further improved by using the WABEpressure boundary condition. Mesh refinement study using P finite elements confirms that, with both diver-gence damping and WABE condition, our scheme is 2nd order accurate up-to the boundary, which is optimalfor the elements used. The numerical results agree with the analysis. To further validate our scheme, bench-mark problems such as lid-drive cavity and flow-past-a-cylinder are also considered; we have shown thatsolutions 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. Inparticular, the flexibility of using the standard Lagrange (piecewise-polynomial) finite elements of equal rder for both velocity and pressure, which violates the classic LBB stability condition, makes it much easierto couple our fluid solver with a structure solver for solving FSI problems. Furthermore, the ability tosuppress the boundary-layer errors in the pressure solution using WABE boundary condition ensures thatinformation can be accurately transferred across the fluid-structure interface, which could be crucial foran FSI algorithm to maintain high-order accuracy. In the future, we will investigate extending our AMPFSI schemes using this finite element INS algorithm to develop an accurate and efficient partitioned FSIalgorithm within the finite element framework.In terms of improving the INS algorithm itself, we will use isoparametric elements to achieve higherorder accuracy. In addition to the WABE boundary condition, we will also investigate the possibility ofaddressing the pressure boundary-layer issue using p -refinement by increasing the polynomial degree for thebasis functions on the boundary nodes, thus we could obtain more accurate approximation of ∇ × u h onthe boundary. This p -refinement strategy is also computationally efficient since the degrees of freedom areonly increased on the boundary with all the other basis functions in the interior being unchanged. Analysisfor this new finite element framework is also under investigation; we are interested in deriving some energyestimates for our scheme under more general assumptions.
10. Acknowledgment
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 − ˆ u j , V j = v j − ˆ v j and P j = p j − ˆ p j into the numerical scheme(18) and expanding the exact solutions at grid j , we arrive at the following error equationsfor j > M ˙ U j = − ikM P j − νk M U j + νD + D − U j + O (cid:0) h (cid:1) ,M ˙ V j = − D P j − νk M V j + νD + D − V j + O (cid:0) h (cid:1) , − k M P j + D + D − P j = αikM U j + αD V j + O (cid:0) h (cid:1) . Similarly, for the the TN condition, we have D + P + νikD + U = k (cid:18) P + 16 P (cid:19) h + αik (cid:18) U + 16 U (cid:19) h + α V − V O ( h ) , and, for the WABE condition, we have D + P + νikD + U = − (cid:18)
23 ˙ V + 13 ˙ V (cid:19) − νk (cid:18) V + 13 V (cid:19) + O (cid:0) h (cid:1) . Realizing that
M f j = f j + O (cid:0) h (cid:1) , and representing O ( h r ) as h r F with some O (1) function F , wederive the error equations (21) and their boundary conditions (22), which completes the proof.Note that replacing M f j with f j is in the same spirit of mass lumping, a technique frequently usedin finite-element methods. In addition, to see the second order accuracy of the WABE condition, it ismost convenient to expand the exact solutions about the point y = 1 / • Lemma 1 roof. If α = 0, the general solution to the pressure equation in the eigenvalue problem (23) is foundto be ˜ P j = C p e − ξy j , where ξ satisfies 4 h sinh (cid:18) ξh (cid:19) = k and ξ > . (A.1)Substituting the pressure solution into the velocity equations, we have (cid:0) s + νk (cid:1) ˜ U j − νD + D − ˜ U j = − ikC p e − ξy j , (A.2) (cid:0) s + νk (cid:1) ˜ V j − νD + D − ˜ V j = 1 h sinh( ξh ) C p e − ξy j . (A.3)We note that both velocity equations have the same homogeneous part, which is solved by˜ U hj = ˜ V hj = e − γy j . Here γ satisfies 4 h sinh (cid:18) γh (cid:19) = sν + k and (cid:60) ( γ ) > . (A.4)Note that γ = γ ( s ) depends on s . The particular solutions of the velocity equations have the forms˜ U pj = A u e − ξy j and ˜ V pj = A v e − ξy j . Substituting them into (A.2) and (A.3), respectively, we get A u (cid:20)(cid:0) s + νk (cid:1) − ν h sinh (cid:18) ξh (cid:19)(cid:21) = − ikC p ,A v (cid:20)(cid:0) s + νk (cid:1) − ν h sinh (cid:18) ξh (cid:19)(cid:21) = 1 h sinh( ξh ) C p . Using (A.1), we have A u = − ikC p s and A v = 1 hs sinh( ξh ) C p . The general solutions of (A.2) and (A.3) are then given by˜ U j = ˜ U pj + C u ˜ U hj = A u e − ξy j + C u e − γy j , ˜ V j = ˜ V pj + C v ˜ V hj = A v e − ξy j + C v e − γy j . Implementing the no-slip boundary condition, ˜ U = ˜ V = 0, we have C u = − A u and C v = − A v . Therefore, we have found the solution given in (25) • Lemma 2
Proof. If q ( s ) is real, there is a real number c such that q ( s ) = 1 s (cid:0) e − ξh − e − γh (cid:1) = c. Then we have e − γh = e − ξh − cs, (A.5) nd squaring (A.5) implies e − γh = c s − cse − ξh + e − ξh (A.6)From (A.1) and (A.4), we have e − ξh − (2 + h k ) e − ξh + 1 = 0 (A.7) e − γh − [2 + h ( s/ν + k )] e − γh + 1 = 0 (A.8)After inserting (A.5) and (A.8) into (A.6) to eliminate e − γh and e − γh , we have[2 + h ( s/ν + k )]( e − ξh − cs ) − c s − cse − ξh + e − ξh . Simplifying the above equation using (A.7), we arrive at c s − cse − ξh + (2 + h k ) cs − h s/νe − ξh + ch s /ν = 0 . This is a quadratic equation for s with real coefficients. Obviously, s = 0 is a root; it follows that theother root must also be real, which proves the lemma. • Lemma 3
Proof. If s is a root for q ( s ), then 0 = q ( s ) = (cid:0) e − ξh − (cid:1) + νk q ( s ) implies q ( s ) is real. Accordingto Lemma 2, s must be real. So it suffices to consider s > s >
0, we solve (A.7) and (A.8) and obtain e − ξh = 12 (cid:104) (2 + h k ) − (cid:112) h k + h k (cid:105) ,e − γh = 12 (cid:104) (2 + h ( s/ν + k )) − (cid:112) h ( s/ν + k ) + h ( s/ν + k ) (cid:105) . We note that both (A.7) and (A.8) have two roots that are reciprocal of each other. Because ofthe regularity conditions at infinity, the roots with magnitude less than one are kept in the aboveexpressions. Therefore, we have q ( s ) = 1 s (cid:0) e − ξh − e − γh (cid:1) = 12 s (cid:16) − h s/ν − (cid:112) h k + h k + (cid:112) h ( s/ν + k ) + h ( s/ν + k ) (cid:17) . Then the derivative of q ( s ) can be obtained explicitly as following: q (cid:48) ( s ) = − N − N N where N = 2 h s + 4 h k ν + h k ν + h k s,N = ν (cid:112) h k + h k (cid:112) h ( s/ν + k ) + h ( s/ν + k ) ,N = 2 νs (cid:112) h ( s/ν + k ) + h ( s/ν + k ) . It is easily seen that N − N = 4 h s > . Notice that N > N >
0, so we have N − N >
0. Since N >
0, we conclude that q (cid:48) ( s ) < q ( s ) is decreasing, and we have q ( s ) < q (0) , where q (0) := lim s → + q ( s ) = h k + 2 h ν √ h k + h k − h /ν. The limit is evaluated using L’Hospital’s rule. Furthermore, q (cid:48) ( s ) = νk q (cid:48) ( s ) < q ( s )
0, and this proves the lemma. eferences [1] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comp. 22 (1968) 745–762.[2] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations. Theory and algorithms,Vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, Germany, 1986.[3] S. A. Orszag, M. Israeli, M. O. Deville, Boundary conditions for incompressible flows, J. Sci. Comp.1 (1) (1986) 75–111.[4] W. D. Henshaw, H.-O. Kreiss, L. G. M. Reyna, A fourth-order accurate difference approximation forthe incompressible Navier-Stokes equations, Comput. Fluids 23 (4) (1994) 575–593.[5] W. E, J.-G. Liu, Projection method I: Convergence and numerical boundary layers, SIAM J. Numer.Anal. 32 (4) (1995) 1017–1057.[6] W. D. Henshaw, N. A. Petersson, A split-step scheme for the incompressible Navier-Stokes equations, in:M. M. Hafez (Ed.), Numerical Simulation of Incompressible Flows, World Scientific, 2003, pp. 108–125.[7] H. Johnston, J.-G. Liu, Accurate, stable and efficient Navier-Stokes solvers based on explicit treatmentof the pressure term, J. Comput. Phys. 199 (2004) 221–259.[8] F. H. Harlow, J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluidwith free surface, J. Comput. Phys. 8 (12) (1965) 2182–2189.[9] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations, J.Comput. Phys. 59 (1985) 308–323.[10] A. J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys.2 (1967) 12–26.[11] S. Abdallah, Numerical solutions for the pressure poisson equation with Neumann boundary conditionsusing a non-staggered grid, I, J. Comput. Phys. 70 (1) (1987) 182–192.[12] G. E. Karniadakis, M. Israeli, S. A. Orszag, High-order splitting methods for the incompressible Navier-Stokes equations, J. Comput. Phys. 97 (1991) 414–443.[13] J. A. Wright, W. Shyy, A pressure-based composite grid method for the Navier-Stokes equations, J.Comput. Phys. 107 (1993) 225–238.[14] K. Z. Korczak, A. T. Patera, An isoparametric spectral element method for solution of the Navier-Stokesequations in complex geometry, J. Comput. Phys. 62 (1986) 361–382.[15] J. C. Strikwerda, Finite difference methods for the Stokes and Navier-Stokes equations, SIAM J. Sci.Stat. Comput. 5 (1) (1984) 56–68.[16] W. Heinrichs, Splitting techniques for the unsteady Stokes equations, SIAM J. of Numer. Anal. 35(1998) 1646–1662.[17] B. Cockburn, G. Kanschat, D. Sch¨otzau, A locally conservative LDG method for the incompressibleNavier-Stokes equations, Math. Comp. 74 (2005) 1067–1095.[18] J. Wang, X. Ye, New finite element methods in computational fluid dynamics by H (div) elements,SIAM J. Numer. Anal. 45 (3) (2007) 1269–1286.[19] W. E, J.-G. Liu, Projection method II: Godunov-Ryabenki analysis, SIAM J. Numer. Anal. 33 (4)(1996) 1597–1621.[20] J. van Kan, A second-order accurate pressure-correction scheme for viscous incompressible flow, SIAMJ. Sci. Stat. Comput. 7 ((1986)) 870 – 891.
21] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible Navier-Stokesequations, J. Comput. Phys. 85 (1989) 257–283.[22] D. L. Brown, R. Cortez, M. L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2001) 464–499.[23] J.-G. Liu, J. Liu, R. L. Pego, Stability and convergence of efficient Navier-Stokes solvers via a commu-tator estimate, Commun. Pur. Appl. Math. 60 (10) (2007) 1443–1487.[24] J.-G. Liu, J. Liu, R. L. Pego, Stable and accurate pressure approximation for unsteady incompressibleviscous flow, J. Comput. Phys. 229 (9) (2010) 3428 – 3453.[25] J. W. Banks, W. D. Henshaw, D. W. Schwendeman, An analysis of a new stable partitioned algorithmfor FSI problems. Part I: Incompressible flow and elastic solids, J. Comput. Phys. 269 (2014) 108–137.[26] L. Li, W. D. Henshaw, J. W. Banks, D. W. Schwendeman, G. A. Main, A stable partitioned FSIalgorithm for incompressible flow and deforming beams, J. Comput. Phys. 312 (2016) 272–306.[27] J. W. Banks, W. D. Henshaw, D. W. Schwendeman, An analysis of a new stable partitioned algorithm forFSI problems. Part II: Incompressible flow and structural shells, J. Comput. Phys. 268 (2014) 399–416.[28] D. A. Serino, J. W. Banks, W. D. Henshaw, D. W. Schwendeman, A stable added-mass parti-tioned (AMP) algorithm for elastic solids and incompressible flow: Model problem analysis, preprintarXiv:1812.03192, submitted for publication (2019).[29] D. A. Serino, J. W. Banks, W. D. Henshaw, D. W. Schwendeman, A stable added-mass partitioned(AMP) algorithm for elastic solids and incompressible flow, preprint arXiv:1812.05208, submitted forpublication (2019).[30] J. W. Banks, W. D. Henshaw, D. W. Schwendeman, Q. Tang, A stable partitioned FSI algorithm forrigid bodies and incompressible flow. Part I: Model problem analysis, J. Comput. Phys. 343 (2017)432–468.[31] J. W. Banks, W. D. Henshaw, D. W. Schwendeman, Q. Tang, A stable partitioned FSI algorithmfor rigid bodies and incompressible flow. Part II: General formulation, J. Comput. Phys. 343 (2017)469–500.[32] J. W. Banks, W. D. Henshaw, D. W. Schwendeman, Q. Tang, A stable partitioned FSI algorithm forrigid bodies and incompressible flow in three dimensions, J. Comput. Phys. 373 (2018) 455–492.[33] P. Causin, J. F. Gerbeau, F. Nobile, Added-mass effect in the design of partitioned algorithms forfluid–structure problems, Comput. Method. Appl. Mech. Engrg. 194 (2005) 4506–4527.[34] M. A. Fern´andez, Coupling schemes for incompressible fluid–structure interaction: implicit, semi-implicit and explicit, SeMa Journal 55 (1) (2011) 59–108.[35] J. Degroote, Partitioned simulation of fluid–structure interaction, Archives of Computational Methodsin Engineering 20 (3) (2013) 185–238.[36] M. Bukaˇc, S. ˇCaniˇc, B. Muha, A partitioned scheme for fluid–composite structure interaction problems,Journal of Computational Physics 281 (2015) 493 – 517.[37] G. Guidoboni, R. Glowinski, N. Cavallini, S. ˇCani´c, Stable loosely-coupled-type algorithm for fluid–structure interaction in blood flow, J. Comput. Phys. 228 (18) (2009) 6916–6937.[38] S. ˇCani´c, B. Muha, M. Bukaˇc, Stability of the kinematically coupled β -scheme for fluid–structureinteraction problems in hemodynamics, ArXiv e-print 1205.6887.
39] M. A. Fern´andez, Incremental displacement-correction schemes for the explicit coupling of a thin struc-ture with an incompressible fluid, Comptes Rendus Mathematique 349 (78) (2011) 473–477.[40] M. A. Fern´andez, J. Mullaert, M. Vidrascu, Explicit Robin-Neumann schemes for the coupling ofincompressible fluids with thin-walled structures, Comput. Method. Appl. Mech. Engrg. 267 (2013)566–593.[41] N. A. Petersson, Stability of pressure boundary conditions for Stokes and Navier-Stokes equations, J.Comput. Phys. 172 (1) (2001) 40–70.[42] G. Carey, Derivative calculation from finite element solutions, Comput. Method. Appl. Mech. Engrg.35 (1) (1982) 1 – 14.[43] W. D. Henshaw, A fourth-order accurate method for the incompressible Navier-Stokes equations onoverlapping grids, J. Comput. Phys. 113 (1) (1994) 13–25.[44] W. D. Henshaw, H.-O. Kreiss, Analysis of a difference approximation for the incompressible Navier-Stokes equations, Research Report LA-UR-95-3536, Los Alamos National Laboratory (1995).[45] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, John Wileyand Sons Inc., 1995.[46] P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa Publish-ers, New Mexico, 1998.[47] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokesequations and a multigrid method, J. Comput. Phys. 48 (1982) 387 – 411.[48] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. Fl. 27 (1998)421 – 433.[49] M. Sch¨afer, S. Turek, F. Durst, E. Krause, R. Rannacher, Benchmark computations of laminar flowaround a cylinder, In: Hirschel E.H. (eds) Flow Simulation with High-Performance Computers II. Noteson Numerical Fluid Mechanics (NNFM) 48.[50] V. John, Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder,Int. J. Numer. Meth. Fl. 44 (2004) 777 – 788.39] M. A. Fern´andez, Incremental displacement-correction schemes for the explicit coupling of a thin struc-ture with an incompressible fluid, Comptes Rendus Mathematique 349 (78) (2011) 473–477.[40] M. A. Fern´andez, J. Mullaert, M. Vidrascu, Explicit Robin-Neumann schemes for the coupling ofincompressible fluids with thin-walled structures, Comput. Method. Appl. Mech. Engrg. 267 (2013)566–593.[41] N. A. Petersson, Stability of pressure boundary conditions for Stokes and Navier-Stokes equations, J.Comput. Phys. 172 (1) (2001) 40–70.[42] G. Carey, Derivative calculation from finite element solutions, Comput. Method. Appl. Mech. Engrg.35 (1) (1982) 1 – 14.[43] W. D. Henshaw, A fourth-order accurate method for the incompressible Navier-Stokes equations onoverlapping grids, J. Comput. Phys. 113 (1) (1994) 13–25.[44] W. D. Henshaw, H.-O. Kreiss, Analysis of a difference approximation for the incompressible Navier-Stokes equations, Research Report LA-UR-95-3536, Los Alamos National Laboratory (1995).[45] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time Dependent Problems and Difference Methods, John Wileyand Sons Inc., 1995.[46] P. J. Roache, Verification and Validation in Computational Science and Engineering, Hermosa Publish-ers, New Mexico, 1998.[47] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokesequations and a multigrid method, J. Comput. Phys. 48 (1982) 387 – 411.[48] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. Fl. 27 (1998)421 – 433.[49] M. Sch¨afer, S. Turek, F. Durst, E. Krause, R. Rannacher, Benchmark computations of laminar flowaround a cylinder, In: Hirschel E.H. (eds) Flow Simulation with High-Performance Computers II. Noteson Numerical Fluid Mechanics (NNFM) 48.[50] V. John, Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder,Int. J. Numer. Meth. Fl. 44 (2004) 777 – 788.