(Approximate) Incremental Projection Method for Navier-Stokes equations

I am trying to implement an incremental projection method for the 2D incompressible Navier-Stokes. The type of projection method I am trying is

$$\frac{u^{*} - u^{n}}{dt} = - \nabla p^{n} - u \cdot \nabla u + \nu \nabla^{2} u\\ u^{*}|_{\partial D} = u^{n+1}|_{\partial D}\\ \nabla^{2} \phi^{n+1} = \frac{\rho}{dt} \nabla \cdot u^{*}\\ \nabla \phi^{n+1}|_{\partial D} = \frac{1}{dt} \left(u^{*} - u^{n+1}\right)\\ \frac{u^{n+1} - u^{*}}{dt} = - \nabla \phi^{n+1}\\ p^{n+1} = p^{n} + \phi^{n+1} (5)$$

In terms of the velocity-pressure arrangement I am using a collocated formulation using the Rhie-Chow interpolation to have stable coupling. I however have run into a problem in terms of the practical implementation of the pressure update. The way I update pressure is as follows

Pnew(1:nx,1:ny) = P(1:nx,1:ny) + Phi(1:nx,1:ny)
P(1:nx,1:ny) = Pnew(1:nx,1:ny)


I am testing my solver on the lid driven cavity. My velocity profile is correct, however my pressure profile is completely wrong as it contains wiggles. Here is a picture:

I have tested my solver on a pressure-free projection method(where there is no $\nabla p^{n}$ in the intermediate veloctiy step and it works fine.

EDIT: (1/17/18)

I forgot to mention but I am using 2nd-order central discretization for both convection and diffusion. My integration technique is currently forward Euler for debugging purposes.

• There are so many reasons why this could be wrong. There are also many reasons why it could be right -- have you investigated whether the pressure solution you get actually converges, and then whether it converges to the correct pressure, and then whether it converges at the right rate? Commented Jan 17, 2018 at 20:54
• The pressure profile in my incremental projection method I know to be wrong because it does not agree with the results of my pressure-free projection method in which I have validated with data. With that said, I do find it strange that while the pressure profile has wiggles the velocity profile is correct. In terms of the rate of convergence, in time I am getting O(h^2,dt) which corresponds to the discretization methods I am using. I have edited my post to contain more details on my numerical schemes. Commented Jan 17, 2018 at 22:05
• Bump. Any thoughts? Commented Jan 23, 2018 at 12:14