I'm trying to solve variable density and viscosity Navier-Stokes equation using lagged pressure projection method. I'm solving for cavity problem as a test case now (once I get projection right, I will be solving a multiphase problem). The problem is I get oscillations developed in the solution as I march in time. .

I'm using collocated grid for solving this problem.

I'm solving variable density, viscosity N-S equation: $$ \dfrac{u^* - u^n}{\Delta t} + \left(u_i\frac{\partial u_i}{\partial x_j} \right)^n = \left(-\frac{1}{\rho}\frac{\partial p}{\partial x_i} \right)^n + \frac{1}{\rho} \left[\frac{\partial}{\partial x_j}\left( \mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right) \right]^n $$

and on boundary, $U^* = U^{n+1}$.

I'm approximating convective terms with 1st order godunov as of now, diffusion term and pressure gradient with central difference.

Then I'm solving the projection step,

$$ \frac{\nabla.u^*}{\Delta t} = \nabla.(\frac{1}{\rho}\nabla\phi) $$

using boundary condition $$ \frac{\partial \phi}{\partial n} = 0 $$

I'm discretizing $\nabla.u^*$ as $\dfrac{u^*_{i+1,j} - u^*_{i-1,j}}{2\Delta x} + \dfrac{v^*_{i,j+1} - v^*_{i,j-1}}{2\Delta y}$ and applying boundary condition using ghost node. and for pressure, I'm discretizing as $\dfrac{\left(\frac{1}{\rho}\frac{\partial p}{\partial x} \right)_{i+1/2,j} - \left(\frac{1}{\rho}\frac{\partial p}{\partial x} \right)_{i-1/2,j}}{h} + \dfrac{\left(\frac{1}{\rho}\frac{\partial p}{\partial y} \right)_{i,j+1/2} - \left(\frac{1}{\rho}\frac{\partial p}{\partial y} \right)_{i,j-1/2}}{h}$

For solving the singular poisson equation, I'm removing the nullspace from the RHS ($\nabla.u^*$), operator matrix ($\nabla.(1/\rho$)). I'm also removing the nullspace from the solution vector.

I'm using petsc to solve the poisson equation and using GMRES ksp solver.

I,m updating velocity using $U^{n+1} = U^* - \Delta t*(\nabla\phi)$ and $P^{n+1} = P^* + \phi$. I'm fine with first order accurate pressure.

However, I get the oscillations in my solution. And then the solver diverges if kept iterating. I strongly feel there is something wrong with the boundary conditions, especially on $U^*$ as the solution is very sensitive to what I specify $U^*$ on boundaries. I also tried with $U^* = 0$ but no improvement.

Is there anything wrong in the method I'm using here? It'll great help if you can point it out.

Thank you.

See results here: Results

Found the solution:

I'm aware of the fact that the staggered grid is the most straight forward solution to this checkerboard problem but I wanted to stick to the collocated grid as it is easier to code and less prone to error, plus having velocities at cell centers will greatly reduce the level set formulation efforts compared to staggered grid.

Now, instead of lagged-pressure projection method, I've changed it to the (Kim and Moin's) pressure-free projection method which is working great as of now. However, I've coded it for constant density and viscosity Navier-Stokes I'm hoping it to work for the variable density and viscosity as well.

While implementing, now there is no pressure gradient term in the momentum equation and I used following boundary conditions on intermediate velocity field: $$ U^*.\hat{n} = U^{n+1}.\hat{n} $$ $$ U^*.\hat{t} = (U^{n+1} + \Delta t*\nabla\phi).\hat{t} $$ where $\hat{n}$ and $\hat{t}$ are normal and tangent vectors. And poisson equation has the neumann B.C. as before.

As I have to get going with my project, I have not thoroughly investigated the error in lagged-pressure projection method But I will do so as soon as I get some time.

Thank you.

Here's the solution (at t = 10sec with Re400) with pressure-free projection method: quiverPressure

  • 2
    $\begingroup$ You need to find out patterns to debug this. For example, what happens when you make the mesh size smaller? Do the oscillations become stronger or weaker? What happens if you make the time step smaller? Observing patterns in the error is often a useful step in finding out where the problem lies. $\endgroup$ Aug 11, 2014 at 22:09
  • $\begingroup$ i am actually working on this question, I appreciate you posted it, thank you I understand this in theory, however as I try to code it up I have two problems, 1. how do you remove the null space from the singular matrix, 2. as for the boundary condition for intermediate U^*, the first, we do not know pressure yet, how to specify the boundary value as U^(n+1)-dt*grad(p), we have't calculated it yet. thank you really hope you can reply me. $\endgroup$
    – user19751
    Mar 25, 2016 at 1:41
  • $\begingroup$ @sunny, 1) I used PETSc library to remove null space (and solve). 2) U can used pressure from previous time step (called as lagged pressure projection method, refer projection method by Bell) or U can neglect the pressure from U* altogether (call pressure-free projection method, refer fractional step method by Kim and Moin). Note, the correction terms are different in both the methods. Choose whichever suits your needs. $\endgroup$
    – Pranav
    Apr 4, 2016 at 20:21
  • $\begingroup$ Are you sure it is odd-even decoupling? I assume that accepted answer is not what you wanted to accept, not a big deal, slightly confusing. I also read in your question you answered yourself, and the solution that you find is using non-incremental scheme. I can agree with that part, but are you sure it is stable unconditionally? I believe not. You may find it interesting section 2.6 and table II in paper from [Guermond and Quartapelle](math.tamu.edu/~guermond/PUBLICATIONS/guerm $\endgroup$
    – trblnc
    May 16, 2016 at 13:16
  • $\begingroup$ it's not unconditionally stable because I've used explicit method. $\endgroup$
    – Pranav
    May 17, 2016 at 5:16

2 Answers 2


What you are observing is odd-even decoupling between the velocity and pressure. This is a well-known problem on collocated grids. At any point on the grid, the gradient of the pressure is calculated using $p_{i+1}$ and $p_{i-1}$:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + (u_i^n.\nabla)u_i^n = -\frac{p_{i+1}-p_{i-1}}{2\Delta x}+\nu\nabla^2 u_i^n$$

Hence the velocity $u_i$ only depends on the pressure at the adjacent grid points, and not on $p_i$ itself. This can cause the numerical pressure to behave strangely. As you can see in your results, a "checkerboard" pattern is observed, where only the pressures at alternate grid points influence each other.

This problem can be solved by using a staggered grid, where the velocities and pressures are calculated at different points in the domain. The most commonly used staggered grid on a Cartesian mesh is the one introduced by Harlow and Welch [1], but it will make your code more complex. If you want to make minimal changes to your code, then use the same collocated grid for the velocities, but compute the pressures at the center of each cell. Take care to set up your pressure equation carefully, and interpolate the velocities as required when you calculate the right-hand side $\nabla\cdot u^*$.

Here is another page that discusses this phenomenon: http://www.imm.dtu.dk/~pcha/Projekter/Shock/main/node17.html

Regarding the boundary conditions for the intermediate velocity field, rather than $u^*=0$, you should be using $u^*=u^{n+1}$. In the problem you have considered, this makes a difference on the top surface where the lid is moving. For more information about the boundary conditions using a projection method, read the paper by Perot [2] and follow the references.

[1] Harlow, Francis H., and J. Eddie Welch. "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface." Physics of fluids 8.12 (1965): 2182.

[2] Perot, J. Blair. "An analysis of the fractional step method." Journal of Computational Physics 108.1 (1993): 51-58.


I think your problem lies in the boundary conditions for $\mathbf{u^*}$. Since your grid is collocated, $\mathbf{u^*}$ is stored at cell centers; and the boundary lies half a cell distance from the place where $\mathbf{u^*}$ is actually computed. Thus a higher-order BC for $\mathbf{u^*}$ should do the trick. Also, I would like to mention that the Kim-Moin projection method does not yield a high order of accuracy in pressure [1]. [1] will also help with many other order of accuracy issues you might face.

[1] Brown D. L., Cortex R. and Minion M. L. "Accurate projection methods for the incompressible Navier-Stokes equations"

  • $\begingroup$ For Kim and Moin's pressure-free projection method, we dont really need pressure to advance in time. However the pressure can be recovered if needed. Also, for my present project (which I have now completed successfully), I do not need pressure as no sources depends on it. I only need accurate velocity field to advect the interface. Also I did not follow your suggestion of using higher order B.C. for u*. I don't need to compute u* at boundaries as there will be B.C. specified. (Btw, I'm using finite difference and not Finite Volume.) $\endgroup$
    – Pranav
    Aug 29, 2014 at 4:30

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.