Is the projection method of integrating Navier-Stokes equations exact?

Take the incompressible flow equations:

$$ \frac{\partial\mathbf{u}}{\partial t} = -\mathbf{u}\cdot \nabla \mathbf{u} -\nabla p + \mu \Delta \mathbf{u} \\ \nabla\cdot\mathbf{u} = 0 $$ along with some meaningful boundary conditions. Using a projection operator $\mathbf P$ defined as $$ \mathbf{P} = I -\nabla (\Delta)^{-1}\nabla $$ the earlier equations becomes an ODE like so: $$\frac{\partial\mathbf{u}}{\partial t} = \mathbf{P}(-\mathbf{u}\cdot \nabla \mathbf{u} + \mu \Delta \mathbf{u}) $$ because $\mathbf{u}$ is solenoidal.

Is there any approximation is involved here or are the two sets of equations exactly same? And, can one discretize and integrate the ODE-looking equation with arbitrary accuracy in space and time using their choice of discretization and integrator because it is an ODE? Is finding (analytical) $\mathbf{P}$ the tricky part here?

I know what might be done in practice in writing a finite difference code e.g. Adams-Bashforth + Crank-Nicholson. A discrete projection operator is obtained by solving a Poisson equation

  • 1
    $\begingroup$ One of the complications is that you have to define what $(-\Delta)^{-1}$ represents. It's obviously the inverse of the Laplace operator, but corresponding to what boundary conditions? $\endgroup$ Sep 4, 2015 at 22:04

2 Answers 2


The projection scheme is indeed analytic, provided that the initial conditions are incompressible. You can derive it by taking the divergence of the evolution equation for $u$.

Finding the projection P can be tricky; it's easier with some methods like spectral / pseudospectral methods. But one has to be careful that numerical error doesn't move the solution off of the incompressible manifold; these errors can add up. Boundary conditions can also be tricky, because P has to accommodate for those as well. Also, the resulting equation is a PDE; the RHS has derivatives in space, so to advance it in time one must perform some sort of spatial discretization to transform it into a system of ODEs.

  • $\begingroup$ Is the special warning about the numerical error because in those methods the divergence free constraint is not explicitly enforced at each step?. I realize that even after the projection it is still a PDE. Dont know what was going through my head. $\endgroup$ Sep 3, 2015 at 18:14
  • $\begingroup$ Yes, that's correct; if the divergence-free condition is not explicitly enforced, some divergence might enter due to round-off error. If this may experience secular growth, you're going to have a bad time. $\endgroup$
    – Malcolm
    Sep 17, 2015 at 12:50

Note: I originally posted an answer that I was not 100% pleased with, so I have revised it heavily. 1/20/2017

The projection method is not in an exact approximation to the full system in general. There are a few reasons for this, but the most immediate has to do with the boundary conditions. Let's look at the continuous setting first. Consider the unsteady Navier-Stokes problem in a sufficiently regular domain $\Omega$ with appropriate boundary conditions, given by:

\begin{align} \begin{split} \label{UnsNS} \frac{\partial \boldsymbol{u}}{\partial t} -\nu\nabla \cdot (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^T) + ( \boldsymbol{u} \cdot \nabla) \boldsymbol{u} + \triangledown p = \boldsymbol{f} \quad \textit{in } \Omega \times \lbrack t_0, t_n \rbrack \\ \nabla \cdot \boldsymbol{u} = 0 \quad \textit{in } \Omega \times \lbrack t_0, t_n \rbrack \\ \mathcal{B}\left( \boldsymbol{u}, p \right) = \boldsymbol{g} \quad \textit{on } \partial \Omega \times \lbrack t_0, t_n \rbrack \\ \boldsymbol{u}(\boldsymbol{x},0) = \boldsymbol{u}_0 \quad \textit{at} \quad t=t_0 \end{split} \end{align} Where $\mathcal{B}$ is an appropriate trace operator describing the boundary conditions. For simplicity, let's assume that $\boldsymbol{u}=0$ on $\partial \Omega$. At a step $n$ the standard Chorin-Temam projection method consists of three steps:

  1. We introduce an intermediate velocity field $\tilde{\boldsymbol{u}}^n$ and separate the velocity and pressure computations in two parts. We compute $\tilde{\boldsymbol{u}}^n$ as follows:

\begin{align*} \frac{1}{\Delta t} (\tilde{\boldsymbol{u}}^n - \boldsymbol{u}^{n-1} ) - \nu \Delta \tilde{\boldsymbol{u}}^n + (\tilde{\boldsymbol{u}}^n \cdot \nabla ) \tilde{\boldsymbol{u}}^n = \boldsymbol{f}^n \\ \tilde{\boldsymbol{u}}^n =0 \, \, \text{on} \, \, \partial \Omega \end{align*} Note we enforce boundary conditions here.

  1. Next, we use the intermediate velocity field to compute the pressure: $$\frac{1}{\Delta t} (\boldsymbol{u}^{n} - \tilde{\boldsymbol{u}}^{n} ) = -\nabla p^n$$ By exploiting the fact that $\nabla \cdot \boldsymbol{u}^n = 0$, we can apply the divergence operator to both sides and formally eliminate $\boldsymbol{u}^n$. We now have a Poisson problem for the Pressure: $$ -\Delta t \Delta p^{n} = - \nabla \cdot \tilde{\boldsymbol{u}}^n $$
  2. Apply the Helmholtz Principle to get the final velocity $\boldsymbol{u}^n$: $$\boldsymbol{u}^n = \tilde{\boldsymbol{u}}^n - \Delta t \nabla p^n$$

There are a couple of things to notice here. The first is that we only have a second order operator for the velocity in step (1), and accordingly we can only control the trace operator (and hence the boundary conditions) during this step. The next thing to notice is that during step (2), we have a second-order operator for the pressure. This implicitly introduces a boundary condition on the pressure; in this case we have imposed homogenous Neumann boundary conditions for pressure. Such boundary conditions may not be accurate for the true pressure.

Because our final velocity computed in step (3) relies on step (2), our boundary conditions are no longer guaranteed to be satisfied (as we enforced them during the first step). Some people have suggested that because of this the intermediate velocity field $\tilde{\boldsymbol{u}}^n$ should be regarded as the 'true' velocity field; however this is equally problematic because the intermediate field is not divergence free in general. So we have two choices: a velocity field that respects incompressibility but not the boundary conditions, and a velocity field that respects the boundary conditions but not incompressibility.

Of course, in the numerical setting there are other reasons why this solution is not "exact"; in particular the splitting error induced by separating the velocity and pressure by means of the intermediate velocity field is nonzero. You can observe this easily by looking at the problem in the algebraic setting. Our standard N-S matrix is given by: \begin{align} \begin{split} \label{UnsAlgProb} \begin{bmatrix} A & B^{T} \\ B & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{u}^{n} \\ \boldsymbol{p}^{n} \end{bmatrix} = \begin{bmatrix} \boldsymbol{f} + \frac{1}{\Delta t} \boldsymbol{u}^{n-1} \\ \boldsymbol{0} \end{bmatrix} \end{split} \end{align} However, for Chorin-Temam we are solving the following block LU system:

\begin{align} \begin{split} \label{projA} \begin{bmatrix} A & 0 \\ B & -BB^T \end{bmatrix} \begin{bmatrix} I & \Delta t B^T \\ 0 & I \end{bmatrix} \begin{bmatrix} \boldsymbol{u}^{n} \\ \boldsymbol{p}^{n} \end{bmatrix} = \begin{bmatrix} \boldsymbol{f} + \frac{1}{\Delta t} \boldsymbol{u}^{n-1} \\ \boldsymbol{0} \end{bmatrix} \end{split} \end{align} Or, similarly, reinterpreting the method as the algebraic projection method (noting that the mass matrices $M$ correspond to the identity operator) we have a similar block LU system:

\begin{align} \begin{split} \label{algprojA} \begin{bmatrix} A & 0 \\ B & -BM^{-1}B^T \end{bmatrix} \begin{bmatrix} I & \Delta t M^{-1}B^T \\ 0 & I \end{bmatrix} \begin{bmatrix} \boldsymbol{u}^{n} \\ \boldsymbol{p}^{n} \end{bmatrix} = \begin{bmatrix} \boldsymbol{f} + \frac{1}{\Delta t} M\boldsymbol{u}^{n-1} \\ \boldsymbol{0} \end{bmatrix} \end{split} \end{align}

In both cases, the difference between the 'exact' system and the projected one is nonzero. This can be verified easily by performing the matrix multiplication.

I will also note that if one replaces the matrix $\Delta t M^{-1} B^{T}$ with $\Delta t A^{-1} B^{T}$ in the (1,2) block of the upper triangular matrix for the algebraic projection method, one can regain control of the boundary conditions (as $A$ corresponds to a second-order operator). This resolves our issue with boundary conditions, though the splitting error remains nonzero.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.