[ODE] Approximate Hydrodynamic Forces

Tue Apr 15 22:55:02 2003

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