Switching from 'RtAero' to TwrBs in an offshore test

I noticed that in the FAST OutData, for an offshore test, ‘TwrBsFxt’, ‘TwrBsFyt’, ‘TwrBsFxt’, ‘TwrBsMxt’, ‘TwrBsMyt’, ‘TwrBsMzt’ also contain inertial loads, and that moreover they are integral with the tower and not with the inertial reference system. My goal is to have the loads at the base of the tower (actually in the origin of the inertial frame of reference, (0,0,0)) due exclusively to the aerodynamic loads on the blades, excluding the inertial loads.
Since there are no such outputs in the OutData, I wrote the formulations to get them starting from ‘RtAeroFxh’, ‘RtAeroFyh’, ‘RtAeroFzh’, ‘RtAeroMxh’, ‘RtAeroMyh’, ‘RtAeroMzh’.
Please note that my simulations consider rigid bodies, with the following Elastodyn setting:

------- ELASTODYN v1.03.* INPUT FILE ------------------------------------------- NREL 5.0 MW Baseline Wind Turbine for Use in Offshore Analysis. Properties from Dutch Offshore Wind Energy Converter (DOWEC) 6MW Pre-Design (10046_009.pdf) and REpower 5M 5MW (5m_uk.pdf) ---------------------- SIMULATION CONTROL -------------------------------------- True Echo - Echo input data to "<RootName>.ech" (flag) 3 Method - Integration method: {1: RK4, 2: AB4, or 3: ABM4} (-) "default" DT - Integration time step (s) ---------------------- ENVIRONMENTAL CONDITION --------------------------------- 9.80665 Gravity - Gravitational acceleration (m/s^2) ---------------------- DEGREES OF FREEDOM -------------------------------------- False FlapDOF1 - First flapwise blade mode DOF (flag) False FlapDOF2 - Second flapwise blade mode DOF (flag) False EdgeDOF - First edgewise blade mode DOF (flag) False TeetDOF - Rotor-teeter DOF (flag) [unused for 3 blades] False DrTrDOF - Drivetrain rotational-flexibility DOF (flag) True GenDOF - Generator DOF (flag) False YawDOF - Yaw DOF (flag) False TwFADOF1 - First fore-aft tower bending-mode DOF (flag) False TwFADOF2 - Second fore-aft tower bending-mode DOF (flag) False TwSSDOF1 - First side-to-side tower bending-mode DOF (flag) False TwSSDOF2 - Second side-to-side tower bending-mode DOF (flag) True PtfmSgDOF - Platform horizontal surge translation DOF (flag) True PtfmSwDOF - Platform horizontal sway translation DOF (flag) True PtfmHvDOF - Platform vertical heave translation DOF (flag) True PtfmRDOF - Platform roll tilt rotation DOF (flag) True PtfmPDOF - Platform pitch tilt rotation DOF (flag) True PtfmYDOF - Platform yaw rotation DOF (flag) ---------------------- INITIAL CONDITIONS -------------------------------------- 0 OoPDefl - Initial out-of-plane blade-tip displacement (meters) 0 IPDefl - Initial in-plane blade-tip deflection (meters) 0 BlPitch(1) - Blade 1 initial pitch (degrees) 0 BlPitch(2) - Blade 2 initial pitch (degrees) 0 BlPitch(3) - Blade 3 initial pitch (degrees) [unused for 2 blades] 0 TeetDefl - Initial or fixed teeter angle (degrees) [unused for 3 blades] 0 Azimuth - Initial azimuth angle for blade 1 (degrees) 8.1334 RotSpeed - Initial or fixed rotor speed (rpm) 0 NacYaw - Initial or fixed nacelle-yaw angle (degrees) 0 TTDspFA - Initial fore-aft tower-top displacement (meters) 0 TTDspSS - Initial side-to-side tower-top displacement (meters) 0 PtfmSurge - Initial or fixed horizontal surge translational displacement of platform (meters) 0 PtfmSway - Initial or fixed horizontal sway translational displacement of platform (meters) 0 PtfmHeave - Initial or fixed vertical heave translational displacement of platform (meters) 0 PtfmRoll - Initial or fixed roll tilt rotational displacement of platform (degrees) 0 PtfmPitch - Initial or fixed pitch tilt rotational displacement of platform (degrees) 0 PtfmYaw - Initial or fixed yaw rotational displacement of platform (degrees) ---------------------- TURBINE CONFIGURATION ----------------------------------- 3 NumBl - Number of blades (-) 63 TipRad - The distance from the rotor apex to the blade tip (meters) 1.5 HubRad - The distance from the rotor apex to the blade root (meters) 0 PreCone(1) - Blade 1 cone angle (degrees) 0 PreCone(2) - Blade 2 cone angle (degrees) 0 PreCone(3) - Blade 3 cone angle (degrees) [unused for 2 blades] 0 HubCM - Distance from rotor apex to hub mass [positive downwind] (meters) 0 UndSling - Undersling length [distance from teeter pin to the rotor apex] (meters) [unused for 3 blades] 0 Delta3 - Delta-3 angle for teetering rotors (degrees) [unused for 3 blades] 0 AzimB1Up - Azimuth value to use for I/O when blade 1 points up (degrees) 0 OverHang - Distance from yaw axis to rotor apex [3 blades] or teeter pin [2 blades] (meters) 0 ShftGagL - Distance from rotor apex [3 blades] or teeter pin [2 blades] to shaft strain gages [positive for upwind rotors] (meters) 0 ShftTilt - Rotor shaft tilt angle (degrees) 0 NacCMxn - Downwind distance from the tower-top to the nacelle CM (meters) 0 NacCMyn - Lateral distance from the tower-top to the nacelle CM (meters) 1.75 NacCMzn - Vertical distance from the tower-top to the nacelle CM (meters) 0 NcIMUxn - Downwind distance from the tower-top to the nacelle IMU (meters) 0 NcIMUyn - Lateral distance from the tower-top to the nacelle IMU (meters) 2.23336 NcIMUzn - Vertical distance from the tower-top to the nacelle IMU (meters) 1.96256 Twr2Shft - Vertical distance from the tower-top to the rotor shaft (meters) 87.6 TowerHt - Height of tower above ground level [onshore] or MSL [offshore] (meters) 10 TowerBsHt - Height of tower base above ground level [onshore] or MSL [offshore] (meters) 0 PtfmCMxt - Downwind distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters) 0 PtfmCMyt - Lateral distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters) -89.9155 PtfmCMzt - Vertical distance from the ground level [onshore] or MSL [offshore] to the platform CM (meters) -0 PtfmRefzt - Vertical distance from the ground level [onshore] or MSL [offshore] to the platform reference point (meters) ---------------------- MASS AND INERTIA ---------------------------------------- 0 TipMass(1) - Tip-brake mass, blade 1 (kg) 0 TipMass(2) - Tip-brake mass, blade 2 (kg) 0 TipMass(3) - Tip-brake mass, blade 3 (kg) [unused for 2 blades] 56780 HubMass - Hub mass (kg) 115926 HubIner - Hub inertia about rotor axis [3 blades] or teeter axis [2 blades] (kg m^2) 534.116 GenIner - Generator inertia about HSS (kg m^2) 240000 NacMass - Nacelle mass (kg) 2.60789E+06 NacYIner - Nacelle inertia about yaw axis (kg m^2) 0 YawBrMass - Yaw bearing mass (kg) 7.46633E+06 PtfmMass - Platform mass (kg) 4.22923E+09 PtfmRIner - Platform inertia for roll tilt rotation about the platform CM (kg m^2) 4.22923E+09 PtfmPIner - Platform inertia for pitch tilt rotation about the platform CM (kg m^2) 1.6423E+08 PtfmYIner - Platform inertia for yaw rotation about the platform CM (kg m^2) ---------------------- BLADE --------------------------------------------------- 17 BldNodes - Number of blade nodes (per blade) used for analysis (-) "NRELOffshrBsline5MW_Blade.dat" BldFile(1) - Name of file containing properties for blade 1 (quoted string) "NRELOffshrBsline5MW_Blade.dat" BldFile(2) - Name of file containing properties for blade 2 (quoted string) "NRELOffshrBsline5MW_Blade.dat" BldFile(3) - Name of file containing properties for blade 3 (quoted string) [unused for 2 blades] ---------------------- ROTOR-TEETER -------------------------------------------- 0 TeetMod - Rotor-teeter spring/damper model {0: none, 1: standard, 2: user-defined from routine UserTeet} (switch) [unused for 3 blades] 0 TeetDmpP - Rotor-teeter damper position (degrees) [used only for 2 blades and when TeetMod=1] 0 TeetDmp - Rotor-teeter damping constant (N-m/(rad/s)) [used only for 2 blades and when TeetMod=1] 0 TeetCDmp - Rotor-teeter rate-independent Coulomb-damping moment (N-m) [used only for 2 blades and when TeetMod=1] 0 TeetSStP - Rotor-teeter soft-stop position (degrees) [used only for 2 blades and when TeetMod=1] 0 TeetHStP - Rotor-teeter hard-stop position (degrees) [used only for 2 blades and when TeetMod=1] 0 TeetSSSp - Rotor-teeter soft-stop linear-spring constant (N-m/rad) [used only for 2 blades and when TeetMod=1] 0 TeetHSSp - Rotor-teeter hard-stop linear-spring constant (N-m/rad) [used only for 2 blades and when TeetMod=1] ---------------------- DRIVETRAIN ---------------------------------------------- 100 GBoxEff - Gearbox efficiency (%) 97 GBRatio - Gearbox ratio (-) 8.67637E+08 DTTorSpr - Drivetrain torsional spring (N-m/rad) 6.215E+06 DTTorDmp - Drivetrain torsional damper (N-m/(rad/s)) ---------------------- FURLING ------------------------------------------------- False Furling - Read in additional model properties for furling turbine (flag) [must currently be FALSE) "unused" FurlFile - Name of file containing furling properties (quoted string) [unused when Furling=False] ---------------------- TOWER --------------------------------------------------- 20 TwrNodes - Number of tower nodes used for analysis (-) "NRELOffshrBsline5MW_OC3Hywind_ElastoDyn_Tower.dat" TwrFile - Name of file containing tower properties (quoted string)

The loads at the nacelle (i.e. in a reference system integral with the tower with z_nac = z_tow and x_nac = x_rot) are:
Fx_nac=RtAeroFxh;
Fy_nac=RtAeroFyh.cos(Azimuthpi/180)-RtAeroFzh.sin(Azimuthpi/180);
Fz_nac=RtAeroFyh.sin(Azimuthpi/180)+RtAeroFzh.cos(Azimuthpi/180);
Mx_nac=RotAeroMxh;
My_nac=RtAeroMyh.cos(Azimuthpi/180)-RtAeroMzh.sin(Azimuthpi/180);
Mz_nac=RtAeroMzh.cos(Azimuthpi/180)+RtAeroMyh.sin(Azimuthpi/180);

Hence, to obtain the loads in (0,0,0), defined TransMat the Euler matrix containing the angular positions (PtfmRoll, PtfmPitch, PtfmYaw) in the inertial reference frame:
F = TransMat * [Fx_nac,Fy_nac,Fz_nac]
M = TransMat * [Mx_nac,My_nac,Mz_nac] + (TransMat * [0, 0, Hrot]) x F
where Hrot is the vertical distance between inertial reference frame (0,0,0) and nacelle reference frame (0,0,Hrot).
Gravity force and gravity restoring moments are not considered in these formulas as they are taken into account using the center of mass coordinates of the entire system.
In this way, the vector F of the 3 forces and the vector M of the 3 moments in the inertial reference system (0,0,0) due to the turbine loads are obtained.
I wonder if you agree with my procedure.
Best regards,
Lorenzo.

Dear Lorenzo,

I agree with your calculation considering that you’ve set ShftTilt = OverHang = 0. And, presumably, Hrot = TowerHt+Twr2Shft. And I presume you want (0,0,0) to translate with the platform transnational DOFs rather than be fixed to the inertial frame?

Then, the only question is how you’ve defined TransMat. You mention an Euler matrix, but ElastoDyn would not use this for the platform rotations, which are assumed to be only moderate in size. Instead, TransMat should equal the transpose of the direction cosine transformation matrix (DCM) from global to local coordinates given in Eq. (2) in my 2009 Wind Energy Paper: onlinelibrary.wiley.com/doi/abs/10.1002/we.347 (this should be equivalent to the Euler matrix for platform rotations).

Best regards,

Dear Jason,
As for TransMat, I have previously verified that my TransMat matrix gives the same results as the matrix in Eq. 2 of the article you mentioned.

As for your final thoughts, you perfectly understand my goal. I would like to obtain the resulting loads in the 6 DOFs acting in (0,0,0) on the structure, so as to find the accelerations of the platform with respect to the inertial reference system (I want to obtain [‘PtfmTAxi’ ‘PtfmTAyi’ ‘PtfmTAzi’ ‘PtfmRAxi’ ‘PtfmRAyi’ ‘PtfmRAzi’]).
Defined:

  • F_tow = [F, M] the vector of the loads in the 6 DOFs discharged by the turbine, defined in the previous post;
  • F_hydro = [‘HydroFxi’, ‘HydroFyi’, ‘HydroFzi’, ‘HydroMxi’, ‘HydroMyi’, ‘HydroMzi’] the hydrodynamic forces obtained from the OutData of FAST;
  • F_moor = [Fx [1] + Fx [2] + Fx [3] + Fx [4] + Fx [5] + Fx [6], Fy [1] + Fy [2] + Fy [3] + Fy [ 4] + Fy [5] + Fy [6], Fz [1] + Fz [2] + Fz [3] + Fz [4] + Fz [5] + Fz [6], Mx_moor, My_moor, Mz_moor] with Mx_moor, My_moor, Mz_moor moments of transport of forces from Map ++ in (0,0,0), using the variable arms at any time through TransMat as already done for the transport of the turbine loads to move from the nacelle to (0,0,0).
  • F_nonlinear = given by omega x (omega x r_cm * mass) and omega x I dot omega, with omega = [‘PtfmRVxi’, ‘PtfmRVyi’, ‘PtfmRVzi’], r_cm = [x_g, y_g, z_g]instantaneous position of the center of mass with respect to (0,0,0), obtained in the same way in which the arms for the transport moments for turbine and moorings were obtained (using TransMat), mass = total mass and I = 3x3 matrix of inertia, variable according to the inclination of the system
  • F_gravity & restoring = mass * g * [0, 0, 1, y_g, -x_g, 0] where the weight of the system and gravity restoring moments are considered.

Finally, I considered M = mass matrix in (0,0,0) variable as the inclination of the system varies, the correctness of which we have already verified in a previous topic.
To obtain the accelerations a of the platform, I use the formula:
M * a = F_tow + F_moor + F_hydro - F_nonlinear - F_gravity&restoring

Comparing with [‘PtfmTAxi’ ‘PtfmTAyi’ ‘PtfmTAzi’ ‘PtfmRAxi’ ‘PtfmRAyi’ ‘PtfmRAzi’] I found some differences, due according to some analyzes to an incorrect trend of Mx.
So I was wondering if the loads I use are correct or if they are in erroneous reference systems (for example, I’m not sure F_hydro act on platform or is fixed to the inertial frame).

Best regards,
Lorenzo.

Dear Lorenzo,

The total hydrodynamic loads integrated to (0,0,0) translate with the displaced platform, but are oriented in a coordinate system that is parallel to the inertial frame coordinate system.

The rest of your question follows what has been discussed in the following forum topic: http://forums.nrel.gov/t/hydrodynamic-implementation/2284/1. I don’t think I can add more to what has already been discussed there.

Best regards,

Dear Jason
My question is if all the loads reported by me in the previous post are in the same coordinate system parallel to the inertial frame coordinate system, and if multiplying these by inv(M) I can get the accelerations [‘PtfmTAxi’ ‘PtfmTAyi’ ‘PtfmTAzi’ ‘PtfmRAxi’ ‘PtfmRAyi’ ‘PtfmRAzi’]).
Best regards,
Lorenzo

Dear Lorenzo,

Yes, all of these terms appear to be using the same coordinate system that translates, but does not rotate, with the platform.

From the forum topic I linked above, it looked to me that load summation resultant and M*a produced nearly identical results, with some small difference in Mz, which i suspected was caused by numerical round off, considering that the value of Mz was so much less than the others. Numerical round off could be playing a role for small differences in the other components as well.

Best regards,

Dear Jason,
I thank you for the great help you are giving me.
When you talk about our topic linked above, you are referring to the image below, which however refers to a test with CompAero = 0, so we are talking about a test in the absence of aerodynamic loads. In this case, I agree, the error is caused by numerical round off.


By doing a complete offshore test (with aerodynamics as well), I get the following graph (also present in the topic linked above, posted by my colleague Riccardo), where you can see a relevant difference on Mx.

This difference on Mx cause errors in SWAY and ROLL displacements:


I have been looking for the cause of this difference for weeks, but everything looks ok, as you confirmed:

  • I considered all the loads: hydrodynamics, aerodynamics, moorings, gravity & restoring, non linear loads.
  • All loads are applied in the same coordinate system, which moves with the structure but has axes parallel to the inertial reference system.
  • The mass matrix is shifted from the center of mass to the SWL, and rotated then rotated to align it with the axes of the inertial reference system.

Through previous topics, we have been able to confirm all 3 of these points.
However through
a = TotalLoads * inv(M)
I cannot get a = [‘PtfmTAxi’ ‘PtfmTAyi’ ‘PtfmTAzi’ ‘PtfmRAxi’ ‘PtfmRAyi’ ‘PtfmRAzi’] due to the error on Mx, from which the errors in y and rx on displacements (Figure3).
I really can’t find the reason why there is such an error on Mx, while in the other directions the trend is similar to that of FAST. I hope for your intuition.

Best regards,
Lorenzo.

Dear Lorenzo,

If my understanding is correct, everything now matches without aerodynamic loads, but with aerodynamic loads enabled, then Mx is not correct. One problem I can see is that you’ve enabled the generator azimuth DOF in ElastoDyn. This means that the aerodynamic torque (Mx) is not directly transferred to the tower base. Instead, the generator torque (computed by ServoDyn) is what gets transferred to the tower base. Any difference between the generator torque and aerodynamic torque result in acceleration/deceleration of the generator/rotor rather than roll of the platform.

Best regards,

Dear Jason,
You are right: without aerodynamics loads everything matches.
As you observed, by carrying out an onshore test I noticed that using the Control Torque I obtained a right Mx_tow=TwrBsMx (output in OutData), while with RotAeroMx the outputs did not return.
So, I also tried using Mx_nac=CntrTrq. However, both using RotAeroMx and CntrTrq I still have this error on the resulting Mx, the difference is very small between the two cases (as you can see in the following figures, for y and rx).
WhatsApp Image 2020-11-29 at 17.40.33.jpegWhatsApp Image 2020-11-29 at 17.41.19.jpeg
At this point I don’t know if the fact that everything goes back in the absence of aerodynamic loads means that the problem is in the Mx at the base of the tower: in the absence of aerodynamic loads, in fact, there are lower displacements and speeds, which could make less visible a error already present, which is instead enhanced in the presence of aerodynamic loads.
As already said, however, I have hydrodynamic loads and moorings identical to the FAST outputs, the restoring moment is a simple formula, which I have already verified several times, finally the non-linear loads have a very low contribution compared to the other loads.
Excluding these terms, and taking into account that the mass matrix should be correct, the problem if it is not in the aerodynamic loads is in some further contribution to those listed, which however does not seem to me to exist.
Best regards,
Lorenzo

Dear Lorenzo,

What do you mean by CntrTrq, which is not an OpenFAST variable name? I would set Mx_nac = GenTq, which is an output from the ServoDyn module.

Best regards,

Dear Jason,
Yes, I’m referring to GenTq, sorry for the incorrect reference.
Despite this, the error on Mx persists, as mentioned in the previous post.
Best regards,
Lorenzo

Dear Lorenzo,

When you disable the platform DOFs, does your hand calculation of the tower-base moments (based on aerodynamic applied loads, generator torque, and mass/gravity) match the tower-base load calculated by ElastoDyn?

Best regards,

Dear Jason,
Yes, in this case (that is an onshore test) my tower base loads are the same of the outputs TwrBsFx, TwrBsFy, TwrBsFz, TwrBsMx, TwrBsMy, TwrBsMz in the OutData.

My model has the same outputs of FAST both when I do not consider the aerodynamic loads (considering only hydrodynamics, I find same hydrodynamics and Mooring loads and same displacements) and when I do an onshore test (considering only aerodynamics, I find same Generated Power and Tower base loads).
But when I put the two things together, which by themselves work perfectly as FAST, i.e. I do an offshore test, I don’t get the same results as FAST in terms of Mx and so in terms of y and rx.

Best regards,
Lorenzo

Hi Lorenzo,

Have you tried isolating each DOF separately (to avoid complications associated with coupled platform motions)? Do you get the answer you expect in surge only, sway only, heave only, roll only, pitch only, or yaw only?

Best regards,

Dear Jason,
I think I have found the problem, which is the mass matrix.
To obtain the accelerations according to the inertial coordinate system, having the loads oriented in this reference system, the mass matrix must also be referred to the inertial reference system.
Starting from the mass matrix in the center of mass (M_cg), I pass to the mass matrix on the SWL (M_swl), through the procedure already discussed in OC3-Hywind RAOs - #21 by Jason.Jonkman, with
M_swl = inv(TransMat ^ T )* M_cg * inv(TransMat).
In TransMat I entered the instantaneous position of the center of mass in the inertial coordinate system.
Limiting myself to using this transformation alone, which made me obtain the erroneous Mx trend discussed above, I however noticed that for the lower right quadrant (the one containing the inertia matrix) it was wrong, as it was oriented according to the system of axes fixed to the structure, instead of in the inertial reference system.
I therefore thought that, in order to obtain the exact matrix in the inertial frame, it was first necessary to rotate the inertia matrix in the inertial frame and then transport it from the cg to the swl. The equation used was
I_cg (inertial_frame) = R * I_cg (fixed_to_struct_frame) * R ^ T,
with R being the Euler matrix. Obtained the new M_cg matrix, equal to M_cg but with the inertia matrix rotated, I then carried out the transport with the usual formula M_swl = inv(TransMat ^ T) * M_cg * inv(TransMat).
Now I get almost overlapping momentum patterns (Figure 1), but not precise enough to guarantee overlapping y and rx (Figure 2).


I use, for the comparison in Figure 1:

  • Hydrodynamics loads directly from OutData (HydroFx & Co.)
  • Mooring loads directly from OutData (Fx[1] … & Co. forces in x, y, z direction for every line, from which moments Mx, My, Mz for each line is calculated using the instant position of the fairlead basing on the positions of the platform, where the positions are the ones in the OutData)
  • Aerodynamics loads obtained starting from RotAeroFx & Co. in OutData, through the transforms discussed at the beginning of this topic
  • Gravity restoring moment and Non linear loads calculated using positions and speeds in OutData

So, each load used to create the Total loads is a load that is taken directly from the OutData of FAST or derived from the outputs in OutData.
For this reason these differences should not exist, since M*a = Total Loads is an identity.
However, the problems could be 2, always located in the mass matrix:

  • The construction of the mass matrix as described above is incorrect
  • Since I have put Overhang, Precone and other geometric quantities of the turbine to 0, the center of mass and the inertia matrix, and therefore the initial mass matrix may have changed. However I use Mcg and the position of the center of mass given to you in the post linked above, so this could be the error here.

I hope I have clarified the current problem.
Best regards,
Lorenzo.

Dear Lorenzo,

OK, I’m glad to hear that you are honing in on the issue. I agree with your steps, but note that R should equal the transpose of the direction cosine transformation matrix (DCM) from global to local coordinates given in Eq. (2) in my 2009 Wind Energy Paper: onlinelibrary.wiley.com/doi/abs/10.1002/we.347, rather than the Euler matrix (these are equivalent for small rotations).

Best regards,

Dear Jason,
Using the matrix in equation (2) in your 2009 paper instead of the Euler matrix my results correspond to the FAST outputs. I didn’t think this replacement could affect that much. Thank you so much for the help you have given me.
Best regards,
Lorenzo.

A long-time coming; hooray!