I'm using a semi-implicit staggered grid(MAC H. Harlow and J. Eddie.) finite difference approximation of incompressible Navier-stokes equations . The convection and diffusion terms are discretized using AB2 and Crank-Nicolson methods respectively, and for the pressure velocity decoupling a first order projcetion method is used where a Poisson equation is solved using neumann boundary conditions and a constant (0) pressure in one node.

All linear systems of equations are solved using a Bicgstab method with ilu preconditioning. However, in a simple lid-driven cavity problem I face a numerical instability which in fact makes all variables diverge after 6 or 7 timesteps (bicgstab solver works just fine until the values become very large) ,while the CFL number is below 1.

Space discretization scheme is similar to that of H. Harlow and J. Eddie.

I rechecked everything many times and yet I'm not able to find the problem.