• Create Account

We need 7 developers from Canada and 18 more from Australia to help us complete a research survey.

Support our site by taking a quick sponsored survey and win a chance at a \$50 Amazon gift card. Click here to get started!

[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.

4 replies to this topic

#1trs79  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:

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
// 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 );
}

....

// 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

#2RulerOfNothing  Members   -  Reputation: 1320

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.

#3trs79  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?

#4luca-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

#5trs79  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