Jump to content

  • Log In with Google      Sign In   
  • Create Account

SPH simulation - pressure force doens't work


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
1 reply to this topic

#1 Sepii   Members   -  Reputation: 200

Like
0Likes
Like

Posted 09 September 2012 - 08:37 AM

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 : http://www.youtube.com/watch?v=82efx2eDOhM
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[i].density += particles[ j ].mass * W6;
		}
		//Finding the maximum density to render a density map
		if( particles[i].density > density_max )
			density_max = particles[i].density;
		particles[i].pressure = k*( particles[i].density - density );
		particles[i].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[i].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[i].pressure / ( density_i * density_i );
				P += particles[j].pressure / ( particles[j].density * particles[j].density );
				Vector2 gradW = grad_W_spike( particles[i].position - particles[j].position, h );
				Vector2 force_pressure = -gradW*P;
				particles[i].add_force( force_pressure*particles[j].mass*particles[i].density  );
				particles[j].add_force( -force_pressure*particles[i].mass*particles[j].density );
			}
		}

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

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

}


Sponsor:

#2 Helicobster   Members   -  Reputation: 195

Like
1Likes
Like

Posted 15 September 2012 - 12:28 AM

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.

Edited by Helicobster, 15 September 2012 - 12:31 AM.





Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS