Update to body drag coefficients

[RyanDavies19]

fixes #205 and closes #206.

Fixes the way that drag coefficients are handled so that the viscous drag force follows the following formula:

F \\
=
0.5 * \rho_w * \\
OrMat\\
\begin{bmatrix}
CdA_x & 0 & 0 \\
0 & CdA_y & 0 \\
0 & 0 & CdA_z  \\
\end{bmatrix}
OrMat^T *\\
(\vec v_i  * \\
|\vec v_i |)

The same equation is applied for the rotational drag coefficients.

[RyanDavies19]

@AlexWKinley can you look this over? My biggest question is the order of OrMat and the transpose correct, or should it be the other way (because we are going from the body frame to the global frame, and I believe OrMat transforms from global to body).

[AlexWKinley]

@RyanDavies19 As far as I can tell, this math is correct. After some thought (and relearning a little linear algebra) I'm pretty certain you're correct about the order of multiplication.

As a sort of sanity check (and sort of just for my own interest) I threw together this animation using Mathematica. Where we have a representative body and then an arrow representing the direction and relative magnitude of the drag force given a current in the +x direction. While it's a little hard to parse, having looked through it fairly carefully I feel pretty confident (acknowledging that my CdA values are very approximate).

https://github.com/FloatingArrayDesign/MoorDyn/assets/129765925/ddb3ecd2-f0d8-4620-af37-69657e3bdc96

Thanks for your and Matt's work on this, and for putting together this PR.

[RyanDavies19]

@AlexWKinley Great thank you for checking that! Nice animation by the way. I will merge this then

[AlexWKinley]

@RyanDavies19 Hate to be a bother, but I was doing some other work to do with drag, and realized that our new drag formulation is either defined imprecisely, or implemented slightly incorrectly.

In particular the interpretation of $\vec v_i |\vec v_i |$ I think is slightly incorrect. We currently have it implemented as vi.cwiseProduct(vi.cwiseAbs()).head<3>().
But I believe it's actually supposed to be `vi.head<3>() vi.head<3>().norm()`.

This is backed up by OrcaFlex's 6d buoy drag formulation https://www.orcina.com/webhelp/OrcaFlex/Content/html/6Dbuoytheory,Lumpedbuoyaddedmass,dampinganddrag.htm (equation 8).

The difference in these forms can be seen if we consider the case where vi = [1, 1, 0] and the rotation matrix is just the identity matrix. If Cdx and Cdy are equal, then we expect the drag to be a vector in the direction [1, 1, 0] with a magnitude (ignoring the 0.5 rho_w) of $Cdx \sqrt{1^2 + 1^2}^2 = Cdx 2$. But our current implementation gives Cdx * [1 * 1, 1 * 1, 0] which has a magnitude of $Cdx \sqrt{2}$.

If OrcaFlex's documentation is anything to go by, the same thing applies to the rotational degrees of freedom. But I'll admit I've given that case much less specific thought.

[RyanDavies19]

@AlexWKinley Good catch, from what I am seeing there on the OrcaFlex site rotational and translational drag are the same formulations. I left in the cwise approach because that had been how it was set up before, but to stay true to the equations we need to use the vector approach you laid out. I can make another PR with this change, and I will check the drag in the other objects too.

[AlexWKinley]

For the other objects that are cylindrical (lines and rods) I believe they intentionally deviate from that formulation slightly. In particular they separate relative velocity normal and relative velocity tangential. This is also something that OrcaFlex does (see standard formulation) https://www.orcina.com/webhelp/OrcaFlex/Content/html/Linetheory,Hydrodynamicandaerodynamicloads.htm

Thanks for your help with this.

[RyanDavies19]

@AlexWKinley They take a different approach and it aligns well with Orcaflex. I will also note that rods and lines both use the $v_i * v_i.norm()$ approach instead of component wise, which is more indication that was a bug in bodies.