Body kinematics when origin not at center of mass

[AlexWKinley]

When the origin of a body is not at the center of mass (including any attached objects), the kinematics of bodies are not correct.

As a somewhat extreme demonstration of this, we can set up two kinematically identical situations, one where the body is at the center of mass, and one where the body is far from the COM.

--------------------- MoorDyn Input File -------------------------------------------------------
Testing body kinematics when body origin != COM
---------------------- ROD TYPES ------------------------------------
TypeName      Diam     Mass/m    Cd     Ca      CdEnd    CaEnd
(name)        (m)      (kg/m)    (-)    (-)     (-)      (-)
buoy          1       1.0e1     1.2    0.0     1.2      0.0
---------------------------- BODIES -----------------------------------------------------
ID   Attachment  X0     Y0    Z0     r0      p0     y0     Mass  CG*   I*      Volume   CdA*   Ca
(#)     (-)      (m)    (m)   (m)   (deg)   (deg)  (deg)   (kg)  (m)  (kg-m^2)  (m^3)   (m^2)  (-)
1       free      0     0      20      0.0       0      0      0     0    0        0       0      0
2       free      0     10      20     0.0       0      0      0     0    0        0       0      0
---------------------- RODS ----------------------------------------
ID     RodType   Attachment  Xa    Ya    Za    Xb    Yb    Zb    NumSegs   RodOutputs
(#)    (name)    (#/key)     (m)   (m)   (m)   (m)   (m)   (m)   (-)       (-)
1      buoy     Body1         0     0     -2     0   0.0    2     10        -
2      buoy     Body2         0     -10   -2     0  -10.0   2     10        -

In this case we have the same 4 meter tall rod, with it's center at (0, 0, 20).
Body 1 starts at (0, 0, 20), and body 2 starts 10 meters to the side in the +y direction (with the rod a corresponding -10 meters to the side to be back at the origin).

In this simulation, gravity is set to 0, and we give the two bodies velocities equivalent to their rod rotating at 0.2 rad/sec around the x axis.
For body 1 this means an initial velocity vector of
0.0, 0.0, 0.0, 0.2, 0.0, 0.0
and for body 2 this means an initial velocity vector of
0.0, 0.0, 10 * 0.2, 0.2, 0.0, 0.0.
The 10*0.2 m/s velocity in the +z direction is due to the fact that the body is 10 meters from the center of rotation.

Running this simulation and viewing the results in paraview this is what we see:

body_center_of_mass.mp4

This is obviously incorrect, since both situations should be equivalent in terms of rod position/rotation.

In terms of solutions, I don't have much beyond avoiding this sort of scenario. I'm pretty certain that the problem is related to missing non-inertial force on the body (likely at least the centrifugal force). But I don't have a good grasp on how that would need to be implemented.
I will note that in more practical situations, things like fluid drag or buoyancy forces often do a pretty good job of mitigating most of the error caused by the incorrect kinematics.

Additional complications include

• Added mass turns masses from scalars to matrices makes deriving the necessary equations more complicated.
• I've found very little literature on handling kinematics of an object in a global coordinate frame when the object origin is not at the center of mass. Most things I've found either uses a body fixed coordinate frame, or require the object origin to be at the center of mass.

I don't necessarily have a strong feeling about how much effort should be put into addressing this issue. If there was sufficient reference theory and desire from other, I likely could put a fair amount of time into this. But I figured I'd bring this up to determine if it was a known limitation, and if there had been any other work (perhaps in OpenFast) to handle issues like this.

[RyanDavies19]

@AlexWKinley does this also happen if you define the body's center of gravity in the input file (should be defined in the bodies reference frame)?

[AlexWKinley]

Assuming that this would be the correct way to do that

---------------------------- BODIES -----------------------------------------------------
ID      Attachment  X0     Y0      Z0      r0        p0     y0      Mass  CG*        I*        Volume   CdA*   Ca
(#)     (-)         (m)    (m)     (m)     (deg)     (deg)  (deg)   (kg)  (m)        (kg-m^2)  (m^3)    (m^2)  (-)
1       free        0      0       20      0.0       0      0       0     0          0         0        0      0
2       free        0      10      20      0.0       0      0       0     0|-10|0    0         0        0      0

then it does still happen with the CG defined at the rods center of gravity. Which I think makes sense since the body cg is only used to account for the body's own mass, which is zero in this situation.

Edit:
If you want to experiment with this more, my full setup is
You'll need to have VTK enabled, but you can also strip out all the vtk parts if needed.
If the body deviates from the correct path (a circle around the rod), the test prints out the error in the position.

tests/Mooring/body_tests/floating_bodies.txt
tests/bodies.zip

[sanguinariojoe]

I am dealing with a model now which is failing because of this. I think I will spend a bit of time fixing it, so I assign this bug to myself.

If anyone has already done some advances on this topic, please let me know

[AlexWKinley]

@sanguinariojoe

I haven't had much time or pressing need to look into this, so I haven't done much work. I'll lay out what thinking I have done, and you're free to work off of that however.

So far as I can tell, the type of solution most likely to successful is to determine the center of mass of a body along with any of the connected objects, and compute the kinematics of that center of mass, then from there compute the kinematics of each individual object. That approach ideally avoids much of the complexity related to determining kinematics based on a point not at the center of mass, potentially at the cost of more transformations between reference frames.

What I don't know how to handle in that case is the inertial component of added mass, especially for partially submerged objects. In particular the question of where the center of mass is seems like it would be complicated by mass being a rank 2 tensor (a matrix) as opposed to a scalar. I've done some looking but haven't managed to find details as to how anyone handles this. I think it's likely some sort of solution exists, but I don't have a good enough understanding of the math/theory to confidently come up with something. It's also possible that a reasonable approximation can be accomplished by ignoring the effect of added mass on the center of mass, but I would worry that such an assumption would break down in some cases.

If you are going to be working on this, it would be great to lay out the theory of whatever approach you want to use, before too much time gets spent on writing/reviewing code. Let me know if there's any way I could be of assistance, I'm always happy to answer questions or discuss ideas.

[sanguinariojoe]

Well, for the time being I found a bug. I have also added your test @AlexWKinley.

[sanguinariojoe]

OK, I have been checking out a bit and I am pretty sure the problem is that we are neglecting the centrifugal forces. Effectively, as long as the body is rotating the rod shall exert a force towards itself.

That's why on the demo video above, the body is just rotating and moving upwards while the rod "orbits" around it, without doing any effect.

Do we agree @AlexWKinley, @RyanDavies19 and @mattEhall ?

If we accept the diagnosis, I suppose the smartest way to go is adding new methods to the Rod:: and Point:: APIs to compute centrifugal forces:

vec Rod::getCentrifugalForce(vec rBody, vec angVel) const;
inline vec Rod::getCentrifugalForce(vec rBody, vec6 vBody) const { return getCentrifugalForce(rBody, vBody.tail<3>()); }

Note: Those functions shall be called after ::getNetForceAndMass(), so the mass matrices are known.

What do you think guys?

[AlexWKinley]

I think conceptually this is on the right path. I'm not sure that the force we're thinking about is technically the centrifugal force, given that we're not really using a rotating reference frame at any point. But the math for both of them are basically the same, and I don't have a better term.

One thing to be weary of is that with the way MoorDyn is setup I don't think we can consider this "force" in the same way we consider other forces. Because MoorDyn bodies use a 6x6 mass matrix, I don't think treating the centrifugal force as an ordinary external force works. The bottom left quadrant of the mass matrix create a coupling between linear forces and angular acceleration, which I don't believe is something that applies to the centrifugal force.

In my provided example, the 6x6 mass matrix of the body that is displaced from the rod, at the very first timestep is:

$$\begin{bmatrix} 40 & 0 & 0 & 0 & 0 & 400 \\\ 0 & 40 & 0 & 0 & 0 & 0 \\\ 0 & 0 & 40 & -400 & 0 & 0 \\\ 0 & 0 & -400 & 4055.83 & 0 & 0 \\\ 0 & 0 & 0 & 0 & 55.8333 & 0 \\\ 400 & 0 & 0 & 0 & 0 & 4005 \end{bmatrix}$$

The centrifugal force is m * r * w^2 = 40 * 10 * 0.2^2 = 16N in the -x direction.

We expect the acceleration to be 16N / 40kg = 0.4m/s^2 in the -x direction. Since that is what would correct our kinematics due to not being at the center of mass.
But if we were to see what force is required to achieve that acceleration (F = ma) we get

$$\begin{bmatrix} 40 & 0 & 0 & 0 & 0 & 400 \\\ 0 & 40 & 0 & 0 & 0 & 0 \\\ 0 & 0 & 40 & -400 & 0 & 0 \\\ 0 & 0 & -400 & 4055.83 & 0 & 0 \\\ 0 & 0 & 0 & 0 & 55.8333 & 0 \\\ 400 & 0 & 0 & 0 & 0 & 4005 \end{bmatrix} \begin{bmatrix} -0.4 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} -16 \\ 0 \\ 0 \\ 0 \\ 0 \\ -160 \end{bmatrix}$$

which is a different force than what our centrifugal force would give.
Perhaps there is a straightforward way to account for this when computing our "centrifugal force", but I haven't dived into it enough to be sure, especially in cases where you have multiple objects contributing mass.

The other thing is that this approach still potentially needs to contend with the idea of determining center of mass including hydrodynamic added mass. In cases where a rod is partially submerged, the part that is under water will effectively have a greater mass, shifting the center of mass, at least for certain axes of rotation. I don't know how important this effect actually is, and the idea of added mass is already a simplification that may be inaccurate for rods large enough for this effect to matter. For my use cases, I'm generally using rods as surface piercing buoys, where the body is at the attachment point of the bouy (to make getting the position information I want slightly easier). So I may have more concern for behavior of rods in surface piercing than others do.

But the main thing of concern is reconciling the 6x6 mass matrix with this "centrifugal force". If I get more time in the near future, I can try and work out more of the theory. But until then, I'm happy to continue providing thoughts and feedback. And if you've got any code you want me to test, I can do that too.

[RBergua]

This is a very interesting topic! As part of the USFLOWT project, we use very large buoyant cans that we model in MoorDyn. You can find image of the OpenFAST model below for reference:

As part of this modeling effort, we performed some verification cases at the beginning of the project. The buoyant cans are able to rotate around the leg tip (we use CoupledPinned bodies in MoorDyn).

The definition that we used for this verification case is attached here:
MoorDyn_v2_one_can_test3.txt

We performed a code-to-code verification for the free-decay motion of these cans (the leg tip location from SubDyn is fixed in place, effectively enabling only rotations of the cans by means of the CoupledPinned approach and imposing an initial angle of 10 degrees). When accounting for the hydrostatics and the drag contribution of the cans, the agreement was relatively good. However, when we accounted for the added mass, clear differences were observed. See below for reference:

I'm concerned about the rod being partially submerged and how the added mass is computed. Are you updating the mass matrix at every time step according to the submerged portion of the rod?

Thanks for the support!

[RyanDavies19]

@AlexWKinley @sanguinariojoe I will add that kinematics between MD-F and MD-C have been identical in all verification testing I have done so my thought was this but could be causing the USFLOWT bugs as well.

[sanguinariojoe]

Tomorrow I am giving it another thought. Anyway, adding centrifugal/centripetal forces is the common way to go in multibody dynamics, afaik. See for instance: https://fab.cba.mit.edu/classes/865.18/design/optimization/dynamics_1.pdf I really think we are missing the second term in eq. (24). I first thought that it is possible to just simply add it later as a force and let getderivstate() do the job for us. But I think Alex is right and we need to do that on a different way

…

On Mon, 1 Jul 2024, 19:54 RyanDavies19, ***@***.***> wrote: @AlexWKinley <https://github.com/AlexWKinley> @sanguinariojoe <https://github.com/sanguinariojoe> I will add that kinematics between MD-F and MD-C have been identical in all verification testing I have done so my thought was this but could be causing the USFLOWT bugs as well. — Reply to this email directly, view it on GitHub <#207 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMXKKC3R7C5RDK5CGKBLXLZKGJV7AVCNFSM6AAAAABGH4FSNWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMBQG4YTMMBSGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[AlexWKinley]

@sanguinariojoe

While that linked document does seem like a high quality resource on dynamics, and probably have some value. I will note that the term you see there in the Newton-Euler equations is not actually present in our situations.

That term is part of Euler's equations, and serves to handle the differences between a global and rotating reference frame (see the derivation section for an explanation).

Because MoorDyn uses a global reference frame for everything, we don't need such a term. But we do need to somehow compensate for the fact our positions are not at the center of mass/rotation.

[sanguinariojoe]

In my provided example, the 6x6 mass matrix of the body that is displaced from the rod, at the very first timestep is:

The centrifugal force is m * r * w^2 = 40 * 10 * 0.2^2 = 16N in the -x direction.

We expect the acceleration to be 16N / 40kg = 0.4m/s^2 in the -x direction. Since that is what would correct our kinematics due to not being at the center of mass. But if we were to see what force is required to achieve that acceleration (F = ma) we get

Nope, you have a typo @AlexWKinley ... Both go on the -y direction... And thus the force resulting from the acceleration is the centripetal force.

Let me try to get the full picture considering a point instead of a rod, for the sake of simplicity:

There are no added masses nor drag froces, and the body, b, is massless, having all the mass, m, concentrated on the point p.

The kinematics say that

The mass matrix (3x3) of the point is just simply

which let's us assemble the mass matrix on the body

with H the alternator matrix of q. Then the Newton equation is

The upper part can be expanded right away as the centripetal force

The down part we can rearrange it as

However, by definition,

So, as expected:

---

I am making some tests and let you know guys.

---

I think we should open a parallel issue for the @RBergua case, so we focus on the @AlexWKinley case now and move to that afterwards

[sanguinariojoe]

OK, I added centripetal forces to the computation of bodies with attached rods and points here, and modified the test here.

When orbiting around the rod it is nailing the kinematics, while when orbiting around the point it tends to accelerate (at least during the beginning of the simulation). The orbital distance is well captured in all the cases.

Now I wonder... What about the attached lines?

---

P.S. I added a body with an attached line like this:

--------------------- MoorDyn Input File -------------------------------------------------------
Testing body kinematics when body origin != COM
----------------------- LINE TYPES ------------------------------------------
TypeName   Diam    Mass/m     EA         BA/-zeta    EI         Cd     Ca     CdAx    CaAx
(name)     (m)     (kg/m)     (N)        (N-s/-)     (N-m^2)    (-)    (-)    (-)     (-)
main       1       1.0e1      0.0        -0.0        0          0.0    0.0    0.0     0.0
---------------------- ROD TYPES ------------------------------------
TypeName      Diam     Mass/m    Cd     Ca      CdEnd    CaEnd
(name)        (m)      (kg/m)    (-)    (-)     (-)      (-)
buoy          1        1.0e1     0.0    0.0     0.0      0.0
---------------------------- BODIES -----------------------------------------------------
ID      Attachment  X0     Y0      Z0      r0        p0     y0      Mass  CG*        I*        Volume   CdA*   Ca
(#)     (-)         (m)    (m)     (m)     (deg)     (deg)  (deg)   (kg)  (m)        (kg-m^2)  (m^3)    (m^2)  (-)
1       free        0      10      20      0.0       0      0       0     0|-10|0    0         0        0      0
2       free        0      10      20      0.0       0      0       0     0|-10|0    0         0        0      0
3       free        0      10      20      0.0       0      0       0     0|-10|0    0         0        0      0
----------------------- POINTS ----------------------------------------------
Node      Type       X         Y          Z        M        V          CdA    CA
(-)       (-)        (m)       (m)        (m)      (kg)     (m^3)      (m^2)  (-)
1         Body1      0.0       -10        0        40.0     0          0      0
2         Free       0.0       -12        20       0.0      0          0      0
3         Body3      0.0       -8         0        0.0      0          0      0
---------------------- RODS ----------------------------------------
ID     RodType   Attachment  Xa    Ya    Za    Xb    Yb     Zb    NumSegs   RodOutputs
(#)    (name)    (#/key)     (m)   (m)   (m)   (m)   (m)    (m)   (-)       (-)
1      buoy      Body2       0     -10   -2    0     -10.0  2     10        -
---------------------- LINES ----------------------------------------
ID   LineType   AttachA  AttachB  UnstrLen  NumSegs  LineOutputs
(#)   (name)     (#)      (#)       (m)       (-)     (-)
1     main       2        3         4         10      -
-------------------------- SOLVER OPTIONS---------------------------------------------------
2        writeLog     - Write a log file
0.0      g            - No gravity
1e-3     dtM          - time step to use in mooring integration
3.0e6    kb           - bottom stiffness
3.0e5    cb           - bottom damping
70       WtrDpth      - water depth
3.0      ICDfac       - factor by which to scale drag coefficients during dynamic relaxation IC gen
0.015    threshIC     - threshold for IC convergence
0.0      TmaxIC       - threshold for IC convergence
0.01     dtIC         - Time lapse between convergence tests (s)
0        Currents     - Whether or not to pull in currents
0        WaveKin      - Whether or not to pull in waves
--------------------------- need this line -------------------------------------------------

and it worked reasonably OK. Probably a more formal analysis is worthy (and probably we are getting surprises), but for the time being I would say the patch is good enough, and nothing special shall be done for the attached lines. We have a technical debt

---

P.S.II Here I attached two points instead of one to the body and the orbital motion is way better

[sanguinariojoe]

Fixed on #231