An efficient method for the incompressible Navier-Stokes equations on irregular domains with no-slip boundary conditions, high order up to the boundary.

An efficient method for the incompressible Navier-Stokes equations on irregular domains with no-slip boundary conditions, high order up to the boundary.


Common efficient schemes for the incompressible Navier-Stokes equations, such as projection or fractional step methods, have limited temporal accuracy as a result of matrix splitting errors, or introduce errors near the domain boundaries (which destroy uniform convergence to the solution). In this paper we recast the incompressible (constant density) Navier-Stokes equations (with the velocity prescribed at the boundary) as an equivalent system, for the primary variables velocity and pressure. We do this in the usual way away from the boundaries, by replacing the incompressibility condition on the velocity by a Poisson equation for the pressure. The key difference from the usual approaches occurs at the boundaries, where we use boundary conditions that unequivocally allow the pressure to be recovered from knowledge of the velocity at any fixed time. This avoids the common difficulty of an, apparently, over-determined Poisson problem. Since in this alternative formulation the pressure can be accurately and efficiently recovered from the velocity, the recast equations are ideal for numerical marching methods. The new system can be discretized using a variety of methods, including semi-implicit treatments of viscosity, and in principle to any desired order of accuracy. In this work we illustrate the approach with a 2-D second order finite difference scheme on a Cartesian grid, and devise an algorithm to solve the equations on domains with curved (non-conforming) boundaries, including a case with a non-trivial topology (a circular obstruction inside the domain). This algorithm achieves second order accuracy in the norm, for both the velocity and the pressure. The scheme has a natural extension to 3-D.


A critical issue in the numerical solution of the incompressible Navier-Stokes equations is the question of how to implement the incompressibility constraint. Equivalently, how to recover the pressure from the flow velocity, given the fact that the equations do not provide any boundary condition for the pressure. This has been an area of intense research, ever since the pioneering MAC scheme [18] of Harlow and Welch in 1965. Of course, one can avoid the problem by simultaneously discretizing the momentum and the divergence free equations, as in the difference scheme proposed by Krzywicki and Ladyzhenskaya [29], which can be shown to converge — while avoiding the need for any pressure boundary conditions. Approaches such as these, however, do not lead to efficient schemes.

Generally the dilemma has been that of a trade-off between efficiency, and accuracy of the computed solution near the boundary. However, many applications require both efficiency, and accuracy. For example, to calculate fluid solid interactions, both the pressure and gradients of the velocity are needed at the solid walls, as they appear in the components of the stress tensor. Furthermore, these objectives must be achievable for “arbitrary” geometries, not just simple ones with symmetries that can be exploited. Unfortunately, these requirements are not something that current algorithms are generally well suited for, as the brief review below is intended to show.1 However, we believe that algorithms based on a Presure Poisson Equation (PPE) reformulation of the Navier-Stokes equations — reviewed towards the end of this introduction — offer a path out of the dilemma. The work presented in this paper is, we hope, a contribution in this direction.

Projection methods

are very popular in practice because they are efficient. They achieve this efficiency by: (i) Interpreting the pressure as effecting a projection of the flow velocity evolution into the set of incompressible fields. That is, write the equations in the form , where is the appropriate projection operator, is the kinematic viscosity, is the flow velocity vector, and is the vector of applied body forces. (ii) Directly evolving the flow velocity. The question is then how to compute .

In their original formulation by Chorin [7] and Temam [48], the projection method was formulated as a time splitting scheme in which: First an intermediate velocity is computed, ignoring incompressibility. Second, this velocity is projected onto the space of incompressible vector fields — by solving a Poisson equation for pressure. Unfortunately this process introduces numerical boundary layers into the solution, which can be ameliorated (but not completely suppressed) for simple geometries — e.g. ones for which a staggered grid approach can be implemented [8].

The development of second order projection methods [2] provided greater control over the numerical boundary layers and accuracy in the pressure [4]. These are the most popular schemes used in practice. However, particularly for moderate or low Reynolds numbers, the effects of the numerical boundary layers can still be problematic [16]. Non-conforming boundaries add an extra layer of difficulty. The search for means to better control these numerical artifacts is an active area of research.

The numerical boundary layers in projection methods reflect in the known convergence results for them (e.g. [40]). Convergence is stated only in terms of integral norms, with the main difficulties near the boundary. There point-wise convergence (and even less convergence of the flow velocity gradient) cannot be guaranteed — even if the solution is known to be smooth. Hence the accurate calculation of wall stresses with these methods is problematic. Guermond, Minev and Shen [16] provide further details on convergence results, as well as an extensive review of projection methods and the improved pressure-correction schemes.

Two other methods for solving the Navier-Stokes equations are the immersed boundary [31], and the vortex-streamfunction [3] methods. These also decouple the calculation of the velocity and of the pressure. The immersed boundary method does so by introducing Dirac forces to replace the domain walls, which makes obtaining high order implementation of the boundary conditions difficult. The vortex-streamfunction formulation decouples the equations, but has dimensional limitations. An interesting variation of the vortex-streamfunction approach, using only local boundary conditions, is presented in reference [17].

Closely related to the immersed boundary methods are the penalty (alternatively: fictitious domain or domain embedding) methods — e.g. see [1]. These methods, effectively, replace solid walls in the fluid by a porous media with a small porosity . In the limit , this yields no slip and no flow-through at the solid walls. Two important advantages of this approach are that complicated domains are easy to implement, and that the total fluid-solid force can be computed using a volume integral, rather than an integral over the boundary of the solid. Unfortunately, the parameter introduces boundary layers which make convergence slow and high accuracy computations expensive, since cannot be selected independently of the numerical grid size.

Finally, we mention the algorithms based on a Pressure Poisson Equation (PPE) reformulation of the Navier-Stokes equations [19], which is the class of methods within which the work presented in this paper falls. In this approach the incompressibility constraint for the flow velocity is replaced by a Poisson2 equation for the pressure. This then allows an extra boundary condition — which must be selected so that, in fact, incompressibility is maintained by the resulting system. This strategy was first proposed by Gresho and Sani [15], who pointed out that adding as a boundary condition yields a system of equations that is equivalent to the Navier-Stokes equations. Unfortunately, their particular PPE formulation incorporates no explicit boundary condition that can be used to recover the pressure from the velocity, by solving a Poisson problem — for a more detailed discussion of this, see remark ? in this paper. In [19] the issue is resolved at the discrete numerical level, where they demonstrate high order schemes. For instance [19] demonstrates a fourth-order in space and 2nd-order in time implementation using overlapping grids. Subsequent work at the continuum level was later introduced by Henshaw et. al. [21] and Johnston and Liu [24]. Recently, work in PPE formulations have led to interesting improvements and analysis of projection methods [32]. In this paper — in equations (Equation 14Equation 15) — we present another PPE system, also equivalent to the Navier-Stokes equations, which allows an explicit recovery of the pressure given the flow velocity. A comparison with the one in [24] can be found in remark ? in this paper. Subtle issues can arise with semi-implicit approaches (pressure treated explicitly and viscous term implicitly) leading to time step restrictions [39]. In Section 5 we show that semi-implicit implementations of our scheme do not have time step restrictions of the diffusive type.

PPE reformulations of the Navier-Stokes equations, such as the one in equations (Equation 14Equation 15), or in reference [24], have important advantages over the standard form of the equations. First, the pressure is not implicitly coupled to the velocity through the momentum equation and incompressibility. Hence it can be directly (and efficiently) recovered from the velocity field by solving a Poisson equation. This allows the velocity field to be marched in time using the momentum equation, with the pressure interpreted as some (complicated) function of the velocity. Second, no spurious boundary layers are generated for neither the velocity, nor the pressure. This because:

  • There are no ambiguities as to which boundary conditions to use for the pressure — hence errors induced by not-quite-correct boundary conditions do not occur. Note that these errors should not be confused with the truncation errors that any discretization of the equations will produce. Truncation errors are controlled by the order of the approximation and, for smooth solutions, are uniformly small.

  • Incompressibility is enforced at all times.

Hence pressure and velocities that are accurate everywhere can be obtained, in particular: near the boundaries. Finally, PPE formulations allow, at least in principle, for the systematic generation of higher order approximations.

It follows that PPE strategies offer the promise of a resolution of the dilemma alluded to in the second paragraph of this introduction. They retain many of the advantages that have made projection methods popular, while not suffering from the presence of numerical boundary layers. On the negative side, the boundary conditions for PPE systems tend to be more complicated than the simple ones that the standard form of the equations has.

This paper is organized as follows. In § Section 2 we discuss the Pressure Poisson Equation (PPE) formulation of the Navier-Stokes equations. In § Section 3 we present a PPE formulation, using an alternative form of the boundary conditions, that allows a complete splitting of the momentum and pressure equations. Namely: the pressure can be recovered from the flow velocity without boundary condition ambiguities. In § Section 4 we address a numerical question (secular growth of the error under some conditions) and introduce a corrected formulation, involving a feedback parameter , where this potential problem is corrected. This system, which is equivalent to the Navier-Stokes equations, is displayed in equations (Equation 14Equation 15). A discussion and comparison with the alternative PPE formulation by Johnston and Liu [24] is also included in this section. In § Section 5 we introduce and examine a second order in time, semi-implicit in viscosity scheme. We show the scheme is unconditionally stable with no diffusive time step restriction. In § Section 6 we describe a second order solver for the new system introduced in § Section 4, for irregular domains embedded within a cartesian staggered grid. In § Section 7, we implement and test our proposed schemes. This section includes convergence plots, indicating second order uniform convergence all the way to the boundary ( norm), of the pressure, the flow velocity, and of the derivatives of the flow velocity. Finally, § Section 8 has the conclusions. The contents of the two appendices is as follows. In Appendix A we present an extended system of equations, valid for arbitrary flows, which for smooth solutions has the Navier-Stokes equations as an attractor. In Appendix B, for completeness, we display formulas for the boundary condition for general conforming boundaries in curvilinear coordinates.

2The Pressure Poisson Equation.

In this section we introduce the well-known pressure Poisson equation (PPE), and use it to construct a system of equations (and boundary conditions) equivalent to the constant-density (hence incompressible) Navier-Stokes equations, with the velocity prescribed at the boundaries. Specifically, consider the incompressible Navier-Stokes equations in a connected domain , where or , with a piece-wise smooth boundary . Inside , the flow velocity field satisfies the equations

where is the kinematic viscosity,3 is the pressure, are the body forces, is the gradient, and is the Laplacian. Equation (Equation 1) follows from the conservation of momentum, while ( ?) is the incompressibility condition (conservation of mass).

In addition, the following boundary conditions apply


is the outward unit normal on the boundary, and is the area (length in 2D) element on . Equation (Equation 3) is the consistency condition for , since an incompressible fluid must have zero net flux through the boundary.

Finally, we assume that initial conditions are given

Next we introduce the Pressure Poisson Equation (PPE). To obtain the pressure equation, take the divergence of the momentum equation (Equation 1), and apply equation ( ?) to eliminate the viscous term and the term with a time derivative. This yields the following Poisson equation for the pressure

Two crucial questions are now (for simplicity, assume solutions that are smooth all the way up to the boundary)

  • Can this equation be used to replace the incompressibility condition ( ?)?

    Since — given (Equation 5) — the divergence of (Equation 1) yields the heat equation

    for and , it would seem that the answer to this question is yes — provided that the initial conditions are incompressible (i.e.: for ). However, this works only if we can guarantee that, at all times,

    for .

  • Given the flow velocity , can (Equation 5) be used to obtain the pressure ?

    Again, at first sight, the answer to this question appears to be yes. After all, (Equation 5) is a Poisson equation for , which should determine it uniquely — given appropriate boundary conditions. The problem is: what boundary conditions? Evaluation of (Equation 1) at the boundary, with use of (Equation 2), shows that the flow velocity determines the whole gradient of the pressure at the boundary, which is too much for (Equation 5). Further, if only a portion of these boundary conditions are enforced when solving (Equation 5) — say, the normal component of (Equation 1) at the boundary, then how can one be sure that the whole of (Equation 1) applies at the boundary?

The issue in item Section 2a can be resolved easily, and we do so next. We postpone dealing with the issue in item Section 2b till the next section, § Section 3. Since the addition of an equation for the pressure allows the introduction of one extra boundary condition, we propose to replace the system in (Equation 1), ( ?), and (Equation 2) by the following Pressure Poisson Equation (PPE) formulation:

for , with the boundary conditions

for  — where, of course, the restriction in (Equation 3) still applies. The extra boundary condition is precisely what is needed to ensure that the pressure enforces incompressibility throughout the flow (see item Section 2a)

3Theoretical Reformulation

We still need to deal with the issue raised in item Section 2b. In particular, in order to implement the ideas in remark ?, we need to split the boundary conditions for (Equation 8?), in such a way that: (i) there is a specific part of the boundary conditions that is used with (Equation 8) to advance the velocity field in time, given the pressure. (ii) The remainder of the boundary conditions is used with ( ?) to solve for the pressure — at each fixed time, given the flow velocity field . This is the objective of this section.

The conventional approach in projection or fractional step methods is to associate (Equation 9) with equation (Equation 8). This is reasonable for evolving the heat-like equation (Equation 8). Unfortunately, it has the drawback of only implicitly defining the boundary conditions for the Poisson equation ( ?). Specifically, the correct pressure boundary conditions are those that guarantee for , and such boundary conditions cannot easily be known a priori as a function of . Hence one is left with a situation where the appropriate boundary conditions for the pressure are not known. This leads to errors in the pressure, and in the incompressibility condition, which are difficult to control. In particular, errors are often most pronounced near the boundary where no local error estimates can be produced, even for smooth solutions. By the latter we mean that the resulting schemes cannot be shown to be consistent, all the way up to the boundary, in the classical sense of finite differences introduced by Lax [30].

In this paper we take a different approach, which has a similar spirit to the one used by Johnston and Liu [24] — see remark ?. Rather than associate all the components of (Equation 9) with equation (Equation 8), we enforce the tangential components only, and complete the set of boundary conditions for (Equation 8) with ( ?). Hence, when evolving equation (Equation 8) we do not specify the normal velocity on the boundary, but — through the divergence condition ( ?) — specify the normal derivative of the normal velocity. Finally, the (as yet unused) boundary condition on the normal velocity, for , is employed to obtain an explicit boundary condition for equation ( ?). We do this by requiring that the pressure boundary condition be equivalent to for  — which then guarantees that the normal component of (Equation 9) holds, as long as the initial conditions satisfy it. This objective is easily achieved: dotting equation (Equation 8) through with , and evaluating at the boundary yields the desired condition. The equations, with their appropriate boundary conditions, are thus:


Again, for smooth (up to the boundary) enough solutions of the equations: the incompressible Navier-Stokes equations (Equation 1?), with boundary conditions as in (Equation 2), are equivalent to the system of equations and boundary conditions in (Equation 10Equation 11).

For the sake of completeness, we display now the calculation showing that the boundary condition splitting in (Equation 10Equation 11) recovers the normal velocity boundary condition . To start, dot the first equation in (Equation 10) with the normal , and evaluate at the boundary. This yields

Next, eliminate from this last equation — by using the boundary condition for the pressure in (Equation 11), to obtain

This is a trivial ODE for the normal component of the velocity at each point in the boundary. Thus, provided that initially, it holds for all time.

An important final point to check is the solvability condition for the pressure problem. Given the flow velocity at time , equation (Equation 11) is a Poisson problem with Neumann boundary conditions for the pressure. This problem has a solution (unique up to an additive constant) if and only if the “flux equals source” criteria

applies. This is satisfied because:

  • ,

  • ,

  • ,

where we have used Gauss’ theorem, incompressibility, and (Equation 3).

4Modification for Stability

For the system of equations in (Equation 10Equation 11), it is important to notice that

  • The tangential boundary condition on the flow velocity, for , is enforced explicitly.

  • The incompressibility condition, for , is enforced “exponentially”. By this we mean that any errors in satisfying the incompressibility condition are rapidly damped, because satisfies (Equation 6Equation 7). Thus this condition is enforced in a robust way, and we do not expect it to cause any trouble for “reasonable” numerical discretizations of the equations.

  • By contrast, the normal boundary condition on the flow velocity, for , is enforced in a rather weak fashion. By this we mean that errors in satisfying this condition are not damped at all by equation (Equation 12). Thus, this condition lacks the inherent stability provided by the heat equation. In practice, numerical errors add (effectively) noise to equation (Equation 12), resulting in a drift of the normal velocity component. This can have de-stabilizing effects on the behavior of a numerical scheme. Hence it is a problem that must be corrected.

In this section we alter the PPE equations (Equation 10Equation 11) to address the problem pointed out in item Section 4c. We do this by adding an appropriate “stabilizing” term. The idea is similar in nature to the feedback controller introduced in [12] for immersed boundary methods. Our goal here is to develop a pair of differential equations — fully equivalent to (Equation 10Equation 11) — which are suitable for numerical implementation. In order to resolve the issue in item Section 4c, we add a feedback term to the equations, by altering the pressure boundary condition. Specifically, we modify the equations from (Equation 10Equation 11) to


where is a numerical parameter — see § Section 4.1. This system is still equivalent to the incompressible Navier-Stokes equations and boundary conditions in (Equation 1Equation 2), for smooth (up to the boundary) enough solutions . The heat equation (Equation 6Equation 7) for still applies, while the equation for the evolution of the normal velocity at the boundary changes from (Equation 12) to:

Thus, if initially, it remains so for all times. In addition, this last equation shows that this new system resolves the issue pointed out in item Section 4c.

Finally, we check what happens to the solvability condition for the pressure problem, given the system change above. Clearly, all we need to do is to modify equation (Equation 13) by adding — to its left hand side, the term

where we have used incompressibility, and equation (Equation 3). It follows that solvability remains valid.

4.1Selection of the parameter .

For numerical purposes, here we address the issue of how large should be, by using a simple model for the flow’s normal velocity drift. Notice that no precision is needed for this calculation, just order of magnitude. In actual practice, one can monitor how well the normal velocity satisfies the boundary condition, and increase if needed. In principle one should be wary of using large values for , since this will yield stiff behavior in time. However, the calculation below shows that does not need to be very large, and does not depend on the grid size .

It seems reasonable to assume that one can model how the numerical errors affect the ODE (Equation 16) for , by perturbing the coefficients of the equation, and adding a forcing term to it. Hence we modify equation (Equation 16) as follows

where characterizes the size of the errors (determined by the order of the numerical method), while and are functions encoding the numerical errors. What exactly they are depends on the details of the numerical discretization, but for this calculation we do not need to know these details. All we need is that4

where and are some positive constants, with , and  — but not necessarily close to one.

The solution to (Equation 18) is given by

where is the initial value, , and . The crucial term is , since the first term decreases in size, and starts at the initial value. However

Within the framework of a numerical scheme, the normal boundary velocity may deviate from the prescribed velocity by some acceptable error . Thus we require or less, which — given (Equation 21) — will be satisfied if

For the second order numerical scheme in § Section 6, it is reasonable to expect that , and to require that . Then (Equation 22) reduces to . Of course, we do not know (a priori) what is; this is something that we need to find by numerical experimentation — see the first paragraph in this § Section 4.1. For the numerical calculations reported in § Section 7, we found that values in the range gave good results.

For this scheme, comparing the time step restriction imposed by feedback to the one imposed by diffusion , results in a stiffness ratio that scales as

Hence for low to moderate Reynolds numbers, can be choosen quite large without effecting the stiffness.

5Stability of Semi-implicit Schemes

In the case of moderate to low Reynolds number flows, the stiff viscosity term requires very small time steps when treated explicitly. Hence, there is a large practical interest in treating implicitly while keeping the associated pressure explicit. In this section we write down and analyze semi-implicit schemes for the PPE splitting (Equation 10Equation 11). Inspired by many of the ideas from [24], we then show the subsequent schemes are unconditionally stable.

To analyze the stability of the semi-implicit schemes, we require several properties of the Hodge-Helmholtz decomposition. Given , then has a unique orthogonal decomposition as

where is determined by with boundary conditions . Hence is divergence free with zero normal boundary component. Moreover the component is orthogonal in the norm to every gradient field (i.e., is the standard inner product on ). We may therefore write and where and are complementary orthogonal projections (ie. , and ). As a result, the following identities hold for any vector fields and


To motivate the stability proofs for semi-implicit time discretizations, note that the linear Navier-Stokes equations can also be written in an equivalent projection form

Dotting equation (Equation 23) through by5 and integrating yields

Here the first line follows directly from integration by parts, noting that the boundary integrals vanish. The last equation can be written in a more compact form:

From this it follows that the curl and divergence are bounded in the sense. As a first example of a semi-implicit scheme, we analyze the backward Euler discretization for the Stokes equation:


We note that the scheme (Equation 25)–(Equation 26) preserves the divergence free condition on the velocity exactly. Assume for . Taking the divergence of (Equation 25), and setting yields

Since does not have a negative eigenvalue, for all is the only solution to (Equation 27). Consequently, the pressure equation for automatically satisfies the Neumann consistency condition and is therefore well defined. As a direct consequence of preserving the divergence condition, the pressure as defined by (Equation 26) is equivalent to the projection:

A stability proof for equations (Equation 25)–(Equation 26) closely follows the one for periodic channel flow in [24]. Specifically, dot both sides of equation (Equation 25) by , and integrate (this is the discrete analog of the steps required to derive equation (Equation 24)). The first two terms on the left hand side of equation (Equation 25) become

In addition, we have


For brevity, we also introduce the energy where

Finally, combining everything, we have

so that

For simply connected domains, the energy bound (Equation 29) implies the scheme (Equation 25)–(Equation 26) is unconditionally stable. Unfortunately the bound (Equation 29) is not sufficient to prove stability in periodic domains, or domains with holes. The reason here is that for any vector field that has in , and in . In simply connected domains, the only such vector fields are . On the other hand, in periodic domains or domains with holes, there exist vector fields , which satisfy and on and . Here the energy bound (Equation 29) does not control the growth of these modes. We note however that such vector fields have non zero flux at the boundary. Hence, we expect that the modified scheme (Equation 14)–(Equation 14) with addition of the term in the boundary condition for the pressure will stabilize the growth of such modes.

5.1Stability of a second order scheme

Following a procedure analogous to the one in [24], we discuss the stability of a second order Crank-Nicholson scheme where we treat the pressure with a second order Adams-Bashforth extrapolation:

where is given by (Equation 26). Proceeding with a normal mode analysis6, we set where satisfies the boundary conditions in (Equation 30), and is an eigenvalue of the time stepping operator. Upon substitution we obtain

In analogy with the steps in the previous section, we dot equation (Equation 31) through by and integrate by parts to obtain

Equation (Equation 32) is now a quadratic of the form for the eigenvalues where

Since the coefficients satisfy the following inequalities and , it follows [46] that the eigenvalues lie within the unit circle. The preceding argument indicates that the normal modes remain bounded for the second order Crank-Nicholson scheme (Equation 30). Hence the scheme is stable for domains with simple geometries7.

6A First-Order Explicit Scheme on Irregular Domains

In this section we outline an explicit numerical scheme for solving the coupled differential equations (Equation 14Equation 15) on a two-dimensional irregular domain. To achieve such a scheme, we decouple the pressure and velocity fields, and explicitly treat each term in the time evolution of (Equation 14Equation 15). Specifically, since both the right hand side and boundary conditions of equation (Equation 15) depend solely on , we may view the pressure as a computable functional of the velocity, . The computation of requires the solution of a Poisson equation with Neumann boundary conditions. With this in mind, the momentum equation then has the form , where has a complicated, yet numerically computable form:

We now use an explicit forward Euler scheme to discretize the time evolution for (Equation 14Equation 15), paired with an appropriate discretization in space described later in this section. This yields the scheme

with boundary conditions and for , where the pressure is given by

with the boundary condition

for . Here, starting with the initial data , a superscript is used to denote a variable at time , where is the time step.

6.1Space grid and discretization.

To discretize the equations in space, we use finite differences over a cartesian, square (), staggered grid. The pressure values are stored at the nodes of the grid, while the horizontal and vertical components of the velocity are stored at the mid-points of the edges connecting the grid nodes (horizontal component on the horizontal edges and vertical component on the vertical edges).

When handling an arbitrary curved boundary, we do not conform the boundary to the grid, but rather we immerse it within the regular mesh — see Figure 1. Then, to numerically describe the domain boundary, we identify a set of points in , say for  — see item Section 6c below. These points are located at distances apart, so that the resolution of the boundary is comparable with that of the numerical grid.

Figure 1: This plot shows the staggered grid, and the boundary. The numerical pressure values correspond to the graph nodes, while the velocities (arrows) correspond to the edge midpoints. Here the circles (\circ\/) and squares (\Box\/) denote ghost pressure points, and boundary velocities, respectively. These are used to implement the boundary conditions in the Poisson and momentum equations, respectively. The diamonds (\Diamond\/) denote the points (x_b\/,\,y_b)_i\/, used to represent the boundary \partial\/\Omega\/.
Figure 1: This plot shows the staggered grid, and the boundary. The numerical pressure values correspond to the graph nodes, while the velocities (arrows) correspond to the edge midpoints. Here the circles () and squares () denote ghost pressure points, and boundary velocities, respectively. These are used to implement the boundary conditions in the Poisson and momentum equations, respectively. The diamonds () denote the points , used to represent the boundary .

For the heat equation , the stability restriction for the standard scheme using a 5 point centered differences approximation for the Laplacian, and forward Euler in time, is

where and is the space dimension. Since the method described here uses exactly the same approach to advance the velocity flow field , we expect the same restriction (with, perhaps, a different constant ) to apply. For the 2D numerical calculations presented in § Section 7, we found that the algorithm was stable with , while generally produced unstable behavior. Below we separately address the numerical implementation of equations (Equation 34) and (Equation 35).

6.2Poisson Equation

To solve the pressure Poisson equation (Equation 35) with Neumann boundary conditions, we use the ghost point idea with a method by Greenspan [14] to implement the boundary conditions. As discussed above, we embed the domain within a cartesian square grid, and classify the computational points in the grid into inner and ghost points (see next paragraph). Then: (i) for the inner points we discretize the Laplace operator using the standard 5 point centered differences stencil. (ii) For the ghost points we obtain equations from an appropriate discretization of the Neumann boundary condition. In particular, if there are inner points, and ghost points, then the discretized pressure is represented by a vector where .

The computational domain, , for the pressure is comprised by the inner points and the ghost points, defined as follows:

  • The inner points are the points in the cartesian grid located inside .

  • The ghost points are the points in the cartesian grid that are either outside or on , and which are needed to complete the 5 point stencil for some inner point.

  • Construct the set used to track the boundary as follows: For every pressure ghost point in the grid: , , select the closest point in the boundary . We will use these points to produce equations that approximate the Neumann boundary conditions (one equation per point). Thus, with this choice for , .

It follows that, on the computational domain , the Poisson equation (including the boundary conditions) is discretized by equations inside , and equations derived from the boundary conditions. Specifically

  • We use the standard 5 point centered differences stencil for the Laplacian, to obtain one equation for each of the inner points. Similarly, we use a second order approximation for the forcing terms on the right hand side of the equation.

  • We construct one boundary condition equation for each pressure ghost point , . This is done by building an extrapolation operator which acts on 6 nearby points from the computational domain to obtain a second order approximation to the normal derivative at the corresponding point  — see item Section 6c above. The 6 points are comprised of the 5 point stencil used to approximate the Laplacian, plus one of the four next closest points to the center point. The extrapolation operator, which contains 6 coefficients, is then build to linearly interpolate the normal derivative using local values of the pressure. Hence each boundary point provides one equation coupling six values of the pressure.

    We use a similar extrapolation operator to approximate the terms on the right hand side of the Neumann boundary conditions.

    Hence, each boundary condition equation involves the ghost point, points inside the domain, and (possibly) other nearby ghost points.

Thus the Poisson equation can be written in a matrix form with the following structure

Here and are the discrete matrix representations for the derivatives () in the Neumann boundary condition, and Laplacian (), respectively. The terms and are the discrete representations of the source terms and applied boundary conditions for the Poisson equation:

Thus the pair of equations

are the discrete analog of equation (Equation 15).

6.3Momentum Equation

The main numerical difficulty with implementing (Equation 34) is produced by the boundary conditions. Specifically: on a cartesian staggered grid the boundary conditions and couple the “horizontal”, , and “vertical”, , components of the flow velocity — with the exception of the special case of a boundary aligned with the grid, where there is no coupling. Hence, in general, the implementation of the boundary conditions requires the solution (at every time step) of a linear system of equations that couples all the boundary velocities (these are defined below). It may be possible to avoid the difficulties with the boundary conditions discussed in this section by using a different spatial discretization. For example, overlapping grids which conform to the domain boundaries, such as those found in [19], may avoid the coupling between velocity boundary components.

The computational domain for the velocities, , is defined in terms of the edges in the cartesian grid with which the velocities are associated (see Figure 1). To define , it is convenient to first introduce the extended set of pressure nodes, , from which is easily constructed:

  • A pressure node in the cartesian grid belongs to if and only if it either belongs to , or if it is connected (by an edge in the grid) to a ghost pressure node that lies on .

  • A velocity is in if and only if: (a) Its corresponding edge connects two points in . (b) At least one of the two points is in (or ).

The computational domain is, in turn, sub-divided into inner and boundary velocity edges

  • An edge in is an inner velocity if and only if includes the four other edges needed to compute the Laplacian (either or ) at the edge mid-point, using the 5 point centered differences approximation.

  • An edge in is a boundary velocity if and only if it is not an inner velocity. Let be the number of boundary velocities.

The solution of the momentum equation (Equation 14) is thus performed as follows:

  • At the start of each time step the right hand side in (Equation 34) is approximated (to second order, using centered differences) at the inner velocity locations, which can then be updated to their values at the next time.

  • Next, using the boundary conditions, the values of the boundary velocities at the next time step are constructed from the inner velocities. This “extension” process is explained below.

Let be the vector of boundary velocities. Then is determined by two sets of equations, corresponding to the discretization on of and . In our approach the divergence free criteria is enforced “point-wise”, while the condition on the tangential velocity is imposed in a least squares sense.

Implementation of the divergence free boundary condition. At first sight, this boundary condition appears to be the hardest to implement, since it is a Neumann condition (essentially) prescribing the value of the normal derivative of the normal component of the velocity, in terms of the tangential derivative of the tangential velocity. However, because (for the exact solution) everywhere, an implementation of this condition which is second order consistent (in the classical sense of finite differences introduced by Lax [30]) is easy to obtain, as follows:

First: identify the pressure nodes in , which have at least one boundary velocity as an adjacent edge. These pressure nodes lie either on , or inside , and are all within a distance of  — definitely no further away than . Second: for each of these points, use centered differences to approximate the flow divergence at the point, and set . This provides equations that couple the (unknown) boundary velocities to the (known) inner velocities. In matrix form this can be written as

where is the portion of the discrete divergence operator acting on the unknown boundary velocities, while is the associated flux derived from the known inner velocities. is a rectangular, very sparse, matrix whose entries are and  — note that can be eliminated from these equations.

Implementation of the tangential velocity boundary condition. It is easy to see that , since most of the pressure nodes in the selected set will connect with two boundary boundary velocities. Thus (Equation 38) above provides approximately one half the number of equations needed to recover from the inner velocities. It would thus seem natural to seek for additional equations using the other boundary condition. Namely, find points on , and at each one of them write an approximation to using nearby points in . Unfortunately, this does not work. It is very hard to do the needed approximations in a fashion that is robust relative to the way is embedded in the rectangular grid. Our attempts at this simple approach almost always lead to situations where somewhere along an instability was triggered.

To avoid the problem stated in the previous paragraph, we over-determine the implementation of the tangential velocity boundary condition, and solve the resulting system in the least squares sense. The boundary condition is replaced by the minimization of a (discrete version) of a functional of the form

where is some (strictly positive) weight function. This approach yields a robust, numerically stable, approximation — fairly insensitive to the particular details of how is embedded within the cartesian grid.

Figure 2: This plot illustrates the implementation of the momentum equation boundary conditions. The circles (\circ\/) indicate the points at which the velocity divergence is set to zero. The six arrows indicate the horizontal velocity components in the patch, \mathcal{P}_{j\,\nu}^u\/, used to extrapolate u\/ to the three boundary points indicated by the diamonds (\Diamond\/). The squares (\Box\/) denote the boundary velocities — which are part of the patch \mathcal{P}_{j\,\nu}^u\/.
Figure 2: This plot illustrates the implementation of the momentum equation boundary conditions. The circles () indicate the points at which the velocity divergence is set to zero. The six arrows indicate the horizontal velocity components in the patch, , used to extrapolate to the three boundary points indicated by the diamonds (). The squares () denote the boundary velocities — which are part of the patch .

In fact, we do not implement the minimization of a functional as in (Equation 39), but a simpler process which is (essentially) equivalent. To be specific, we start with the set of points in  — see item Section 6c. For each , where , we identify a local horizontal, , and vertical, , velocity “patch”. These patches — see Figure 2 — have the following properties

  • Each patch contains both inner and boundary velocities.

  • Each patch contains 6 velocities, in an appropriate structure, so that these velocities can be used to extrapolate (with second order accuracy) values of the corresponding velocity ( or ) to nearby points along the boundary. For example: the 5 point stencil used to approximate the Laplacian, plus one of the four next closest points to the center point.

  • The union of all the patches contains all the boundary velocities.

These patches are used to (linearly) interpolate/extrapolate the velocities to nearby points along the boundary. For example, the horizontal velocity at point is extrapolated using patch by

where is the velocity at , and are precomputed coefficients which linearly extrapolate to second order.

Each patch is used to extrapolate the velocity to the three “closest” boundary points to the patch. In this fashion, for every , , three different approximations to the tangential velocity at the boundary follow. These involve different (linear) combinations of velocities at the nearby edges in , including both boundary and inner velocities. In this fashion, a set of additional linear equations — beyond those in equation (Equation 38) — for the boundary velocities follow.

Putting everything together, an overdetermined system of equations for the boundary velocity vector follows, which can be written in the form

Here , and the second equation is to be solved in the least squares sense, subject to the constraint imposed by the first equation. One way to do this is as follows: write

where is a matrix whose columns are a basis for the kernel of  — hence , is a particular solution of (Equation 40), and is a constant vector parameterizing the space of (numerical) divergence free boundary velocities. Substituting the ansatz (Equation 41) into ( ?) yields

which is a constraint-free least squares problem for the vector . Thus we can write

which, together with (Equation 41), gives the solution .

6.4Comparison with the Projection Method

The scheme proposed here — which involves the implementation of the equations in (Equation 34Equation 35), appears to be quite similar to fractional step methods, such as Chorin’s [7] original projection method. Specifically: advancing one time step requires both the evolution of a diffusion equation for the flow velocity, and the inversion of a Poisson equation for the pressure — which is the same as in projection methods. But there are important differences, mainly related to the implementation of the boundary conditions, and of the incompressibility condition. Below we highlight the similarities and differences between the pressure Poisson approach proposed here and, for simplicity, the projection method in its original formulation [7]. A (fairly recent) thorough review of projection methods, exploring the improvements to the approach since [7], as well as their drawbacks, can be found in [16].

First, the starting point for the method here is the discretization of a reformulation of the equations: (Equation 14Equation 15), not the “original” Navier-Stokes equations (Equation 1?). In this equivalent set: (i) there is a natural way to recover the pressure from the flow field at any given time. (ii) The time evolution automatically enforces incompressibility. As a consequence, there is (in-principle) no limitation to the order in time to which the reformulated equations can be numerically discretized. In contrast, the projection method is equivalent to an approximate matrix factorization [16] of the discrete differential operators coupling and . This approximate factorization yields a splitting error in time, which is very hard to circumvent in order to achieve higher accuracy.

Second, the pressure Poisson formulation used here ensures both that, at every point in the time integration: (i) the normal and tangential velocity boundary conditions are accurately satisfied. (ii) The zero divergence condition is accurately satisfied. On the other hand, in the original projection method, the step advancing the flow velocity forward in time enforces the velocity boundary conditions, but cannot guarantee that incompressibility is maintained. In the other step, the pressure is used to recover incompressibility — by projecting the flow velocity field onto the space of divergence free fields. Unfortunately, while removing the divergence, this second step does not necessarily preserve the correct velocity boundary values [16].

Finally the modified equation (Equation 34) implemented here differs from the projection’s method velocity forward step primarily by the boundary conditions imposed. Specifically, the projection method imposes Dirichlet boundary conditions for each component of the velocity field when propagating the flow velocity. Dirichlet boundary conditions have the obvious advantage of decoupling the velocity components. Numerically this means that an implicit treatment of the stiff viscosity operator is straightforward. By contrast, on a regular grid, the boundary conditions in (Equation 34) couple the boundary velocities. Although this coupling does not pose a serious difficulty for explicit schemes, it adds a (potentially serious) difficulty to the implementation of any scheme that treats the viscosity operator in the equations implicitly.


The objective of this section is to study the convergence and accuracy of the schemes proposed in § Section 5 and § Section 6, and to test the methods for a simple physical problem. In particular: to verify the theoretical prediction that the temporal splitting in § Section 5 is second order in time, while the spatial discretization in § Section 6 is second order in space, all the way up to the boundary.

To this end, here we present the results from numerical computations for (i) Externally forced flow on a square domain in § Section 7.1, (ii) Computation of flow in a lid driven cavity in § Section 7.2, and (iii) Externally forced flow on an irregular domain in § Section 7.3. In examples (i) and (iii), the external forcing is selected so that an exact (analytic) solution to the equations can be produced, following the same procedure as in [16].8 In what follows the numerical errors are measured using the discrete grid norm defined by

Here the maximum is over all the indexes and  — corresponding to the appropriate grid coordinates in the computational domain — for the field in the staggered grid. Namely: nodes for the pressure , mid-points of the horizontal edges for , and mid-points of the vertical edges for . Note that neither ghost pressure points, nor boundary velocities, are included within the index set in (Equation 44) — we consider these as auxiliary numerical variables, introduced for the purpose of implementing the boundary conditions, but not part of the actual solution. Finally, for the exact (analytic) solution, grid values are obtained by evaluation at the points . For example, to compute the error in the pressure, first define the discrete field  — where is the numerical pressure field and is the exact continuous field, and then compute the norm above for .

As should be clear from the prior sections, the main issue we aim to address in this paper, is that of how to effectively implement the incompressibility condition and the boundary conditions for the pressure, avoiding the difficulties that projection and fractional step methods have. These are problems that are not related to the nonlinearities in the Navier-Stokes equations, and occur even in the absence of the advective term . Hence, for simplicity, the calculations presented in subsections § Section 7.1 and § Section 7.3 do not contain the advective term. This is the same practice employed in [16]. Subsection § Section 7.2, however, investigates a physical flow and therefore includes the advective term.

7.1Flow on a square domain

Here we test the semi-implicit scheme outlined in § Section 5. Specifically, we present the results of solving the linear Navier-Stokes equations on the unit square , with no-slip and no flux boundary conditions (i.e.: on ), and viscosity . We select the forcing function to yield the following (incompressible) velocity and pressure fields:

We numerically evolve the equations with varying grid and time steps. We then compare the numerical results with the exact fields in (Equation 45?). For instance, Figure 3 shows the error fields for the velocity and pressure at time , for an grid, time step and feedback . Note that this test is essentially the same as the one used in [16]9, to test projection and fractional step methods. Hence the results in this subsection can be used as a comparison basis for the approach in this paper versus projection methods. Unlike the common situation with projection methods, in our case there are neither numerical boundary layers, nor numerical corner layers.

Figure 3: Error fields for the second order semi-implicit numerical scheme. The plots are at time t = 4\,\pi\/, for an 80\times 80\/ grid, with \Delta\,t=0.2\,(\Delta\,x)\/ and \lambda = 100. The horizontal velocity u\/ (top) and the pressure p\/ (bottom) error fields are shown. The error is uniform in size across the domain. There are neither numerical boundary layers, nor numerical corner layers.
Error fields for the second order semi-implicit numerical scheme. The plots are at time t = 4\,\pi\/, for an 80\times 80\/ grid, with \Delta\,t=0.2\,(\Delta\,x)\/ and \lambda = 100. The horizontal velocity u\/ (top) and the pressure p\/ (bottom) error fields are shown. The error is uniform in size across the domain. There are neither numerical boundary layers, nor numerical corner layers.
Figure 3: Error fields for the second order semi-implicit numerical scheme. The plots are at time , for an grid, with and . The horizontal velocity (top) and the pressure (bottom) error fields are shown. The error is uniform in size across the domain. There are neither numerical boundary layers, nor numerical corner layers.
Figure 4: Convergence plot for semi-implicit numerical scheme. The plot corresponds to a fixed 220\times 220 grid and varying time steps (\lambda = 30). The straight line corresponds to second order convergence, while the saturation occurs around \Delta t = \Delta x when the temporal errors become the same order as the spatial errors. The circles (\;\circ\/) and squares (\;\square\/) correspond to the L^\infty norm of the pressure and horizontal velocity respectively.
Figure 4: Convergence plot for semi-implicit numerical scheme. The plot corresponds to a fixed grid and varying time steps (). The straight line corresponds to second order convergence, while the saturation occurs around when the temporal errors become the same order as the spatial errors. The circles () and squares () correspond to the norm of the pressure and horizontal velocity respectively.

To verify the scheme is second order in time, we fix a grid and evolve the solution from to