I have been trying to make a 2D fluid simulation with the SPH technique after reading http://image.diku.dk.../kelager.06.pdf
I could follow everything in the paper/thesis, so I thought that's a piece of cake.
But after several days of trying it still doens't work right. The problem I have is in the pressure term.

All the particles are pulled down because of the gravity, and this works correct.
So then I added collisions with the wall and those were working correct too.
Then I tried to add pressure and it doesn't work. The pressure term is supposed to push the particles away from each other.
Like in this video : [media][/media]
In here the particles create several layers, but with my code I get only one layer.
It's as if the particles don't push each other away. So I tried to increase the density and the gas stiffness term.
But both didn't work.

This is what I got:
http://imageshack.us/photo/my-images/513/sph.png/

I'm a bit clueless as what to do now, because this is my first time simulating a fluid.
I really wanted to make a SPH simulation to eventually do star formation simulations.
Math is no problem for me as I am a 3rd year bachelor in physics and astrophysics.

I've programmed it in c++ ( I started in python but it was really slow ).

void Simulation::timestep( ){

//Calculating density field
for( unsigned int i = 0; i < particles.size(); i++ ){
for( unsigned int j = 0; j < particles.size(); j++ ){
float W6 = W_pol6( particles[ i ].position - particles[ j ].position, h );
particles.density += particles[ j ].mass * W6;
}
//Finding the maximum density to render a density map
if( particles.density > density_max )
density_max = particles.density;
particles.pressure = k*( particles.density - density );
particles.clear_force();
}

std::vector<Plane> planes = container.get_planes();

for( unsigned int i = 0; i < particles.size(); i++ ){

//Adding internal forces
float density_i = particles.density;

for( unsigned int j = (i+1); j < particles.size(); j++ ){
if( j != i ){

//Adding pressure
// P = pi/ri^2 + pj/rj^2
float P = particles.pressure / ( density_i * density_i );
P += particles[j].pressure / ( particles[j].density * particles[j].density );
Vector2 gradW = grad_W_spike( particles.position - particles[j].position, h );
Vector2 force_pressure = -gradW*P;
particles.add_force( force_pressure*particles[j].mass*particles.density );
particles[j].add_force( -force_pressure*particles.mass*particles[j].density );
}
}

//Integrating with leapfrog
particles[ i ].last_velocity = particles[ i ].velocity;
particles[ i ].velocity += ( particles.force / particles.density + gravity ) * dt;
particles[ i ].position += ( particles.velocity ) * dt;

//Collision With Walls
//Collision response by projection
for( unsigned int j = 0; j < planes.size(); j++ ){
if( particles.velocity.dot( planes[j].normal ) < 0 ){
float distance = particles.position.dot( planes[j].normal ) + planes[j].distance;
if( distance - radius < -threshold ){
particles.velocity -= planes[j].normal * planes[j].normal.dot( particles.velocity ) * ( 1.0 + epsilon );
particles.position -= ( planes[j].normal * ( distance - radius ) )*(2.0 );
}
}
}
}

}


Is there any particular reason you're dividing the pressure force by the particle density?

particles[ i ].velocity += ( particles[ i ].force / particles.density + gravity ) * dt;

That means that as particles pile up, the pressure force will decrease. That's the opposite of what you want to happen.

You should instead divide the force by each particle's mass, such that F = ma. Density is not a substitute for mass in this case. Density shouldn't even be considered as a per-particle value - you're only storing it so that you can sample the density field via a smoothing kernel.

Simply removing that reference to density in the integration step should get you some fluid-like behaviour.