Jump to content

  • Log In with Google      Sign In   
  • Create Account


[solved] SPH fluids "boiling", not colliding properly?


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
4 replies to this topic

#1 trs79   Members   -  Reputation: 126

Like
0Likes
Like

Posted 14 December 2011 - 09:14 PM

Hey all,

I'm seeing some strange behavior from some SPH fluid code I'm working on. I modified it to have a cylinder as a container/collision shape and that works fine. However, I then change it so that it slopes inward as it gets taller, and the particles start to move around constantly like the fluid is boiling and they never settle down. Also, the fluid level doesn't get higher as the container slopes inward as I'd expect. See image:

fluid.gif

Here is the code I'm using to calculate positions (it's a modified form of the code at http://www.rchoetzle...hics/fluids.htm)


accel = p->sph_force;
accel *= m_Param[SPH_PMASS];
 	.....
// p is the particle structure, with z being up
// radius is the particle radius
// norm is a 3-vector

   	float boundary = abs((-.4 * p->pos.z) + 9.8);

		diff = 2.0f * radius - (boundary - (float)fastSqrt_2(p->pos.x*p->pos.x + p->pos.y*p->pos.y)) * ss;

		if (((diff < 8 *radius) && (diff > EPSILON)))
		{
			norm.Set ( -p->pos.x, -p->pos.y, 0 );
			norm.Normalize();
			
			adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
			accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
		}

   ....

// Leapfrog Integration ----------------------------
		vnext = accel;							
		vnext *= m_DT;
		vnext += p->vel;						// v(t+1/2) = v(t-1/2) + a(t) dt
		p->vel_eval = p->vel;
		p->vel_eval += vnext;
		p->vel_eval *= 0.5;					// v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5		used to compute forces later
		p->vel = vnext;
		vnext *= m_DT/ss;
		p->pos += vnext;						// p(t+1) = p(t) + v(t+1/2) dt
 

I'm pretty sure the way the normal is being added to the acceleration is causing the "boiling" of the particles, but I'm not sure. Thanks for any help

Sponsor:

#2 RulerOfNothing   Members   -  Reputation: 1149

Like
0Likes
Like

Posted 15 December 2011 - 12:19 AM

Here I am assuming the norm variable holds the normal vector of the container, which would indicate to me that it needs to have a non-zero z value if the container is a frustrum instead of a cylinder.

#3 trs79   Members   -  Reputation: 126

Like
0Likes
Like

Posted 15 December 2011 - 09:44 AM

Here I am assuming the norm variable holds the normal vector of the container, which would indicate to me that it needs to have a non-zero z value if the container is a frustrum instead of a cylinder.


Thanks for the reply, you are correct the norm variable holds the normal vector of the container. Unfortunately when I change it to have a non-zero z value I still get the erratic particle movement, i.e :

norm.Set ( -p->pos.x, -p->pos.y, p->pos.z );
norm.Normalize();

Maybe I need something different for the z value?

#4 luca-deltodesco   Members   -  Reputation: 637

Like
0Likes
Like

Posted 15 December 2011 - 10:34 AM

do you understand what the normal to a surface means?

you've been calculating the normal from the particle position, which is a little strange but I guess it works if you know the centre of the cylinder is at (0,0).

The equation for a cone with it's apex at the origin and axis parallel to z with aperture 2θ is:
f(u) = u.z - |u|.cos θ = 0

Since your cone won't have it's apex at the origin, you can simply offset u.z by -h for apex at height 'h'.

Letting d = sqrt(x2+y2+(z-h)2)
the partial derivatives of f(u+h.z) are:
pf/px = -(x/d).cos θ, pf/py = -(y/d).cos θ, pf/pz = 1 - ((z-h)/d).cos θ

which are the components of your normal (to be normalised)



you can speed up the normalisation as we can determine the magnitude of the vector (pf/px, pf/py, pf/pz) to be sin θ as long as you know the point is on the surface of the cone

#5 trs79   Members   -  Reputation: 126

Like
0Likes
Like

Posted 16 December 2011 - 10:29 AM

do you understand what the normal to a surface means?

I do, but this is my first attempt at coding collision detection, probably why I got confused on the implementation :(

Thanks for the help, it was the z component that was causing problems, I finally ended up using this based on your code:
norm.Set ( -p->pos.x/diff, -p->pos.y/diff, 1.0-(p->pos.z/diff) );

along with the initial boundary slope code:
float boundary = abs((-.4 * p->pos.z) + 9.8);

Thanks again!




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