When using the LBGK equation
f[sub]a[/sub](x+c[sub]a[/sub], t+1) = (1-omega)f[sub]a[/sub](x, t) + omegaf[sub]a,eq[/sub]
for stability it is derived in e.g. the book by Succi that
0 < omega < 2
If, however, 1 < omega < 2 it could (mathematically) happen for certain f that in the next timestep it becomes negative, which is of course unphysical. Is it realistic to expect that this could indeed happen in a simulation or will f always be close enough to f[sub]eq[/sub] to ensure that negative populations can never occur?
the populations can indeed become negative, which is unphysical. This is also a reason why simulations may become unstable for very small lattice viscosities. You should choose the simulation parameters in such a way that the populations are always positive.
As soon as you ensure that your populations are positive you will always obtain the stable simulations with LBGK. However, in reality we don’t know when it becomes unstable. There is a complex dependancy of achieved velocity on relaxation parameter. Usually if you have unstable situation, you need to reduce velocity by adjusting grid, body force, etc.
However, it’s a too oversimplified recipe and doesn’t work too well. In this case you need to consider whole stability picture.