SubDyn: Consistent postprocessing for self-weight reconstruction

[RBergua]

This is still a work in progress. The current development addresses the consistent reconstruction of SubDyn self-weight effects. However, for offshore systems, a similar reconstruction procedure should also be implemeneted for hydro loads to ensure the reported outputs remain consistent.

Feature or improvement description
SubDyn computes beam self-weight using 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 associated with 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 should be applied to any load that enters the system as a distributed load and 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 any correction.

Interestingly, the source code suggests 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

The above implementation behaves as expected for fixed-bottom systems, where the self-weight force and moment contributions remain constant. For floating systems, however, the self-weight components computed during initialization had to be projected onto the current floating-body reference frame before being reported as outputs. In addition, the self-weight bending-moment components had to be recomputed from scratch at every time step using the actual element direction-cosine matrix.

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

I think it would be also good to cleanly separate the cable pretension contribution from the gravity contribution within ElemF.

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.

[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. In the short term, we can potentially recommend aligning SubDyn and HydroDyn mesh nodes and change the SubDyn section load outputs to be at the midpoint of an element instead of at the end nodes. This should provide outputs that are easier to interpret physically.

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.