[ODE] Approximate Hydrodynamic Forces

[Ander Taylor]

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