[ODE] what a drag

[david@csworkbench.com]

Ok, I've had a chance to think this through, and here's what I came up with.

First, I don't think the method at the link below will work because it
calculates the cross sectional area of an Axis Aligned Bounding Box.  If
you tried to create an AABB aligned to the velocity every step, you would
spend more time doing that than this method.  It also assumes a fixed
viewing point, not a direction.  The difference in area returned would be
similar to the difference in the amount of stuff lit when you use point
lights as opposed to directional lights in OpenGL.

======================================================================

The first step in a few of the procedures is to translate the velocity
vector into the body frame of reference and normalize it (you can either
translate the velocity vector into body coordinates or all the body
vectors into world coordinates..... this way should be faster):

dReal *lvel = dBodyGetLinearVel(body);
dVector3 blvel;
dBodyVectorFromWorld(body, lvel[0], lvel[1], lvel[2], blvel);
dNormalize4(blvel);

Normalizing the vector saves you several sqrt's later on.

=======================================================================

For a box:
    ____
   /    /
  /____/| h
 |    | /
 |____|/ l
   w

 z  y
 | /
 |/___x

Now that we're on the same terms with dimensions and axes (excuse my poor
ascii art, but you get the idea), I can explain a little better.  You have
3 dimensions: l, w, and h, that in combinations make up the surface area
of a given face.  Now, take each of those dimensions and turn them into
vectors:

w' = [w,0,0], l' = [0,l,0], h' = [0,0,h].

Now we need the projection of each vector in the direction perpindicular
to blvel.  The equation for that is:

p = v - dot(v,u) * u

where u is a unit vector in the direction of the velocity, v is one of the
dimension vectors, and p is the vector projected perpindicular to u.
let's see what this would come out to be with the width vector (sw=scaled
width, [x,y,z] = linear vel in body coordinates (blvel)):

sw = [w,0,0] - dot([w,0,0],[x,y,z]) * [x,y,z]
sw = [w,0,0] - w*x * [x,y,z]
sw = [w-w*x*x, -w*x*y, -w*x*z]

and the others:

sl = [-l*x*y, l-l*y*y, -l*y*z]
sh = [-h*x*z, -h*y*z, h-h*z*z]

area of the 2d parallelogram projected to a plane perpendicular to the
velocity is now the magnitude of the cross product of the two vectors that
make up the face.

sw x sh = (some algebra) = w*h*[x*y, 1-x^2-z^2, y*z]

little more algebra....

afront = ||sw x sh|| = w*h*sqrt(y^2*(x^2+z^2)+(x^2+z^2-1)^2)

and the others:

atop   = ||sw x sl|| = w*l*sqrt(z^2*(x^2+y^2)+(x^2+y^2-1)^2)
aside  = ||sl x sh|| = l*h*sqrt(x^2*(y^2+z^2)+(y^2+z^2-1)^2)

Note that the above formulas follow a nice pattern, so the benefit from
coding them out as these formulas (as opposed to using dCROSS and taking
the magnitude) is that a lot of repetitve multiplies get cancelled out. 
You can save x^2, y^2, z^2, and the sum of each pair, and just sub in the
variables where needed above.  This should save lots of multiplies, so it
will speed things up considerably.

Finally, the total cross sectional area, i.e. the area of the silhouette
as viewed from the direction of motion, is:

atotal = afront + atop + aside

(After you Normalize blvel, the rest is just derivation of the final
equations.  x is blvel[0], y is blvel[1], and z is blvel[2].... w, l, and
h are the lengths of the sides of the box.  Just sub in to the equations
for afront, atop, and aside, sum them up, and you're done).

====================================================================

A sphere is simple.  It's cross sectional area is always pi*r^2.

atotal = pi*r^2

====================================================================

A capped cylinder is only slightly more complicated.  You can take the
caps out of the equation by realizing that the will always supply pi*r^2
to the total cross sectional area.... i.e. if the capped cylinder was
length 0, it would be a sphere (assuming the caps aren't counted in the
length).  Now we just have the l x 2r rectangle (the main part of the
body) left.  You may say it's not a rectangle once it's rotated away from
directly looking at it, but you can cut the curves off one end and put
them on the other and it becomes a rectangle again....  I wish I could
draw it, but I can't do curves in ASCII....  It's the same argument that
says a parallelogram's area is base x height, even if it's tilted.  So now
that we've established the area is length * 2 *radius, let me say that the
2 * radius part doesn't change, no matter what direction you are looking
from, so the only thing left is to find the scaled length.  Again start by
translating the velocity into body coordinates, as we did in the box. 
Then, borrowing an intermediate step from the box:

sh = [-h*x*z, -h*y*z, h-h*z*z]

since capped cylinders are aligned to the z axis, I chose the scaled
height vector.  Let's change that to length to avoid confusion:

slength = [-length*x*z, -length*y*z, length-length*z*z]

Now all we need is the scalar magnitude of this vector:

sslength = sqrt(slength[0]^2 + slength[1]^2 + slength[2]^2)

and the formula for the cross sectional area of a capped cylinder would be:

atotal = pi*r^2 + 2*r*sslength

====================================================================

For a non-capped cylinder, the main part of the body is still going to be
2*r*sslength as calculated above.  I'm almost positive that the circle on
the end would contribute abs(pi*r^2*dot(blvel,[0,0,1])), but I can't prove
it without thinking over it some more.

atotal = abs(pi*r^2*blvel[2]) + 2*r*sslength

====================================================================

For a triangle, you would take the same approach as the box, but more
generalized.  First, pick a vertex, and find the vectors in world
coordinates to the other two vertices.  Call these v1 and v2.  Then
normalize the linear velocity (in world coordinates this time).  Call this
u.

s1 = v1 - dot(v1,u) * u
s2 = v2 - dot(v2,u) * u

s1 and s2 are now your scaled lengths projected perpendicular to the
velocity.

atotal = ||s1 x s2|| / 2

Should give you the area scaled to the velocity direction.  Go ahead and
use a full dot and cross product here, there won't be any more
optimizations unless you have stringent triangle definitions. i.e. all
triangles have a vertex at the origin of their body and a vertex on the
body x axis.... not likely.

====================================================================

So that gives some fairly simple equations for the cross sectional /
silhouette area of the primitives I could think of.

Finally to get the fluid drag:

drag = - atotal * speed^2 * scalefactor;

where speed^2 = (lvel[0]^2 + lvel[1]^2 + lvel[2]^2)
and scalefactor is the coefficient of fluid resistance (you won't read
about that one in a book unless I got lucky, but hey, it sounds cool :)

Big scalefactor for water, small for air, etc.

Finally, assuming lvel is a copy of the linear velocity so you can change it:

dNormalize4(lvel);
dBodyAddForce(body, drag*lvel[0], drag*lvel[1], drag*lvel[2]);

Last but not least, this might be the most pencil and paper math I've done
since Junior High, so feel free to check up on my math and/or concepts. 
And of course, none of this is proven until it's implemented and looks
good.  If anyone wants to give it a try, go for it..... I'll get around to
it eventually if not, since I want water in my simulation, if I ever get
back to my simulation.

One parting word, the above ideas don't take into account objects which
are in 2 different fluids, i.e. a paddle.  If you use these ideas and want
to implement a paddle, make the handle a separate body from the blade, and
join them with a fixed joint.  Then if the center of the blade body is
below water level, use the water scale factor, otherwise use the air
scalefactor.  I'm not sure I have the mental energy or capability to find
a more accurate way to handle this case.  Maybe a well integrated verlet
integrator / particle engine could simulate water with "heavy" particles
(and not even need this drag code), but large bodies of water would be
difficult to do well and fast.

David

> On Fri, 4 Apr 2003, Pierre Terdiman wrote:
>
>> > The part about computing the cross section sounds hard, unless the
>> geoms are all spheres.  Anyone know of a simple and/or fast way to
>> do that for boxes at arbitrary orientations?
>>
>> Isn't it the same as in the old code from Dieter Schmalstieg, to
>> compute the screen projection area of an arbitrary box......... ?
>
> I'd never heard of him 'til now, but I think you're right.  If anyone
> else is interested, here's an URL:
>
> http://www.acm.org/jgt/papers/SchmalstiegTobler99/
>
> I just had an alternate idea though... Could one create "drag matrix"
> for a geometry object, such that multiplying a velocity vector by that
> matrix would result in an approximate (if not highly accurate) drag
> force vector?   I realize that I have just gone one step beyond my real
> understanding of linear algebra, so please be gentle if this is totally
> wrong. :-)
>
> --
>
> Nate Waddoups
> Redmond WA USA
> http://www.natew.com
>
>
> _______________________________________________
> ODE mailing list
> ODE@q12.org
> http://q12.org/mailman/listinfo/ode