SubDyn output improvements: beam self-weight consistency, end-node forces, and revolute joint fixes

[RBergua]

This PR is ready to merge.

Feature or improvement description
This PR addresses several issues observed in SubDyn:

1. Consistent postprocessing for self-weight reconstruction in beams (@RBergua).
2. Force outputs at member end nodes/joints (@luwang00).
3. Loads and rotations outputs for revolute joints (@luwang00).
4. Deprecated inertial/dynamic load outputs (@luwang00).

The changes only affect the timeseries printout from SubDyn and do not impact SubDyn's internal states or module inputs/outputs. So, the physics of the system (e.g., eigenfrequencies, forces, deflections) are unchanged.

1. Consistent postprocessing for self-weight reconstruction in beams
SubDyn computes beam self-weight using a consistent load vector (i.e., equivalent nodal forces and moments). See the documentation for details: https://openfast.readthedocs.io/en/main/source/user/subdyn/theory.html#self-weight-loads

At present, internal forces are evaluated solely from the nodal displacements:

As a result, the contribution of the equivalent nodal forces and moments from the consistent load vector is not reflected in the reported forces. To obtain consistent beam outputs, the postprocessing step should include a superposition correction that reconstructs the effect of these loads:

This correction applies to any load that enters the system as a distributed load represented through a consistent element load vector during assembly of the global force vector F. Loads applied directly as nodal point forces do not require correction.

Interestingly, the source code shows that this was the behavior up to OpenFAST v2.6.0 for fixed-bottom systems. The current implementation appears to have been 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

This PR extends the fix to floating systems as well. The self-weight bending-moment components are fully recomputed at every time step using the actual rigid-body rotation matrix and element direction-cosine matrix.

Note that the equivalent force correction for self-weight in fixed-bottom and floating systems is not applied. Force reconstruction is instead carried out with the change described in Item 2 below.

2. Improved force outputs at member end nodes/joints
Currently, axial and shear forces at interior nodes of a member are averaged/interpolated from the forces over the two neighboring elements. This generally provides accurate results. However, at member end nodes, we only have a single element connecting to the node, and it is not possible to perform averaging. This can lead to large errors in the reported forces if the member discretization is coarse.

To address this, forces at end nodes are now computed using linear extrapolation from the forces over the first and second elements from the member end, providing significantly more accurate force estimates at the end nodes in most cases. If NDiv=1 (i.e., each member only has a single element), no extrapolation is performed. Therefore, it is recommended to use NDiv>=2 for improved end-node force accuracy.

This change affects both user-requested load outputs and the outputs obtained with the AllOuts option. Seabed reaction load outputs are not modified.

Effectively, we are assuming the internal force over the two elements at the end of a member to be continuous and vary linearly. Correspondingly, the applied force on the two elements is approximated as a distributed and piecewise constant load instead of lumped forces at the nodes. The existing SubDyn formulation for the forces at member interior nodes based on averaging the forces of the two neighbor elements implicitly relies on the same assumption/approximation.

Note that this method for estimating end-node forces can result in significant errors if there is a large lumped force applied at the first interior node next to the end node, such as mooring fairlead forces. This lumped force creates a jump in the internal force at the second node, causing the extrapolation to overshoot. This situation should be avoided when setting up the model. It is recommended to always map loads and motion between mooring fairleads and SubDyn user-defined joints instead of interior nodes.

3. Loads and rotation outputs for revolute joints
Previously, SubDyn was reporting incorrect member-end rotation and loads at revolute joints. This is because revolute joints have more than 6 DOF with 3 translational DOF and 3N rotational DOF, where N is the number of members connecting to this joint. For example, a revolute joint connecting two members will have 9 DOF total, with 3 translational DOF, 3 rotational DOF for the end of one member, and another 3 rotational DOF for the end of the second member. Previously, SubDyn was reconstructing and reporting member-end rotation and loads at revolute joints based on the first 3 rotational DOF always. As a result, the reported values were incorrect for the second member and onward. This bug is now fixed by having SubDyn looking up the appropriate set of rotational DOF based on the member being queried.

4. Deprecated inertial/dynamic load outputs
The previously available inertial/dynamic load outputs (e.g., M1N1FMxe and M1N2MMye) have been removed to avoid confusion, as these outputs are not relevant to structural loading and could be misleading.

Related issue, if one exists
#2325, #855

Impacted areas of the software
SubDyn

Test results, if applicable
The r-test results have been updated to reflect expected outputs. Test 4 is a newly added r-test to verify self-weight outputs for a system that reorients in space (floating system).

Test 1: Vertical cantilever beam (from GitHub 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 2a: Horizontal cantilever beam (from GitHub 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 2b: Same horizontal cantilever beam as Test 2a, but under water
This test was discused in the comments below and has been summarized here.

The system considered is the same as the one shown in Test 2a, but the beam longitudinal axis is placed 10 m below the still water level and the beam is considered massless.

By integrating the hydrostatic pressure, we get a vertical buoyancy force and one axial compression force at the end plate.

The distributed vertical load can be computed as follows:
q = ρ·g·A = 1025·9.80665·(π·((5/2)^2)) = 197.367 kN/m

The hydrostatic pressure on the vertical circular endplate acts primarily as an axial compressive force, but its line of action does not pass through the plate centroid because the hydrostatic pressure increases linearly with depth. As a result, this offset creates a lever arm, which in turn generates a bending moment.

The resultant hydrostatic force acting normal to the endplate is: F = ρ·g·h·A = 1025·9.80665·10·pi·(5/2)^2= 1973.67 kN

For a vertical circular plate, the center of pressure offset from the centroid is: e = (R^2)/(4h) = (2.5^2)/(4·10) = 0.1562 m

The resulting bending moment is: 1973.67·0.1562 = 308.4 kNm

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 consistent.

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#183)

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.

[luwang00]

After further discussion with @RBergua, we have decided to put this PR on hold for now.

While this PR provides a more accurate treatment of structural self-weight when computing the sectional load outputs from SubDyn, it does not address the fact that the hydrostatic and hydrodynamic loads from HydroDyn are also lumped and mapped to the SubDyn structural nodes. These load components are often of similar magnitudes as the structural self-weight. The modification proposed here can lead to further inconsistencies, i.e., distributed structural self-weight and lumped hydrostatic and hydrodynamic loads, rendering the interpretation of SubDyn outputs more challenging.

Unfortunately, it is difficult to formulate and implement similar treatments for the hydrodynamic and hydrostatic loads as for the self-weight. Until a more permanent solution can be developed, we shall postpone the changes proposed here. Lastly, note that as the mesh discretization is increased in both SubDyn and HydroDyn with the HydroDyn mesh being either aligned to the SubDyn mesh or considerably finer, the load output discrepancies observed here will shrink, eventually converging to the expected values.

[RBergua]

I'm leaving this here for future reference. I took the system from the Test 2 shown above and placed it 10 m below the mean sea level. In this case, I removed the beam mass and only included the hydrostatic force (i.e., vertical buoyancy and axial compression). The distributed vertical load can be computed as follows:
q = ρ·g·A = 1025·9.80665·(π·((5/2)^2)) = 197.367 kN/m

The results from the analytical solution and the SubDyn output are:

For a distributed load (q), the consistent load vector for this horizontal beam would be:

In this case, I can reproduce the shear force (246.7 kN). However, I don't seem to reproduce the bending moment. Analytically I have 102.8 kNm while OpenFAST seems to be using 308.4 kNm (there seems to be a factor of 3: 308.4/102.8). For the self-weight calculation shown above, the consistent load vector is aligned with the expected one. But in this case (hydrostatic load), the bending moments seem to be different.

[RBergua]

After discussing this with @luwang00, it appears that the 308.4 kNm bending moment indeed arises from the hydrostatic center of pressure acting on the vertical circular endplate. Although this hydrostatic load acts primarily as an axial compressive force, its line of action does not pass through the plate centroid.

Hydrostatic pressure increases linearly with depth. As a result, the center of pressure on a vertical circular endplate lies below the centroid (where the structural node is located). This offset creates a lever arm, which in turn generates a bending moment.

The resultant hydrostatic force acting normal to the endplate is: F = ρ·g·h·A = 1025·9.80665·10·pi·(5/2)^2= 1973.67 kN

For a vertical circular plate, the center of pressure offset from the centroid is: e = (R^2)/(4h) = (2.5^2)/(4·10) = 0.1562 m

Therefore, the resulting bending moment is: 1973.67·0.1562 = 308.4 kNm

This confirms that the reported 308.4 kNm moment is the moment generated by the hydrostatic resultant force acting through the center of pressure.

This also indicates that the distributed load in the vertical direction due to the hydrostatic pressure integration, do not correspond to a consistent element load vector. They include the force component, but not the bending moment.

[RBergua]

@luwang00 In commit 6857d52 Fg is now handled separately for beam and cable elements in member force recovery:

Beam elements: self-weight force DOFs are zeroed in Fg; only the gravity fixed-end bending moments are kept and recomputed from the current element orientation for floating systems.

Cable elements: Fg is replaced by the initial pretension FCe.

I understand that your force extrapolation that handles hydro loads should be able to handle the remaining beam self-weight force contributions on top of this. Let me know if we are not aligned!

[luwang00]

Additional changes from PR #3367 have been merged in to have a combined pull request. See PR #3367 for more information.

[luwang00]

As part of this PR, Issue #855 should also be addressed.

[luwang00]

• Merge r-test requested. See PR New test for SubDyn self-weight r-test#183