Pressure projection methods

[fb77 (François Bouchut)]

The equation you have to solve for p^{n+1} is a Laplace equation with non homogeneous Neumann boundary condition.
You can look at this post

Uzawa scheme for stokes bingham flow

In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct). You can modify easily your code to do the incompressible case, by using the extra variables p,q as solve bing([ux,uy,p],[vx,vy,q]) and adding -int2d(Th)(p*Divergence(vx, vy)+q*Divergence(ux, uy)) But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmente…

From “We are thus led to the problem…”