[ODE] Approximate Hydrodynamic Forces

Ander Taylor ander_taylor at hotmail.com
Tue Apr 15 23:57:02 2003


Cool, I will try that when I get home.

Thankyou man : )

I like the sound of the Drag Joint.

Ander


David wrote:

I gave the docs another look and it looks like I did get the l, w, and h
out of order... so give this a try:

const dReal *lvel = dBodyGetLinearVelocity(body);
const dReal *R = dBodyGetRotation(body);
dReal ss[3];
dGeomBoxGetLengths(geom, ss);

dReal nx = ss[1]*ss[2] * fabs(R[0]*lvel[0] + R[1]*lvel[1] + R[2]*lvel[2]);
dReal ny = ss[0]*ss[2] * fabs(R[4]*lvel[0] + R[5]*lvel[1] + R[6]*lvel[2]);
dReal nz = ss[0]*ss[1] * fabs(R[8]*lvel[0] + R[9]*lvel[1] +R[10]*lvel[2]);

//either this....
dReal temp = -nx*scale;
dBodyAddForce(body, temp*R[0], temp*R[1], temp*R[2]);
temp = -ny*scale;
dBodyAddForce(body, temp*R[4], temp*R[5], temp*R[6]);
temp = -nz*scale;
dBodyAddForce(body, temp*R[8], temp*R[9], temp*R[10]);

//this might actually work now with the fabs()'s (may need to adjust scale):
dReal temp = -scale*(nx+ny+nz);
dBodyAddForce(body, temp*lvel[0], temp*lvel[1], temp*lvel[2]);


Let me know how that works for you.

David

P.S. You may need a tiny bit more code for a really good simulation... you
need to simulate bouyancy (force of buoyancy = weight of displaced fluid =
density of fluid * volume of object * gravity):

dMass m;
dBodyGetMass(body, &m);
dReal gravity[3];
dWorldGetGravity(world, gravity);
//ss, temp already defined above

temp = -fluidDensity * ss[0] * ss[1] * ss[2];
dBodyAddForce(body, temp*gravity[0], temp*gravity[1], temp*gravity[2]);

//which, unless you have a non-axis-aligned gravity, is probably just:
dBodyAddForce(body, 0, temp*gravity[1], 0);

Note that this method assumes the body is completely submerged in the
fluid.  A more generic solution is not all that complicated, I just
haven't thought it out yet.

Hmmm, now I'm starting to get more ideas about a 6 dof rolling friction
contact joint and a 3 dof drag joint that could be integrated into the
collision system... If only there were more hours in the day.  Good night
all.

 >
 > Hi David,
 >
 > Cheers for the code : )
 >
 > It's almost working, but I am getting some funny behaviour.
 >
 > If I use the first method, applying forces to each face normal, I get
 > very  strange results.
 >
 > Say I make a Box 200L by 200W by 10H, set the viscosity quite high, put
 > it  in my world with gravity set. The box will fall under gravity and be
 > slowed  by the hydrodynamic forces, great!! If I rotate the box by 90
 > degrees on its  Z axis, then the box experiences less resistance,
 > magic!!!
 >
 > But if I make the Box 10L by 200W by 200H, put in my world rotated so a
 > large face is pointing down, it will fall as if the edge is pointing
 > down.  If I rotate it so the edge is pointing down it will then act as
 > if one of  it's large face's is pointing down. You get the picture.
 >
 > Also, if I rotate the first box by 45 degrees it will fall slowly. But,
 > instead of moving in the direction of its lowest edge it will move to
 > the  opposite side.
 >
 > The face normals seem to be getting mixed up!
 >
 > I hope I am making sense.
 >
 > If I use the second method it explodes!
 >
 > I am going to spend more time tonight testing to see if I can work out
 > what  is going wrong .
 >
 > Cheers,
 >
 > Ander
 >
 >
 > David Wrote:
 >
 > I wrote up a fairly lengthy discussion of the concept of finding the
 > area of the 'silhouette' or 2 dimensional cross-sectional area of the
 > basic geoms as viewed from the direction of their velocity in
 > http://q12.org/pipermail/ode/2003-April/004101.html .
 >
 > To go from your direction, i.e. to get the velocity that each side of
 > the cube is going in the direction of it's normal, you can use a fairly
 > similar concept (which may actually turn out to be faster).
 >
 > First, how to get the normals for the faces of the cube:  Note that the
 > normals for the faces are the bodies x, y, and z axes.  Therefore, the
 > normal for those faces are the first, second, and third row of the R
 > matrix.  They will already be normalized and in global coordinates, so
 > that's a couple of computations knocked off my method right there.
 >
 > Now you need the magnitude of the projection of the velocity in the
 > direction parallel to the x, y, and z axes.  This is nothing more than a
 > dot product.  Dot your velocity with each of the x, y, and z axes and
 > you have the magnitude of the velocity (i.e. the speed) in those
 > directions.
 >
 > For the box:
 >
 >      ____
 >     /    /
 >    /____/| h
 >   |    | /
 >   |____|/ l
 >     w
 >
 >   z  y
 >   | /
 >   |/___x
 >
 > Which I beleive is the way ODE's box sides are laid out, your final
 > calculation would be:
 >
 > const dReal *lvel = dBodyGetLinearVelocity(body);
 > const dReal *R = dBodyGetRotation(body);
 > dReal ss[3];
 > dGeomBoxGetLengths(geom, ss);
 > dReal l = ss[0], w = ss[1], h = ss[2];
 > dReal side = l * h * (R[0]*lvel[0] + R[1]*lvel[1] + R[2]*lvel[2]); dReal
 > front = w * h * (R[4]*lvel[0] + R[5]*lvel[1] + R[6]*lvel[2]); dReal top
 > = w * l * (R[8]*lvel[0] + R[9]*lvel[1] + R[10]*lvel[2]);
 >
 > //either this......
 > dReal temp = -side*scale;
 > dBodyAddForce(body, temp*R[0], temp*R[1], temp*R[2]);
 > temp = -front*scale;
 > dBodyAddForce(body, temp*R[4], temp*R[5], temp*R[6]);
 > temp = -top*scale;
 > dBodyAddForce(body, temp*R[8], temp*R[9], temp*R[10]);
 >
 > //or this is probably close enough, and will save a few clocks:
 > dReal temp = -scale*(side+front+top);
 > dBodyAddForce(body, temp*lvel[0], temp*lvel[1], temp*lvel[2]);
 >
 > The scale factor represents the viscosity of the fluid you are moving
 > through, so this same method could be used to simulate air with a low
 > scale or water with a higher scale.  Man, I wish I had thought of this
 > method before I wrote that one up.  Hope it works for you.
 >
 > David
 >
 >
 > _________________________________________________________________
 > Hotmail now available on Australian mobile phones. Go to
 > http://ninemsn.com.au/mobilecentral/hotmail_mobile.asp
 >
 > _______________________________________________
 > ODE mailing list
 > ODE@q12.org
 > http://q12.org/mailman/listinfo/ode



_______________________________________________
ODE mailing list
ODE@q12.org
http://q12.org/mailman/listinfo/ode


_________________________________________________________________
Hotmail now available on Australian mobile phones. Go to  
http://ninemsn.com.au/mobilecentral/hotmail_mobile.asp