SubDyn: Consistent postprocessing for self-weight reconstruction

[RBergua]

This PR is ready to be reviewed.

Feature or improvement description
SubDyn computes the self-weight in beams using equivalent nodal forces and moments (see for reference: https://openfast.readthedocs.io/en/main/source/user/subdyn/theory.html#self-weight-loads).

Currently, SubDyn evaluates internal forces using only the nodal displacements:

However, the postprocessing should include a superposition correction for the equivalent nodal forces and moments. This would return consistent beam outputs.

Interestingly, when looking at the source code, it seems that this was indeed the behavior till OpenFAST v 2.6.0. The changes were introduced in commit 0631a9b. It seems that the reason for these changes were related to the pretensioned cables r-test: https://github.com/OpenFAST/r-test/tree/main/modules/subdyn/SD_Cable_5Joints

The above changes work as expected for fixed-bottom systems (where the self-weight forces and moments are constant). For floating systems, the self-weight components computed at initialization are projected into the current floating reference body before output. The self-weight bending moment components are recomputed from scratch using the current element direction cosine at every time step.

Related issue, if one exists
#2325

Impacted areas of the software
Note that this change only changes the sensor outputs, not the physics of the system (e.g., the eigenfrequencies or deflections of the system are not impacted).

However, all r-test that use SubDyn and request beam outputs will return different values. Previously, the self-weight was not properly accounted for.

Test results, if applicable
I'm thinking about including one new r-tests to ensure consistency in the SubDyn self-weight calculation and outputs: Test 4 detailed below.

Test 1: Vertical cantilever beam (from issue cited above)

The beam cross-sectional properties are as follows:
Hollow cylinder
D = 3 m
t = 0.1 m

Material properties:
E = 210E9 N/m^2
G = 8.077E10 N/m^2
ρ = 7,860 kg/m^3

The distributed weight from the free-end can be computed as follows:
m = ρ·A·g = 7,860·(π·((3/2)^2-((3-2·0.1)/2)^2))·9.80665 = 70.225 kN/m

As commented, this fix only impacts the output sensors. Not the physics of the system. For example, the first bending mode (3.243 Hz) and the second bending moe (18.59 Hz) remain the same.

Test 2: Horizontal cantilever beam (from issue cited above)

The beam cross-sectional properties are as follows:
Hollow cylinder
D = 5 m
t = 0.1 m

Material properties:
E = 210E9 N/m^2
G = 8.077E10 N/m^2
ρ = 7,860 kg/m^3

The uniformly distributed load due to the gravity acceleration can be computed as:
ω = ρ·A·g = 7860·(π·((5/2)^2-((5/2)-0.1)^2))·9.80665 = 118.656 kN/m
The shear force can be computed analytically as: F = ω·L
The bending moment can be computed analytically as: M = (ω·L^2)/2

Test 3: Pretensioned cables (from r-test, commented above)
Schematic representation of the system:

A displacement along the x direction is prescribed to the interface joint:

Accordingly, the pretensioned cable at the left side increases the tension and the pretensioned cable at the right side drops the tension. The axial stiffness of the cable at the left side is half compared to the right side:

As observed, the results are consistant.

Test 4: System with two members (4 beam elements, NDiv = 2) and one concentrated mass with an eccentricity that experiences rigid body motion [floating systems] (this is being added as new r-test for self-weight verification: OpenFAST/r-test#182)

In this case, the angular velocity of the interface joint around the global y-axis is -0.0175 rad/s (-1 deg/s). This results in the -90 deg rotation shown above in a time span of 90 s. The angular velocity is constant, so there is no angular acceleration.

The beam cross-sectional properties are as follows:
Hollow cylinder
D = 5 m
t = 0.1 m

The concentrated mass m = 2,000 kg has an eccentricity of 1 m from the free-end of the beam.

Material properties:
E = 210E9 N/m^2
G = 8.077E10 N/m^2
ρ = 7,860 kg/m^3

The uniformly distributed load due to the gravity acceleration can be computed as:
ω = ρ·A·g = 7860·(π·((5/2)^2-((5/2)-0.1)^2))·9.80665 = 118.656 kN/m

The weight due to the concentrated mass and the beam can be computed analytically as: W = 2000*9.80665+ω·(10-L).

This weight (always vertical due to the gravity acceleration) has to be projected to the local beam orientation at every time step. For example, at t = 0 s (horizontal beam), the load is a pure shear force while at t = 90 s (vertical beam), the load is a pure axial force.

In general terms, the shear and axial forces can be computed according to the orientation angle (α) as follows:
Shear Force = W·cos(α)
Axial Force = W·sin(α)

The bending moments can be computed as:
Bending Moment = 2000·9.80665·(10+1-L)·cos(α) + w·(10-L)·(((10-L)/2).cos(α))

For reference, if the concentrated mass is removed, the loads at the beam free-end (M2N3) are zero as expected.