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:
- 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.
- 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 $$
- 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.
