[ODE] Approximate Hydrodynamic Forces

[David Whittaker]

>
> Hey David,
>
> This is what I got now. I works great.
>
> Could you have a look and let me know if you I should be doing anything
> different?
>
> Cheers,
>
> Ander
>
> note: I have added this to the .Net wrapper that I am using
>
> void BoxGeom::ApplyHydrodynamicForces(dReal linear_viscosity, dReal
> angular_viscosity, dReal density, Vector3 flow)
>     {
>         const dReal *lvel = dBodyGetLinearVel(this->_body->Id);
>         const dReal *avel = dBodyGetAngularVel(this->_body->Id);
>         const dReal *R = dBodyGetRotation(this->_body->Id);
>         dReal ss[3];
>         dGeomBoxGetLengths(this->_id, ss);
>
>         dReal AreaX = ss[1] * ss[2];
>         dReal AreaY = ss[0] * ss[2];
>         dReal AreaZ = ss[0] * ss[1];

You might consider making another dVector3 here that holds the compensated
velocity:
dVector3 compFlow;
compFlow[0] = lvel[0] - flow.x;
compFlow[1] = lvel[1] - flow.y;
compFlow[2] = lvel[2] - flow.z;

Use the stored values instead, and it'll save you 6 subtractions anyway
(and will probably exist completely in the fp stack).

>         dReal nx = (R[0] * (lvel[0] - flow.x) + R[4] * (lvel[1] -
> flow.y) +
> R[8] * (lvel[2]) - flow.z) * AreaX;

check your parenthesis here (lvel[2] - flow.z))

>         dReal ny = (R[1] * (lvel[0] - flow.x) + R[5] * (lvel[1] -
> flow.y) +
> R[9] * (lvel[2] - flow.z)) * AreaY;
>         dReal nz = (R[2] * (lvel[0] - flow.x) + R[6] * (lvel[1] -
> flow.y) +
> R[10] * (lvel[2] - flow.z)) * AreaZ;
>
>         dReal temp = -nx*linear_viscosity;
>         dBodyAddForce(this->_body->Id, temp*R[0], temp*R[4], temp*R[8]);
>
>         temp = -ny*linear_viscosity;
>         dBodyAddForce(this->_body->Id, temp*R[1], temp*R[5], temp*R[9]);
>
>         temp = -nz*linear_viscosity;
>         dBodyAddForce(this->_body->Id, temp*R[2], temp*R[6],
> temp*R[10]);
>
>         nx = (R[0] * avel[0] + R[4] * avel[1] + R[8] * avel[2]) * AreaZ;
> ny = (R[1] * avel[0] + R[5] * avel[1] + R[9] * avel[2]) * AreaX;
> nz = (R[2] * avel[0] + R[6] * avel[1] + R[10] * avel[2]) *
> AreaY;

Take a body that is rotating exactly around it's x axis, with a large
AreaY and a tiny AreaZ (i.e. A 1x1x100 box).  This box won't have as large
an nx as it should.  Perhaps one of the following might work better:

nx = dot(x, avel) * max(AreaY, AreaZ)
ny = dot(y, avel) * max(AreaX, AreaZ)
nz = dot(z, avel) * max(AreaX, AreaY)

or

nx = dot(x, avel) * (AreaY + AreaZ)
ny = dot(y, avel) * (AreaX + AreaZ)
nz = dot(z, avel) * (AreaX + AreaY)

or, probably most physical, but pretty expensive too:

dReal AreaDiagX = ss[0] * dSqrt(ss[1] * ss[1] + ss[2] * ss[2]);
dReal AreaDiagY = ss[1] * dSqrt(ss[0] * ss[0] + ss[2] * ss[2]);
dReal AreaDiagZ = ss[2] * dSqrt(ss[0] * ss[0] + ss[1] * ss[1]);

nx = dot(x, avel) * AreaDiagX;
ny = dot(y, avel) * AreaDiagY;
nz = dot(z, avel) * AreaDiagZ;

(excuse my dot shorthand... easier than typing out the first half again,
but the same).

>         temp = -nx * angular_viscosity;
>         dBodyAddTorque(this->_body->Id, temp * R[0], temp * R[4], temp *
>
> R[8]);
>
>         temp = -ny * angular_viscosity;
>         dBodyAddTorque(this->_body->Id, temp * R[1], temp * R[5], temp *
>
> R[9]);
>
>         temp = -nz * angular_viscosity;
>         dBodyAddTorque(this->_body->Id, temp * R[2], temp * R[6], temp *
>
> R[10]);
>
>         dMass m;
>         dBodyGetMass(this->_body->Id, &m);

When I wrote this up, I was thinking you would use the mass of the body
somewhere in the equation, did a little mid-email research, and determined
you wouldn't need it, but forgot to take it out.  The previous two lines
are useless.

>         dReal gravity[3];
>         dWorldGetGravity(this->_body->get_WorldId(), gravity);
>
>         temp = -density * ss[0] * ss[1] * ss[2];
>         dBodyAddForce(this->_body->Id, temp*gravity[0], temp*gravity[1],
>
> temp*gravity[2]);
>
>     }
>
> }
>
> _________________________________________________________________
> MSN Instant Messenger now available on Australian mobile phones. Go to
> http://ninemsn.com.au/mobilecentral/hotmail_messenger.asp