An accurate and efficient method for the incompressible Navier–Stokes equations using the projection method as a preconditioner
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 without imposing the constraint of incompressibility, yielding an “intermediate” fluid velocity field , and generally . 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 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 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 and [18], [23], and obtaining high-order accuracy for 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 to , for a variety of non-trivial physical boundary conditions. The accuracy of the scheme is also demonstrated at 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 , 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 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 equationswhere are fixed physical coordinates, is the fluid velocity, is the pressure, 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 Cartesian grid with uniform meshwidths . The discretized velocity field is defined in terms of those
Linear solvers
Solving for and in Eqs. (17), (18) requires the solution of the block linear systemwhere and . Note that the choice of sign for and 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 , and 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)
A projection method for locally refined grids
J. Comput. Phys.
(1996)- et al.
A cell-centered adaptive projection method for the incompressible Euler equations
J. Comput. Phys.
(2000) Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries
J. Comput. Phys.
(2003)- et al.
A fourth-order auxiliary variable projection method for zero-Mach number gas dynamics
J. Comput. Phys.
(2008) - et al.
Application of a fractional-step method to incompressible Navier–Stokes equations
J. Comput. Phys.
(1985) - et al.
A second-order projection method for the incompressible Navier–Stokes equations
J. Comput. Phys.
(1989) A computational model of the cochlea using the immersed boundary method
J. Comput. Phys.
(1992)A 2nd-order projection method for the incompressible Navier–Stokes equations in arbitrary domains
J. Comput. Phys.
(1994)- et al.
A conservative adaptive projection method for the variable density incompressible Navier–Stokes equations
J. Comput. Phys.
(1998) - et al.
An adaptive version of the immersed boundary method
J. Comput. Phys.
(1999)