I have been going over the theoretical foundations quite a bit this week. Reading papers and generally taking note of what different people do for boundary conditions. I have also done quite a bit of testing with very small particles and large particles.
I found one very big problem. At boundary and free surfaces, the density is "smoothed". This creates big pressure gradients at the surface and can cause big problems for initial conditions. However this is well known in the SPH literature and is called the particle consistency problem.
One method of fixing it is to consider the Taylor series and you get a reasonably simple expression for density. But the equations become implicit, and thus solving systems of equations are required. Instead I use the correction from the previous time step and thus is just a estimate. This works very well for density, but when I use the corresponding expression for the grad operator, I get nothing but instability. The expression can be 0/0 and the correct limit is 0, but numericaly.....
So at this point I only use the correction for the density and I'm getting some very nice results. But I will do some more work on proper correction of partial derivatives and attempt to get that working.