SubDyn: Consistent postprocessing for self-weight reconstruction

[RBergua]

This PR is not ready to be merged. It's a work in progress.

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

Related issue, if one exists
#2325

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

However, all r-test that use SubDyn and request beam outputs 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. This test will consist in a system made of 3 beams originally vertical and rotated in time domain 90 degrees through the interface joint to ensure the self-weight consistency.

Below I have rerun the systems presented in the issue commented above with the proposed changes.

Test 1: Vertical cantilever beam

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

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.

Next I will check:

1. Loading in floating systems (rigid body motion).