[ODE] Approximate Hydrodynamic Forces

[David Whittaker]

Hmmm, just thought of an interface change that would make this a lot more
general.  If you took an array of planes (i.e. a, b, c, and d of
ax+by+cz=d) as one of the arguments, then an array of linear_viscosity,
angular viscosity, density, and flow, and the number of planes (other
arrays are num_planes+1).  Then each entry in the arrays could correspond
to the fluid parameters below the corresponding plane.  For example, if I
wanted to create a water level at z = 0 (z pointing up):

dReal planes[4] = {0,0,1,0};
//no clue if these parameters are "right", just an example
dReal lin_vis[2] = {10, 0.1};
dReal ang_vis[2] = {10, 0.1};
dReal density[2] = {1, 0.01};
Vector3 flow;
flow.x = flow.y = flow.z = 0;
Vector3 flows[2] = {flow, flow};

box->ApplyHydrodynamicForces(planes, lin_vis, ang_vis, density, flows, 1);

That would make your function start out looking like:

void BoxGeom::ApplyHydrodynamicForces(dReal planes[], dReal lin_vis[], dReal
ang_vis[], dReal densities[], Vector3 flows[], int num_planes)
{
  const dReal *pos = dBodyGetPosition(this->_body->Id);
  dReal linear_viscosity, angular_viscosity, density;
  Vector3 flow;
  int p;
  for (p = 0; p < num_planes; p++)
  {
    if (planes[p*4+0] * pos[0]
      + planes[p*4+1] * pos[1]
      + planes[p*4+2] * pos[2]
      < planes[p*4+3])
    {
      linear_viscosity = lin_vis[p];
      angular_viscosity = ang_vis[p];
      density = densities[p];
      flow = flows[p];
      break;
    }
  }
  if (p == num_planes)
  {
    linear_viscosity = lin_vis[p];
    angular_viscosity = ang_vis[p];
    density = densities[p];
    flow = flows[p];
  }

  ... and now, for the rest of the function ...
}

You could pass null as planes and 0 as num_planes, and the address of the
variables to get the same functionality as before (i.e. the whole world in
one fluid).

You might even be able to make a floating ship, though it would float at
it's midpoint, not necessarily it's correct height (until the partially
submerged dynamics get worked out), unless there were enough bodies
vertically that it was a good approximation.

Just a thought,
David

>>
>> 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
>
>
>
> _______________________________________________
> ODE mailing list
> ODE@q12.org
> http://q12.org/mailman/listinfo/ode