An accurate and efficient method for the incompressible Navier–Stokes equations using the projection method as a preconditioner

https://doi.org/10.1016/j.jcp.2009.07.001Get rights and content

Abstract

The projection method is a widely used fractional-step algorithm for solving the incompressible Navier–Stokes equations. Despite numerous improvements to the methodology, however, imposing physical boundary conditions with projection-based fluid solvers remains difficult, and obtaining high-order accuracy may not be possible for some choices of boundary conditions. In this work, we present an unsplit, linearly-implicit discretization of the incompressible Navier–Stokes equations on a staggered grid along with an efficient solution method for the resulting system of linear equations. Since our scheme is not a fractional-step algorithm, it is straightforward to specify general physical boundary conditions accurately; however, this capability comes at the price of having to solve the time-dependent incompressible Stokes equations at each timestep. To solve this linear system efficiently, we employ a Krylov subspace method preconditioned by the projection method. In our implementation, the subdomain solvers required by the projection preconditioner employ the conjugate gradient method with geometric multigrid preconditioning. The accuracy of the scheme is demonstrated for several problems, including forced and unforced analytic test cases and lid-driven cavity flows. These tests consider a variety of physical boundary conditions with Reynolds numbers ranging from 1 to 30000. The effectiveness of the projection preconditioner is compared to an alternative preconditioning strategy based on an approximation to the Schur complement for the time-dependent incompressible Stokes operator. The projection method is found to be a more efficient preconditioner in most cases considered in the present work.

Introduction

Since its introduction by Chorin [1], [2], the projection method has been widely used as a solver for the incompressible Euler [3], [4], [5], [6], [7], [8] and Navier–Stokes [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29] equations. Generally speaking, the projection method is a fractional-step algorithm for incompressible flow problems which obtains updated values for the fluid velocity u and pressure p in two steps. First, an approximation to the momentum equation is solved over a time interval Δt without imposing the constraint of incompressibility, yielding an “intermediate” fluid velocity field u, and generally ·u0. To obtain an approximation to the updated fluid velocity which does satisfy the constraint of incompressibility, the projection method uses the Hodge decomposition to compute efficiently the projection of u onto the space of divergence-free vector fields. Doing so requires the solution of a Poisson equation. The solution φ of this linear equation is used both to project u to obtain the updated fluid velocity u and also to compute the updated pressure p.

The enduring popularity of the projection method may be attributed to its decomposition of a difficult problem into sub-problems for which efficient solvers are readily available. For instance, the highly influential second-order projection method of Bell et al. [10] requires only a solver for the diffusion equation arising from their implicit treatment of the viscous terms along with a Poisson solver for the discrete projection. Fast solvers employing FFT or multigrid algorithms may be used in both cases, thereby yielding incompressible Navier–Stokes solvers which are essentially algorithmically optimal.

Although the projection method framework yields timestepping schemes for the incompressible Navier–Stokes equations which require relatively simple linear solvers, this approach greatly complicates the specification of physical boundary conditions. Appropriate “artificial” boundary conditions must be prescribed for u and φ [18], [23], and obtaining high-order accuracy for u and p may not be possible in some cases, such as at outflow boundaries [30]. Moreover, there is relatively little theory guiding the choice of approximations used to impose these artificial boundary conditions, and different standard approximations to the same artificial boundary conditions may result in schemes which are unstable or stable but low-order accurate, whereas non-standard approximations can result in schemes which are stable and high-order accurate [21].

In the present work, we view the projection method as an approximate solver for the time-dependent incompressible Stokes equations, and we use this approximate solver as a preconditioner for the iterative solution of those equations. This iterative solver is used with a linearly-implicit staggered-grid discretization of the incompressible Navier–Stokes equations which employs an implicit treatment of the viscous terms and an explicit second-order Godunov (upwind) method for the nonlinear advection terms. The Godunov scheme used in the present work is based on xsPPM7 [31], a recent version of the piecewise parabolic method (PPM) [32].

By using an unsplit scheme instead of a fractional-step algorithm, we are able to prescribe general physical boundary conditions accurately and, in most cases, easily. (There appears to be some ambiguity in the specification of normal traction boundary conditions; we describe herein the approach which we have found to be the most accurate and robust.) Numerical results for forced and unforced analytic test cases presented in Section 5 demonstrate that our scheme attains fully second-order convergence rates over a broad range of Reynolds numbers, from Re=1 to Re3000, for a variety of non-trivial physical boundary conditions. The accuracy of the scheme is also demonstrated at Re30000 for a variety of boundary conditions, although these results appear under-resolved on all but the finest computational grids employed in the present work. When we employ inexact multigrid-preconditioned subdomain solvers, our numerical results demonstrate that the scheme is scalable in the sense that the number of linear solver and subdomain solver iterations are largely insensitive to the grid spacing. In many cases, the number of solver/sub-solver iterations is also demonstrated to be largely independent of Re.

The scheme is also demonstrated to converge to benchmark steady solutions [33], [34], [35] to the lid-driven cavity flow problem for Re=1000,5000, and 7500. Essentially first-order global convergence rates are obtained for the standard lid-driven cavity flow, which possesses well-known corner singularities [34], and fully second-order pointwise convergence rates are obtained for a regularized version of this problem [36]. We also demonstrate that improved accuracy (although not improved order of accuracy) may be obtained for the regularized lid-driven cavity problem by employing an alternative third-order boundary treatment. Because we do not use a fractional-step scheme, it is straightforward to use alternative, higher-order boundary condition implementations.

The projection method-based preconditioner described in the present work is similar to but distinct from a preconditioner based on an approximation to the Schur complement for the time-dependent incompressible Stokes operator. Such approximate Schur complement methods were originally introduced in the context of the Oseen equations [37], [38], and these methods have been thoroughly studied by Elman and co-workers [39], [40], [41]. (Although it seems likely that the projection preconditioner could be extended to treat the Oseen equations, this has not yet been done.) In the results presented in Section 5, we compare the efficiency of the projection preconditioner to a preconditioner based on an approximate Schur complement. In nearly all cases, the projection method is found to be a significantly more efficient preconditioner than the approximate Schur complement method. Although other preconditioning strategies for the incompressible Navier–Stokes equations have been described (see, e.g., [42], [43]), most alternative approaches do not appear to use subdomain solvers.

Although we describe a particular linearly-implicit timestepping scheme for the incompressible Navier–Stokes equations, we believe that our preconditioning approach is widely applicable. In particular, we expect that the basic approach described in the present work could be used to convert an existing fractional-step staggered-grid incompressible flow solver into an unsplit solver. Doing so would require only “wrapping” the existing projection code with a Krylov solver. The additional software required is modest and includes: (1) the implementation of code to apply the time-dependent incompressible Stokes operator to a given vector; (2) the implementation of a (possibly matrix-free) preconditioned Krylov solver for the incompressible Stokes system; and (3) the conversion of the existing staggered-grid projection solver into a preconditioner. Most of the work required by this conversion could be eliminated by using a high-quality numerical software library such as PETSc [44], [45], [46].

Finally, before proceeding, we note that although the present work treats the case of two spatial dimensions, the modifications required to extend the methodology to three spatial dimensions are straightforward. Results from an initial three-dimensional implementation of the scheme are encouraging and will be reported in future work.

Section snippets

The continuous equations of motion

Consider a fixed region ΩR2 which is filled with a viscous incompressible fluid. (To simplify the presentation, in this work the physical domain Ω is taken to be the unit square.) The equations of motion for the fluid are the incompressible Navier–Stokes equationsρut+(u·)u=-p+μ2u+f,·u=0,where x=(x,y)Ω are fixed physical coordinates, u(x,t)=(u(x,t),v(x,t)) is the fluid velocity, p(x,t) is the pressure, f(x,t)=(f1(x,t),f2(x,t)) is an applied body force, ρ is the uniform fluid density, and

Spatial discretization and finite difference approximations

In the present work, we employ a staggered-grid spatial discretization of the incompressible Navier–Stokes equations. (The layout of degrees of freedom used in the present work is the same as the MAC scheme [48], although note that we use spatial and temporal discretizations which are different from those of the MAC scheme.) Briefly, the physical domain Ω is described using a fixed N×N Cartesian grid with uniform meshwidths Δx=Δy=h=1N. The discretized velocity field is defined in terms of those

Linear solvers

Solving for un+1,k+1 and pn+12,k+1 in Eqs. (17), (18) requires the solution of the block linear systemρΔtI-μ2Lx0Gx0ρΔtI-μ2LyGy-Dx-Dy0un+1,k+1vn+1,k+1pn+12,k+1=ρΔtI+μ2Lxun-ρN1n+12,k+f1n+12ρΔtI+μ2Lyvn-ρN2n+12,k+f2n+120,where Nun+12,k=N1n+12,k,N2n+12,k and fn+12=f1n+12,f2n+12. Note that the choice of sign for Dx and Dy in Eq. (31) makes the matrix symmetric for certain choices of boundary conditions, e.g., in the case of periodic boundary conditions.

We solve Eq. (31) via the flexible GMRES

Numerical results

We now present numerical examples which demonstrate the accuracy of the discretization, and which compare the efficiency of the projection method-based preconditioner to the approximate Schur complement preconditioner. Errors and convergence rates are reported in discretized L1,L2, and L norms. When analytic solutions are available, errors and convergence rates are reported with respect to those exact solutions. When analytic solutions are not available, errors are estimated for the purpose of

Discussion

In this work, we have presented an unsplit, linearly-implicit discretization of the incompressible Navier–Stokes equations on a staggered grid along with an effective solution methodology for the resulting system of linear equations. Because we do not use a fractional-step algorithm, it is straightforward to specify physical boundary conditions accurately, and our convergence results demonstrate that the discretization is second-order accurate for a variety of choices of physical boundary

Acknowledgments

The author thanks D.B. Kothe of Oak Ridge National Laboratory for the initial suggestion to investigate the combination of Godunov methods and staggered-grid spatial discretizations, and acknowledges correspondence with W.J. Rider of Sandia National Laboratories on the use of Godunov methods with staggered-grid discretizations during the initial phase of this work. The author also thanks C.S. Peskin and D.M. McQueen of the Courant Institute-New York University for helpful comments on an early

References (62)

  • D.L. Brown et al.

    Accurate projection methods for the incompressible Navier–Stokes equations

    J. Comput. Phys.

    (2001)
  • Z.-L. Li et al.

    The immersed interface method for the Navier–Stokes equations with singular forces

    J. Comput. Phys.

    (2001)
  • R.D. Guy et al.

    Stability of approximate projection methods on cell-centered grids

    J. Comput. Phys.

    (2005)
  • B.E. Griffith et al.

    On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems

    J. Comput. Phys.

    (2005)
  • B. Yang et al.

    A second-order boundary-fitted projection method for free-surface flow computations

    J. Comput. Phys.

    (2006)
  • C. Min et al.

    A second order accurate projection method for the incompressible Navier–Stokes equations on non-graded adaptive grids

    J. Comput. Phys.

    (2006)
  • Z. Zheng et al.

    Runge–Kutta–Chebyshev projection method

    J. Comput. Phys.

    (2006)
  • B.E. Griffith et al.

    An adaptive, formally second order accurate version of the immersed boundary method

    J. Comput. Phys.

    (2007)
  • D.F. Martin et al.

    A cell-centered adaptive projection method for the incompressible Navier–Stokes equations in three dimensions

    J. Comput. Phys.

    (2008)
  • D.V. Le et al.

    An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains

    Comput. Methods Appl. Mech. Eng.

    (2008)
  • J.L. Guermond et al.

    An overview of projection methods for incompressible flows

    Comput. Methods Appl. Mech. Eng.

    (2006)
  • W.J. Rider et al.

    Accurate monotonicity- and extrema-preserving methods through adaptive nonlinear hybridizations

    J. Comput. Phys.

    (2007)
  • P. Colella et al.

    The piecewise parabolic method (PPM) for gas-dynamical simulations

    J. Comput. Phys.

    (1984)
  • U. Ghia et al.

    High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method

    J. Comput. Phys.

    (1982)
  • O. Botella et al.

    Benchmark spectral results on the lid-driven cavity flow

    Comput. Fluid

    (1998)
  • C.J. Roy et al.

    On the generation of exact solutions for evaluating numerical schemes and estimating discretization error

    J. Comput. Phys.

    (2009)
  • D. Silvester et al.

    Efficient preconditioning of the linearized Navier–Stokes equations for incompressible flow

    J. Comput. Appl. Math.

    (2001)
  • H.C. Elman et al.

    A parallel block multi-level preconditioner for the 3D incompressible Navier–Stokes equations

    J. Comput. Phys.

    (2003)
  • H. Elman et al.

    A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations

    J. Comput. Phys.

    (2008)
  • D.A. Knoll et al.

    On Newton–Krylov multigrid methods for the incompressible Navier–Stokes equations

    J. Comput. Phys.

    (2000)
  • Y.-F. Peng et al.

    Transition in a 2-D lid-driven cavity flow

    Comput. Fluid

    (2003)
  • Cited by (0)

    View full text