Sum of all forces on a FOWT substructure

[Victor.Rappe]

Hello all,

I am currently working on a model to map all relevant loads from an openFAST simulation, which considers a rigid substructure, to a shell-based finite element model.

To verify this mapping strategy, I would then consider the sum of all loads and the sum of all moments, which in theory should be equal to zero. However, this is not the case in both my finite element simulations and my openFAST simulations.

The openFAST model has following details:

• Version: custom version of OpenFAST-e0bc0cf9-dirty (number of member level outputs was extended to 9 points on 16 members)
• Aerodyn is disabled
• Still water conditions or periodic waves with H_s=10m and T_p=9s
• Potential flow model: PotMod = 1, ExctnMod = 1, RdtnMod = 0
• No second order loads are considered
• No Morison drag loads are considered

My method for calculating the sum of all forces in a local reference from linked to the floater goes as follows:
1.Calculate the rotation matrix R(roll,pitch,yaw) to transform load vectors from the global to the local reference system

2.Hydrodynamic loads: HydroFxi, HydroFyi, HydroFzi

• From the documentation:

When coupled to FAST, the hydrodynamic loads computed by HydroDyn are applied to the displaced position of the substructure (i.e., the displaced platform in ElastoDyn and/or the displaced substructure in SubDyn), but are based on wave kinematics at the undisplaced position.

• These forces are directly applied in the local coordinate system (CSYS)

3.Tower base: TwrBsFxt, TwrBsFyt, TwrBsFzt

• From the documentation: directed along the xt axis (=local coordinate system?)
• These forces are directly applied in the local coordinate system

4.Mooring line loads: TwrBsFxt, TwrBsFyt, TwrBsFzt

• Calculate the orientation of the mooring line
• Calculate the force components of the mooring line loads via:
  F_(mooring,i, global)=FAIRTENi direction/(|direction|)
• Transform from the global CSYS to the local CSYS:
  F_(mooring,i) =R F_(mooring,i, global)

5.Gravity load

• Transform the gravity load vector to the local CSYS:
  F_gravity=R (m_structure [0,0,-9.81])

6.Inertia loads

• Extract: PtfmTAxi, PtfmTAyi, PtfmTAzi
• Define: F_(inertia,global)=m_structure [PtfmTaxi,PtfmTAyi,PtfmTAzi]
• Transform to local CSYS: F_inertia=R F_(inertia,global)

7.Take the sum of all loads

• Σ_i F_i = HydroF+TwrBsF+∑_j F_(mooring,j) +F_gravity+F_inertia =?= 0

If this is implemented, the following results are found for no wind and no wave conditions:

The forces in y and z appeared to die out relatively quickly, but there was something wrong with F_x. This was solved by not transforming the gravity loads, i.e. the gravity loads are considered to be m_structure [0,0,-9.81] in the local CSYS:

These results are relatively good, however note that F_y and F_z are centered around 0, but that F_x is not.

When waves are introduced (and the potential flow solution becomes important), the results become extremely bad:

The loads are still centered around 0, however the amplitude in F_x and F_z is in the order of 10MN.

If someone has any ideas of what is going wrong and is willing to share them, I’d be very grateful.

Kind regards,
Victor Rappe

[Jason.Jonkman]

Dear @Victor.Rappe,

Just a few comments:

• The total integrated hydrodynamic loads (HydroFxi, etc.) are expressed in global coordinates, so, you’d half to transform these to local coordinates (by pre-multiplying by R) to compare to the other local coordinate system forces.
• I’m not sure I fully understand how you are calculating the mooring loads in local coordinates.
• The gravity load should be transformed to local coordinates for consistency with the other local coordinate system forces.
• The equation F = m * a assumes that the forces and accelerations are computed at the center of mass of the platform. However, the platform acceleration outputs you are using are not at the platform center of mass unless you’ve set PtfmCMzt = PtfmRefzt. Also, your inertia force should have a negative sign, bringing m * a to the left-hand side.

Best regards,

[Victor.Rappe]

Dear @Jason.Jonkman,

Thank you very much for your quick reply.

For calculating the mooring line loads I do the following:

1. Extract the coordinates of the node of the mooring line which is connected to the floater and the coordinates of the next point from the ___.MD.Linei.out files.
   For my simulation, there are 20 elements per mooring line, meaning that I extract the coordinates of node 19 and node 20.
2. Calculate the normalised vector that is tangential to the mooring line at the platform:
   vector = (XYZ_{node19} - XYZ{node20}) / |XYZ_{node19} - XYZ{node20}|
3. Calculate the mooring line loads in the global CSYS via:
   F_{mooring, glob, i} = FAIRTENi * vector
4. Transform to the local CSYS via:
   F_{mooring, loc, i} = R F_{mooring, glob, i}

After implementing your comments I now get following results for a simulation with both (linear) PF and the drag term of the Morison equation implemented:

While the magnitude of F_y seems to be in the realm of numerical noise, F_x and F_z have oscillations with amplitude of around 400kN and 25kN, respectively. When you compare this to the buoyancy force as a characteristic force of the floating body (≈140MN), these loads are small, but I wonder what your opinion is on these results.

Also note that the oscillations of F_y and F_z are not centered around 0. Is this (un)expected behaviour?

Kind regards,
Victor

EDIT:
The updated method now goes like this:

1.Calculate the rotation matrix R(roll,pitch,yaw) to transform load vectors from the global to the local reference system

2.Hydrodynamic loads: HydroFxi, HydroFyi, HydroFzi

• Transfrom to the local CSYS: F_hydro = R [HydroFxi, HydroFyi, HydroFzi]

3.Tower base: TwrBsFxt, TwrBsFyt, TwrBsFzt

• From the documentation: directed along the xt axis (=local coordinate system?)
• These forces are directly applied in the local coordinate system

4.Mooring line loads: TwrBsFxt, TwrBsFyt, TwrBsFzt

• Calculate the global force components of the mooring line at the floating platform F_(mooring,i, global) (see above)
• Transform from the global CSYS to the local CSYS:
  F_(mooring,i) =R F_(mooring,i, global)

5.Gravity load

• Transform the gravity load vector to the local CSYS:
  F_gravity=R (m_structure [0,0,-9.81])

6.Inertia loads

• Extract: a_ref = [PtfmTAxi, PtfmTAyi, PtfmTAzi]
• Extract: omega = [PtfmRVxi,PtfmRVyi,PtfmRVzi]
• Extract: alpha = [PtfmRAxi,PtfmRAyi,PtfmRAzi]
• Define: a_cog = a_ref + alpha x r + omega x (omega x r), with r the vector between the reference point and the COG
• Define: F_(inertia,global)=m_structure a_cog
• Transform to local CSYS: F_inertia=R F_(inertia,global)

7.Take the sum of all loads

• Σ_i F_i = HydroF+TwrBsF+∑_j F_(mooring,j) +F_gravity+F_inertia =?= 0

[Jason.Jonkman]

Dear @Victor.Rappe,

Thanks for clarifying regarding the mooring line forces.

Regarding your change to the inertia forces, I agree with your changes, assuming you are using consistent units for the omega (rad/s) and alpha (rad/s^2); that said, I still don’t see that you flipped the sign on F_Inertia.

Can you clarify how you are calculating the rotation matrix R?

Yes, I would expect the sum of forces to equal zero.

Best regards,

[Victor.Rappe]

Dear @Jason.Jonkman,

Omega and alpha are indeed expressed in rad/s and rad/s², respectively.

The sign of F_inertia was flipped in my script, but I forgot to change it in my previous comment.

The rotation matrix R is calculated as follows:

def coord_transform(vector, roll, pitch, yaw):
        """Transform vectors from the global inertial CSYS to the floater CSYS.
        
        Parameters
        ----------
        vector : np.array
            The vector in the global CSYS (t x 3).

        Returns
        -------
        new_vector : np.array
            The vector in the floater CSYS (t x 3).
        
        """     
        
        # method 1: using a for loop
        # new_vector = np.empty_like(vector)
        # for i in range(len(roll)):
        #     R_x = np.array([[1, 0, 0],
        #                     [0, np.cos(roll[i]), -np.sin(roll[i])],
        #                     [0, np.sin(roll[i]), np.cos(roll[i])]])
        #     R_y = np.array([[np.cos(pitch[i]), 0, np.sin(pitch[i])],
        #                     [0, 1, 0],
        #                     [-np.sin(pitch[i]), 0, np.cos(pitch[i])]])
        #     R_z = np.array([[np.cos(yaw[i]), -np.sin(yaw[i]), 0],
        #                     [np.sin(yaw[i]), np.cos(yaw[i]), 0],
        #                     [0, 0, 1]])
        #     R = R_x @ R_y @ R_z

        #     new_vector[i,:] = R @ vector[i,:]

        # method 2: using numpy broadcasting
        R = np.empty((len(roll), 3, 3))

        R[:, 0, 0] = np.cos(yaw) * np.cos(pitch)
        R[:, 0, 1] = np.cos(yaw) * np.sin(pitch) * np.sin(roll) - np.sin(yaw) * np.cos(roll)
        R[:, 0, 2] = np.cos(yaw) * np.sin(pitch) * np.cos(roll) + np.sin(yaw) * np.sin(roll)

        R[:, 1, 0] = np.sin(yaw) * np.cos(pitch)
        R[:, 1, 1] = np.sin(yaw) * np.sin(pitch) * np.sin(roll) + np.cos(yaw) * np.cos(roll)
        R[:, 1, 2] = np.sin(yaw) * np.sin(pitch) * np.cos(roll) - np.cos(yaw) * np.sin(roll)

        R[:, 2, 0] = -np.sin(pitch)
        R[:, 2, 1] = np.cos(pitch) * np.sin(roll)
        R[:, 2, 2] = np.cos(pitch) * np.cos(roll)

        new_vector = np.einsum('ijk,ik->ij', R, vector)

        return new_vector

It has been verified that the method using numpy broadcasting gives the same results as the for loop.
Roll, pitch and yaw are defined as:

roll = df_HD["PtfmRDxi"].values * np.pi / 180 # rad
pitch = df_HD["PtfmRDyi"].values * np.pi / 180 # rad
yaw = df_HD["PtfmRDzi"].values * np.pi / 180 # rad

Kind regards,
Victor

[Jason.Jonkman]

Dear @Victor.Rappe,

Perhaps the difference is small if the platform rotations are small, but version of OpenFAST you are using makes use of a different rotational convention for the platform rotation DOFs in ElastoDyn and HydroDyn. The rotation matrix used instead is documented in Equation (2-2) of my Ph.D. thesis-turned NREL report: https://www.nrel.gov/docs/fy08osti/41958.pdf.

FYI: This rotation convention is changing in the upcoming release of OpenFAST v4, which is being modified to support large rotation of the floater. I’ll present this work at the upcoming DeepWind 2025 conference.

Best regards,

[Victor.Rappe]

Dear @Jason.Jonkman,

Sorry in advance for the lenghty post, but I still have not arrived at a satisfactory result.

Using your small-angle approximation rotation matrix has made a big difference. However the sum of all forces is still not zero.

Using a newly compiled openFAST v3.5.4, I did a few simulations to try and isolate the problem:

a) Full turbine + platform

1. No wind, no waves, with potential flow (PF) and Morison drag term
   

   

2. No wind, no waves, only Morison drag term
   

   

3. No wind, no waves, only PF
   

   

Although there seem to be small differences between the three cases, the thing that seems to be most important to me is that ΣFz is the only component that is centered around 600N instead of 0N.

b) Only the floating substructure + mooring lines
I followed this post to remove the tower and turbine blades from the simulation. Note that this meant that there was a significant movement in heave, as I did not alter the parameters of the floating structure. This gave me following results:

1. No wind, no waves, PF + Morison drag
   

   

2. No wind, waves, PF + Morison drag
   

   

While removing the tower significantly reduced ΣFx and ΣFy (both are below 1N), an offset can still be seen in ΣFz. When waves were included, loads similar to the simulations with tower were again found.

I currently think that the issue lies in the inertial loads (see the end of the post for the code I used) for the following reasons:

• It cannot be the tower base loads as we still see the same errors when considering waves on the structure without the tower
• I don’t think it is the hydrodynamic loads, as changing between PF and Morison, PF only and Morison does not produce significant differences.
• The overhang of the turbine causes the substructure to pitch, leading to a sinusoidal load in x and z in the local CSYS.
  
  
  
  
  
  • This could explain why ΣFx is extremely small in case b1 (no imbalance in the structure, as there is no tower), but becomes significant again in case b2 (waves causes pitching).
  • For the same case of b1, ΣFz is the only component which is largely different from zero. This could be linked with the heave motion of the structure due to removing the tower:
    
    
    
    Note the similarity between this acceleration and ΣFz in case b1.
    However, this does not explain the offset which can be seen in ΣFz .

The thing I cannot explain is the offset which is always visible in ΣFz.

Given all this I have a few questions:

• Do you have any idea what could still be wrong or see any links that I have not yet seen?
• What order of magnitude of the forces (and moments) do you think is acceptable if we look at the sum of all loads?
  • Is this an absolute number or is this dependent on the scale and mass of the structure that is being considered?
• Do you know of any resources (simulation data, verification data, literature,…) that also considers the sum of all forces? The closest I can find are papers that compare the dynamics calculated in one program to other programs.

I look forward to your response,
Victor

---

Addendum
Below you can find the code I used to calculate the inertial loads. If necessary I can share the full code I use.

# Load in data
df_HD = ... # .out file
df_ML1 = ... # MD.line.out file
df_ML2 = ....
df_ML3 = ...

# Transformation matrix
def coord_transform_Jonkman(vector, theta1, theta2, theta3):
    """Transform vectors from the global inertial CSYS to the floater CSYS, based on PhD of Jason Jonkman. https://www.nrel.gov/docs/fy08osti/41958.pdf 
    
    Parameters
    ----------
    vector : np.array
        The vector in the global CSYS (t x 3).

    Returns
    -------
    new_vector : np.array
        The vector in the floater CSYS (t x 3).
    
    """     
    # Small angle approximation
    R = np.empty((len(roll), 3, 3))

    noemer = (theta1**2 + theta2**2 + theta3**2) * np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) + 1e-40 # to avoid error messages due to division by zero

    R[:, 0, 0] = (theta1**2 * np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) + theta2**2 + theta3**2) / noemer
    R[:, 0, 1] = (theta3 * (theta1**2 + theta2**2 + theta3**2) + theta1 * theta2 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer
    R[:, 0, 2] = (-theta2 * (theta1**2 + theta2**2 + theta3**2) + theta1 * theta3 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer

    R[:, 1, 0] = (-theta3 * (theta1**2 + theta2**2 + theta3**2) + theta1 * theta2 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer
    R[:, 1, 1] = (theta1**2 + theta2**2 * np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) + theta3**2) / noemer
    R[:, 1, 2] = (theta1 * (theta1**2 + theta2**2 + theta3**2) + theta2 * theta3 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer

    R[:, 2, 0] = (theta2 * (theta1**2 + theta2**2 + theta3**2) + theta1 * theta3 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer
    R[:, 2, 1] = (-theta1 * (theta1**2 + theta2**2 + theta3**2) + theta2 * theta3 * (np.sqrt(1 + theta1**2 + theta2**2 + theta3**2) - 1)) / noemer
    R[:, 2, 2] = (theta1**2 + theta2**2 + theta3**2 * np.sqrt(1 + theta1**2 + theta2**2 + theta3**2)) / noemer
    
    new_vector = np.einsum('ijk,ik->ij', R, vector)


roll = df_HD["PtfmRDxi"].values * np.pi / 180 # rad
pitch = df_HD["PtfmRDyi"].values * np.pi / 180 # rad
yaw = df_HD["PtfmRDzi"].values * np.pi / 180 # rad

# Inertial loads
time = df_HD["Time"].values

## Translational acceleration
a_x = df_HD["PtfmTAxi"].values
a_y = df_HD["PtfmTAyi"].values
a_z = df_HD["PtfmTAzi"].values

a_ref = np.array([a_x, a_y, a_z]).T

## Rotational velocity
omega_x = df_HD["PtfmRVxi"].values * np.pi / 180 # rad/s
omega_y = df_HD["PtfmRVyi"].values * np.pi / 180 # rad/s
omega_z = df_HD["PtfmRVzi"].values * np.pi / 180 # rad/s

omega = np.array([omega_x, omega_y, omega_z]).T

## Rotational acceleration
alpha_x = df_HD["PtfmRAxi"].values * np.pi / 180 # rad/s^2
alpha_y = df_HD["PtfmRAyi"].values * np.pi / 180 # rad/s^2
alpha_z = df_HD["PtfmRAzi"].values * np.pi / 180 # rad/s^2

alpha = np.array([alpha_x, alpha_y, alpha_z]).T

## Translational inertia
r = [0,0,-12.0527] # m --> distance from reference point (=(0,0,0)) to center of gravity NOTE: for own updated version of OC4 DeepCWind semisubmersible

a_cog = a_ref + np.cross(alpha, r) + np.cross(omega, np.cross(omega, r))

F_inertia_glob = M * a_cog

F_inertia = coord_transform_Jonkman(F_inertia_glob, roll, pitch, yaw)

[Jason.Jonkman]

Dear @Victor.Rappe,

I’ve certainly seen similar work performed before, but I can’t think of a specific reference to cite you.

I would suggest normalizing the force sum relative to a nominal force, e.g., normalize the force sum in z to the mean load of HydroFzi, which I is likely quite large. Perhaps this normalized error is quite small in all cases, which could imply that it is simply numerical round-off that you are tracking. You could always increase the precision of the OpenFAST outputs by using binary format or an ASCII format with high precision (by changing OutFmt).

Best regards,

[Victor.Rappe]

Dear @Jason.Jonkman,

I changed the output from ES15.7E2 to ES20.11E3 as described in the documentation, but this did unfortunately not have any meaningful impact on my results.

The mean load of HydroFzi (and of the resultant of the total hydrodynamic load) is approximately equal to the buoyancy force, so I took this as the normalisation factor.

For multiple simulations I find values of the mean total resultant force to be consistently below 0.05% of the hydrostatic load.
An example for 8m/s wind, waves with Hs=10m and Tp=9s with heading 0° is given below.

I now tried to do the same for the moments and once again the results are quite bad (see below for some examples).
The method I followed is quite similar to the one for the sum of all forces:

1. Calculate the rotation matrix (small angle approximation)
2. Hydrodynamic loads: M_hydro_glob = [HydroMxi, HydroMyi, HydroMzi]
   • Transform to local CSYS: M_hydro = R M_hydro_glob
   • Calculate the extra moments around the CoG due to F_hydro_glob = [HydroFxi, HydroFyi, HydroFzi] working at (0,0,0):
     • Forces to local CSYS: F_hydro = R F_hydro_glob
     • Add extra moments: M_hydro += r_cog x F_hydro_glob
3. Tower base loads: M_twr = [TwrBsMxt, TwrBsMyt, TwrBsMzt]
   • Calculate the extra moments around the CoG due to F_twr = [TwrBsFxt, TwrBsFyt, TwrBsFzt] working at (0,0,0):
     M_twr += r_cog x F_twr
4. Mooring line loads: TwrBsFxt, TwrBsFyt, TwrBsFzt
   • Calculate F_(mooring,i) = R F_(mooring,i, global) in the local CSYS as described before.
   • Calculate the vector r_i between the CoG and the attachment point of the mooring line i.
   • Calculate the moments due to the mooring line forces:
     M_(mooring,i) = r_i x F_(mooring,i)
5. Gravity loads: no moments around CoG
6. Inertia loads:
   • Define: I = [[Ix,0,0],[0,Iy,0],[0,0,Iz]]
   • Define: omega_local = R [PtfmRVxi,PtfmRVyi,PtfmRVzi]
   • Define: alpha_local = R [PtfmRAxi,PtfmRAyi,PtfmRAzi]
   • Calculate: M_inertia = I · alpha_local + omega_local x (I · omega_local)
7. Take the sum of all loads:
   • Σ_i M_i = M_hydro + M_twr+∑_j M_(mooring,j) - F_inertia =?= 0

When following this method I get following results:

1. No tower, no wind, no waves
   

   
   
   Notice the similarity between ΣMy, the pitch and the hydrostatic restoring moment of the platform:
   
   
   
   
   

2. No tower, no wind, waves, heading = 0°
   

   
   
   Notice the similarity between ΣMy and the hydrostatic restoring moment of the platform:
   
   

3. Tower, no wind, no waves
   

   
   
   Notice the similarity between ΣMy, the pitch and the hydrostatic restoring moment of the platform:
   
   
   
   
   

4. Tower, wind, waves, heading = 0°
   

   
   
   Notice the similarity between ΣMx and the hydrostatic restoring moment:
   
   

5. Tower, wind, waves, heading = 45°
   

   

I suspect that the issue lies somewhere with the roll, pitch and yaw of the platform, but I cannot find what exactly. I personally suspect it has something to do with the hydrostatics of the structure, but as this load is included in HydroMxyzi this should not be the problem.

Lastly, I looked for a moment similar to the buoyancy load to normalise the moments. As the average of the total hydrodynamic restoring moment is close to zero, I was thinking of using M_restoring = rho * g * V0 * GM * sin(theta), with theta either equal to the maximum of the roll or pitch of the structure or a constant of 10° (as a maximum angle that should be achieved for a FOWT).

Kind regards,
Victor

[Jason.Jonkman]

Dear @Victor.Rappe,

I’m glad to hear that the normalized force sums are all now close to zero.

I generally agree with your approach for the moments, but I don’t have anything specific to add.

Best regards,