Analysis for Chorin''s Original Fully Discrete Projection Method and Regularizations in Space and Time.
Over twenty-five years ago, Chorin proposed a computationally efficient method for computing viscous incompressible flow which has influenced the development of efficient modern methods and inspired much analytical work. Using asymptotic error analysis techniques, it is now possible to describe precisely the kind of errors that are generated in the discrete solutions from this method and the order at which they occur. While the expected convergence rate is seen for velocity, the pressure accuracy is degraded by two effects: a numerical boundary layer due to the projection step and a global error due to the alternating or parasitic modes present in the discretization of the incompressibility condition. The error analysis of the projection step follows the work of E and Liu and the analysis of the alternating modes is due to the author. The two are combined to show the asymptotic character of the errors in the scheme. Regularization methods in space and time for recovering full accuracy for the computed pressure are discussed.
Introduction
In 1968, Chorin [3] proposed a computationally efficient method for computing viscous,
incompressible flow. The method was based on the primitive variables, velocity and
pressure, with all unknowns at the same grid points. The discretization was centered
in space (second order in space step h) and implicit in time (first order in time step
with the projection part of the Stokes operator split off for computational efficiency.
The discretization of the incompressibility condition allowed for alternating (parasitic or
checkerboard) modes. The idea of a projection step has been used in modern efficient
methods (e.g. [6]) (since the literature in this field is vast, the references here and in
what follows are only intended to be illustrative not exhaustive). Many authors have
considered the analysis of the projection step [4, 15, 16, 9] and have proposed higher
order corrections [13, 18, 14]. The precise description of the errors from the projection
step in [9] will be used in this work. The effect of parasitic modes on the accuracy of the
computed pressure has also been an area of much interest, especially in the finite element
method [8] and spectal method [2] communities. It is well known that the presence of
parasitic modes can lead to a degradation in the convergence rate of the pressure. A
precise description of this effect for Chorin's discretization is given in [20].
What the work in [9, 20] does is to characterize the errors from the projection step and
the parasitic modes precisely for smooth problems with a smooth discretization in space
and time. These analyses will be combined in this paper to give a precise description of
the errors from Chorin's original fully discrete scheme. Although there are no significant
new difficulties that arise from the interaction of the discretizations in space and time,
the author believes it is worthwhile to give the full error analysis for this historically
important algorithm. A simplified presentation of the boundary layer effects from the
projection step is given.
We describe the numerical order reducing effects on the pressure briefly below. If it
is assumed that convergence of second order in h might be
expected. While this is true for the velocity, it is not true for the pressure. In fact, the
discrete pressure from Chorin's method has an O(k 1=2 ) numerical boundary layer due to
the projection step and an O(h) global error due to the alternating modes. In order to
recover full accuracy for the pressure, we also consider regularization methods in space
and time. Computational evidence for all predictions is given. The reader may wish to
take a short tour through pressure errors from Chorin's original scheme
O(h) alternating terms dominate, to a space regularized scheme only the
boundary layer is left, to a fully regularized scheme (Fig.5) where the errors are
spatially smooth.
Some discussion should be made here about the real and artificial limitations of this
work. First of all, the analysis is presented for the two-dimensional (2D) Stokes equations
with homogeneous Dirichlet boundary conditions but can be extended in a straightforward
way to 3D with nonhomogeneous compatible boundary conditions and to smooth
solutions of the nonlinear Navier-Stokes equations. Secondly, the error expansions presented
assume a great number of compatibility conditions at When these are not
satisfied, convergence is not uniform in some quantities up to (see [12] for computational
results and formal analysis of the behaviour near
Thirdly, the computational issue of how to efficiently implement the projection step is
not addressed in this work. We use a simplified geometry for our numerical tests in which
it is easy to implement an exact projection efficiently which allows us to obtain refined
solutions to verify the predicted error structure. Fourthly, the temporal regularization
discussed in the final section uses unsplit time integration. In the context of split-step
projection methods, it would have been more appropriate to present a pressure increment
scheme [18, 6] but the error analysis of these schemes is not well understood in the fully
discrete case [10]. Finally, because the discrete divergence and gradient operators are
not adjoint, a simple stability result based on energy estimates as used in [4, 9] is not
possible for Chorin's original method when boundaries are present. Thus, the stability
analysis of the method is still an open problem.
In the next section, Chorin's original method is described. Then, computational
results showing the boundary layer and alternating errors in the pressure are presented.
In Section 4 we present the error analysis for the method, describing the alternating and
boundary layer errors and at what order they occur. Finally, in Section 5 we present
analysis and computation of regularized methods.
2 Description of the Scheme
We consider fluid flow in a simplified domain: a two-dimensional (2D) [0; 1] \Theta [0; 1]
channel with fixed walls on the top and bottom boundaries and periodic in the horizontal
direction. The incompressible Stokes equations are given below
where are the velocities, p is the pressure and - is the kinematic viscosity.
Boundary conditions are used. Initial data u 0 is given and it
is assumed that r \Delta u We note that p can only be determined up to an arbitrary
constant. A unique p is recovered if we require
Z
It is well known that any square integrable vector function can be orthogonally decomposed
into a divergence-free part with homogeneous normal boundary conditions and a
part that can be represented as the gradient of a scalar (see e.g. [5]). In this framework,
we can interpret the pressure gradient rp as a term that projects the right hand side of
(1) onto the space of divergence free fields and summarize its action with the projection
operator P:
To describe the discrete scheme we approximate in space on a regular grid with spacing
h and in time with spacing k. It is assumed that 1=h is even for convenience. We use
and P n
to denote approximations of u(ih; jh; nk) and p(ih; jh; nk) respectively. To
proceed, we need to define the approximate "projection", P h , derived by Chorin [3]. We
use discrete divergence D h and gradient operators G h based on long, centered differences,
i.e.
away from boundaries. Near the lower boundary we can use the fact that U
derive
On the boundary, using second order one sided differencing gives
Similar expressions apply on the upper wall. To divide an arbitrary vector W defined in
the interior of the domain into a gradient part G h P and a discrete divergence free part
U (D h \Delta must solve
and then
We summarize this process as
denotes the scalar corresponding to the gradient part of the vector.
This projection approach is convenient because it does not require the specification of
any additional "pressure boundary conditions". Such conditions can be considered to be
implicitly given by (7). We note that P h is not a projection matrix since D h and G h are
not negative adjoint. Also, the matrix D h \Delta G h has four null modes corresponding to the
four null modes of G h , constant vectors on the four subgrids shown in Figure 1. However,
(7) is solvable [1] up to the four null vectors of G h . The four arbitrary constants are
normalized using appropriate trapezoidal or midpoint approximations of (3). From the
structure of G h P it is easy to see that the errors on the four subgrids can be different,
leading to so-called alternating error expansions. The order that these effects enter the
velocity and pressure is described in detail below.
We now turn to a discretization in time. Chorin [3] proposed splitting the diffusion
step and the projection step in the following scheme:
where \Delta h denotes the usual five point approximation of the Laplacian with Dirichlet
data. This scheme gives an uncoupled system for P n+1 and W n+1 , an auxiliary quantity
computed during the diffusion sub-step. The fact that the system is decoupled is the
advantage of using the split-step technique.

Table

1: Normalized pointwise pressure errors e p and velocity errors e u (and estimated
convergence rates in h) from Chorin's scheme.
We note that in Chorin's original work an ADI technique was used to approximate
(10) and an iterative technique was used to approximate the projection step (11). Here
we analyze the underlying exact discretization for simplicity and because more modern
solution techniques are available that can efficiently solve these subproblems.
We consider computational results for this method below, showing the numerical
boundary layers from the projection step and the alternating errors from the parasitic
modes. The detailed analysis of these phenomena is then done in Section 4.
Computational Results
We demonstrate the types of errors discussed above with computational results for the
Stokes equations in the periodic channel. The initial data from [19] is used (a perturbation
of Poiseuille flow) with 1=64. Errors are calculated by comparing the solutions
from Chorin's method with those from the Marker and Cell (MAC) grid with high order
accurate explicit time stepping (the discrete pressures from this scheme have no alternating
or boundary layer effects [11]). Comparisons are made at
When relatively large and k relatively small), the pressure
errors are dominated by the O(h) alternating errors from the parasitic mode effects as
shown in Figure 2. Note that the error alternates in sign in the vertical direction only
and is not confined to a region near the boundary. If the computation had been done
in a box with vertical walls as well as horizontal walls there would also be horizontally
alternating components of the error.
When relatively small and k relatively large) the pressure
errors are dominated by the boundary layer due to the projection step with size and
width O(k 1=2 ). This is seen in the top picture of Figure 3. A contour plot of this same
data has jagged contour lines, showing the presence of (smaller) alternating terms. When
k is reduced to 0:01, the boundary layer is reduced in size and extent as shown in the
lower picture of Figure 3.
The error plots above verify the qualitative description of the errors. To show their
order reducing effects on computed P we perform computations with several
with
. In

Table

1 we see that P converges with first order (in h) and that U
converges with second order (in h).
We develop a asymptotic error expansion for the pressure and computed velocity for
Chorin's original scheme (10)-(12) consisting of regular and alternating errors and numerical
boundary layer terms as described in [20]. We will use the asymptotic descriptions
of P h from [20] and present a simplified derivation of the errors from the split-step
time-stepping first given in [9].
It is convenient to first derive an error expansion for W, the intermediate velocity,
and then derive expansions for U and P from (11) and (12). The update equation for
W is
with boundary conditions
convenient scaling for this analysis. We index the errors by powers of h so O(k) errors
are listed as O(h 2 ) errors.
Part of the errors in W at grid point level n can be described by numerical
boundary layers of the form
where A 2 (x; t) is a smooth function that depends only the exact solution u and - depends
only on -. These errors appear at in the computed velocities W. Here
that (14) has a width of a fixed number of grid points in space and so will shrink
as the computation is refined (the size is thus similar boundary layer
will appear at the upper boundary. From now on, we will consider only the bottom
boundary explicitly. It will be shown below (in Lemma 1) that the projection of such a
boundary layer (14) is zero at O(1): This allows us to determine -. The boundary layer
should satisfy the discrete equations (13) exactly to highest order. Inserting (14) into
(13), using Lemma 1 and collecting terms of O(1) (so the differences in the x direction
can be neglected) we obtain
which reduces to a quadratic equation for -:
This equation has two real positive roots for every - ? 0 that occur in reciprocal pairs.
The root with magnitude less than one we denote by - (the other root describes the
boundary layer at the upper wall). The boundary layer (14) does not satisfy the boundary
conditions for \Delta h . In fact, it is generated by a mismatch in the boundary conditions for
W from the global error terms. The details of this are seen below. We note that numerical
boundary layers are normally associated with finite difference methods with wide stencils
that require additional, artificial boundary values to be specified. This is not the case in
(13). The boundary layer that arises in the projection method can be shown formally to
arise from a singular perturbation of the underlying pressure equations with a mismatch
in boundary conditions [9]. We now show the action of the discrete projection on the
boundary layer (14).
has values at grid point that tend asymptotically to
2 (x), a (2) (x; y) and - a (2) (x; y) are smooth functions determined by A 2 .
These functions are also smooth.
We discuss the notation and meaning of this lemma before turning to the proof.
In general, the superscript in brackets denotes the order a term appears and a subscript
denotes the component for a vector quantity. Vectors appear in bold. The term described
by a (2) is a smooth, global regular error term and the term described by - a (2) is an
alternating error term caused by the decoupled stencil for the pressure. Alternating
errors dominate in the pressure errors in Figure 2. What the expansions show is what
will be computed when the discrete projection operator acts on a boundary layer to
high accuracy (a weak but sufficient stability result for the projection step alone can be
derived). Later, we will write
as shorthand for (16). In following lemmas, we present expansions for P h acting on
boundary layers in the horizontal component and regular and alternating terms. We can
then derive an error expansion for W (as noted in the Introduction, a stability argument
for the scheme is still missing so convergence cannot be shown). Here, the terms in the
error expansion will show the order that the various types of errors appear. Expansions
for U and P follow easily. Here and in what follows we retain only the highest order
terms of each type except when necessary to explain some more subtle point. We return
now to the proof of Lemma 1.
Proof of Lemma 1 We denote To satisfy the interior
equations (7) (for Q not P ) the following conditions for the boundary conditions
apply
Where centered differencing of a boundary layer is like
multiplication by -h \Gamma1 ) and the primes denote differentiation in x. These equations
determine the C functions in terms of A
The interior equations for global terms from (7) are
\Deltaq
since there is no global source term. To determine q (2) and -
q (2) we derive Neumann
boundary conditions for them. In [20] it was shown that the effect of the reduced
stencils near the boundary (5), (6) was equivalent to the following two discrete
boundary conditions for Q:
~
where ~
centered differencing
in the x and y directions. These can be interpreted as pressure Neumann
boundary conditions. Using (17) these relationships are both satisfied at O(1). The
action of ~
B on a boundary layer is like multiplication by
smooth terms ~
to third order and for alternating terms, ~
f i;0 to
first order. We note that D y
f actually approximates \Gamma -
f y since centered differencing
uses adjacent grid points of opposite sign. Putting this together we find that
at second order (22) and (23) give the following relationships at
q (2)
y
and
q (2)
y
q (2)
y
These give solvable Neumann data for (19) and (20) and so determine q (2) and -
q (2) .
All the listed terms in the expansion for Q in (16) have been determined. Further
terms in the expansion can be determined similarly.
Having determined Q we can now determine
At O(1) the boundary layers cancel (recall centered differencing in y of a boundary
layer is like multiplication by -h \Gamma1 and note (17)). We then have
q (2)
where r  due to the effect noted above. 2
We have shown Lemma 1 in some detail so the reader can see the idea of the technical
arguments. However, the important features of Lemma 1 are that a vertical boundary
layer is removed (to highest order) by the action of P h and that the boundary layer in
smaller by a factor of h. Later, we will see that there is a boundary layer of
size O(h 2 ) in W. This leads to a boundary layer of size O(h 3 ) in kP
and so a boundary layer of O(h) in P n+1 .
By taking P h of (15) and bringing all the boundary layers to the left hand side we
obtain the following Corollary. The fact that a (2) and - a (2) in Lemma 1 are pure gradients
is used with Lemmas 3 and 5 to show that the global errors are suppressed to fourth
order, although this is not important.
Corollary 2
What we have created is a "pure gradient" boundary layer that has given normal boundary
data at highest order. This is implicitly done in the spatially continuous analysis in
[9].
Lemmas describing the action of P h on regular and alternating terms (denoted by
are stated below. Proofs of Lemmas 3 and 5 can be found in [20].
Lemma 3 When a is a smooth function, P h a has the following error expansion:
when a is compatible. For incompatible a the error terms of both types will appear at first
order.
A compatible function a is one for which the tangential component of Pa also vanishes
on the boundary. A solution u of the Stokes equations and \Deltau are compatible as well
as pure gradient fields. We need a small refinement of this lemma for the error analysis
below. We note that on an alternating term D h \Delta - a approximates r   \Delta -
a. The modified
projection P   describes the projection onto divergence-* free fields with zero normal
boundary values, which is orthogonal to gradient-* fields. Details are given in [20].
Corollary 4 If a is divergence free then a (1) and a (2) are pure gradients. If a is divergence
free and compatible then a (3) is a pure gradient and - a (3) is a pure gradient-*
field.
Proof: We refer the reader to [20] for the details of the proof of Lemma 3 to make
this rigorous, but the idea is simple. We use
a. The error a (2) comes from two sources: rq (2) (a pure gradient) and
the second order errors from computing G h q instead of rq. When a is divergence
so the second source of error is not present and a (2) is a pure gradient.
Similar reasoning applies to the other statements. 2
Lemma 5 The discrete projection acting on an alternating term gives the following error
expansion:
a h- a (1)
The discrete projection acting on an alternating gradient-* field has an expansion beginning
at second order.
We are now in a position to state and prove the main error expansion result for the
intermediate computed velocities W:
Theorem 6 The intermediate velocities have the following error expansion
where u is the exact solution of the Stokes equations. That is, regular errors and vertical
boundary layers begin at second order and alternating errors and horizontal boundary
layers begin at third order.
Proof: For notational simplicity we assume . The divergence-free
and gradient parts of the regular errors are determined at different levels in
the discrete equations (13) so we divide the error terms explicitly
d
We similarly divide the alternating terms into divergence-* free and gradient-*
fields. We insert (27) into the discrete equations (13), expanding \Delta h in a Taylor
series as well as W n about the time level n+ 1. Regular interior terms are collected
below:
a (2) [u]
a (3) [u] (31)
d
\Deltau (2)
+a (2)
d
[u (2)
d
a (2)
d
+a (2)
a (2)
The terms a represent error terms from the discrete projection operator, with square
brackets to denote their source. The interior equations force P h
so
A (3)
using Corollary 4. To determine the equations for alternating terms we use the fact
that the discrete Laplacian amplifies alternating terms in the following way
(see [20]). The third order alternating terms in (13) give
-
boundary conditions are written below for normal component and tangential
component:
+A (2)
d
We will now discuss all of the terms above in detail. Equations (28), (29) and
boundary conditions (36) show that u is indeed the solution of the Stokes equations
we seek so W is a consistent approximation. Equation (30) then determines u (2)
to be
a (2) [u] +rp (38)
is the exact pressure gradient for the Stokes equations. We
know a (2) [u] is a pure gradient from Corollary 4. Once u (2)
is known, A (2)
2 can be
determined from (37) and tangential boundary conditions for u (2)
d
are known and
can be used with the equations (32) to determine u (2)
d
. We note that u
is given in (31). Continuing to ignore the alternating terms for the moment,
the pattern to determine the regular and boundary layer terms is the following:
1. If u (p)
d
is known, u (p+2)
can be determined from the O(h p+2 ) g expansion (i.e.
is determined from (33)).
2. u (p+2)
determines the vertical boundary layer at order p and the tangential
boundary conditions for u (p+2)
d
and (through the effect of the tangential
boundary layers) u (p+3)
d
3. u (p+2)
d
can now be determined from the O(h p+4 ) d expansion (i.e. u (2)
d
is determined
from (32)).
In the discrete setting an important technical detail is that a (2)
d
d
example, in the second line of (32) there is no term a (2)
d
[u (2)
d
]. This is guaranteed
by Corollary 4. This separated determination of the gradient and divergence-free
components of the error expansion for the space continuous analysis is implicitly
present in [9] but not clearly laid out. This technique easily allows for the implicit
handling of the convection terms, for instance, which is avoided in [9].
We turn our attention now to the alternating errors. Equation (35) implies
- u (3)
d
d
- u (3)
In fact, 1-
a (3)
d
u (3) is a pure gradient-*, we can
use Lemma 5 to justify the missing error terms from P h -
u (3) in (32) and (33).
Higher order alternating error terms are determined statically like (39) and (40)
from alternating errors from the projection of lower order terms. An alternating
divergence-* error appears at fifth order. 2
We now turn to the expansions for U and P .
Theorem 7 The computed
W has an error expansion with regular errors at
second order, alternating errors at third order and no boundary layers. The computed
pressure has alternating errors and boundary layers at first order and regular errors at
second order.
Proof: We take P h of (27). The boundary layers are removed (they were so constructed)
and the following results:
d
d
This verifies the first part of the Theorem. Now asymptotically
by (27) minus (41) and results in
a (2) [u])
with a nonzero fourth order regular error a (4) . When (38), (31) and (40) are used
this becomes
The corresponding scalar has an expansion
(recall how boundary layers scale from Lemma 1). Since
we have chosen the convenient scaling for the analysis, we see that
This verifies the second claim of the Theorem. 2
Regularizations
The alternating errors were generated by the uncoupled stencil for D h and G h . Following
[17] we can use higher order regularizing terms (with corrections at the boundary) in D h
and G h to eliminate these alternating errors. A projection operator based on this idea
is described in [20]. We consider the scheme (10)-(12) with this regularized projection:
Theorem 8 The computed velocities W from the spatially regularized scheme
have an error expansion with regular errors at second order, numerical boundary layers
at third order and no alternating errors. The computed pressure has boundary layers at
first order, regular errors at second order and no alternating errors.
Here, numerical boundary layers that occur due to the wide stencil of D h \Delta G h do enter
the computed velocities U. This theorem can be proven using the asymptotic error
description of the regularized projection in [20] following the general framework of the
proof of Theorem 6. We omit the technical details. The presence of the dominant
boundary layer error in the computed pressure for this scheme and the suppression of the
alternating errors can be seen in Fig. 4 (compare to Fig. 2 for Chorin's original scheme
with the same h and k).
Using the regularized D h and G h as discussed above, we can further eliminate the
dominant boundary layer errors in the pressure by using a non-split-step scheme:
As shown in the theorem below, this scheme suppresses the numerical boundary layers
from the projection step.

Table

2: Normalized pointwise pressure errors e p (and estimated convergence rates in h)
from the fully regularized scheme.
Theorem 9 The computed velocities from the scheme (42),(43) with spatially regularized
D h and G h have an error expansion with regular errors at second order, numerical
boundary layers at fourth order and no alternating errors. The computed pressure has
regular errors at second order, numerical boundary layers at third order and no alternating
errors.
Again, we omit the technical details. We note that the scheme (42),(43) requires the
solution of a coupled system for U n+1 and P n+1 . As mentioned in the introduction, it
would be more computationally efficient to use a pressure increment scheme [18, 6] to
suppress the numerical boundary layers, but the analysis of this approach is not fully
understood in the discrete setting.
Second order convergence for the pressure from the fully regularized scheme is shown
in

Table

2 using the same parameters as the convergence study from Section 3 (compare

Table

1). The smooth nature of the errors in computed pressure is shown in Fig. 5.
6

Summary

We have presented an error analysis for Chorin's original fully discrete method for computing
the incompressible Navier-Stokes equations. The velocities from this scheme converge
with full order O(k)+O(h 2 ). The computed pressures have O(h) global alternating
errors due to the uncoupled approximation used for the incompressibility condition and
layers due to the split-step projection step. These errors can
be removed by using a regularized stencil to approximate the incompressibility condition
and a non-split-step time integration procedure.



--R

"Derivation and solution of the discrete pressure equations for the incompressible Navier Stokes equations,"
Spectral Methods in Fluid Dynamics (Section 11.3)
"Numerical solution of the Navier-Stokes equations,"
"On the convergence of discrete approximations to the Navier-Stokes equations,"
"A Mathematical Introduction to Fluid Dynamics,"
"A second-order projection method for the incompressible Navier-Stokes equations,"
"An efficient second-order projection method for viscous incompressible flow,"
Finite Element Methods for Navier-Stokes Equations
"Projection Method I: Convergence and Numerical boundary Layers,"
"Projection Method II: Rigorous Godunov-Ryabenki Analysis,"
"Second Order Convergence of a Projection Scheme for the Incompressible Navier-Stokes Equations with Boundaries,"
"Discrete Compatibility in Finite Difference Methods for Viscous Incompressible Flow,"
"Application of a fractional-step method to incompressible Navier-Stokes equations,"
"boundary conditions for incompressible flows,"
"On Chorin's projection method for the incompressible Navier-Stokes equations,"
"On error estimates of projection methods for Navier-Stokes equations: first order schemes,"
"High-Order Accurate Schemes for Incompressible Viscous Flow,"
"A second-order accurate pressure-correction scheme for viscous incompressible flow,"
"Finite Difference Vorticity Methods"
"Analysis of the spatial error for a class of finite difference methods for viscous incompressible flow,"
--TR

--CTR
Robert D. Guy , Aaron L. Fogelson, Stability of approximate projection methods on cell-centered grids, Journal of Computational Physics, v.203 n.2, p.517-538, 1 March 2005
Weinan , Jian-Guo Liu, Projection method III: spatial discretization on the staggered grid, Mathematics of Computation, v.71 n.237, p.27-47, January 2002
parasitic modes;numerical boundary layers;projection methods