Develop Prescribed Platform Motion capability with aero-elasticity

[Amarelio]

Good evening,

I'm trying to implement a capability in OpenFAST that will allow a user to prescribe periodic platform motions to the simulations. The objective is to run these simulations with aero-elasticity so I'm trying to see how to implement this in a fully coupled framework instead of in the AeroDyn module alone. I'll show below where in the code I'm standing now and where I'm thinking of making the changes. If you could give me some guidance, it would be amazing.

As far as my understanding goes, ElastoDyn is the module that receives the loads and calculates the displacements, velocities and accelerations of the system's DOF, so the motion prescription should be done here. I also went through the code starting from the OpenFAST driver (FAST_Prog.f90) and it lead me to ElastoDyn. In the beginning of ElastoDyn.f90 (line 557), there's this:

! SEE IF THESE NEED TO BE CALLED (i.e., if UpdateStates was called, these values are already calculated)
   IF ( UpdateValues ) THEN
         ! Update the OtherState data by calculating the derivative...
      !OtherState%HSSBrTrqC = SIGN( u%HSSBrTrqC, x%QDT(DOF_GeAz) )
      CALL ED_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, ErrStat, ErrMsg ) ! sets m%QD2T = dxdt%QDT
      CALL ED_DestroyContState( dxdt, ErrStat2, ErrMsg2 )
      IF (ErrStat >= AbortErrLev) RETURN
   END IF 

and UpdateValues is set to .TRUE..

At the beginning of the function ED_CalcContStateDeriv (ElastoDyn.f90 line 1825), there's this:

         ! set the coordinate system variables:
      CALL SetCoordSy( t, m%CoordSys, m%RtHS, u%BlPitchCom, p, x, ErrStat2, ErrMsg2 )
         call setErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
         IF (ErrStat >= AbortErrLev) RETURN

      CALL CalculatePositions(        p, x, m%CoordSys,    m%RtHS ) ! calculate positions
      CALL CalculateAngularPosVelPAcc(p, x, m%CoordSys,    m%RtHS ) ! calculate angular positions, velocities, and partial accelerations, including partial angular quantities
      CALL CalculateLinearVelPAcc(    p, x, m%CoordSys,    m%RtHS ) ! calculate linear velocities and partial accelerations
      CALL CalculateForcesMoments(    p, x, m%CoordSys, u, m%RtHS ) ! calculate the forces and moments (requires AeroBladeForces and AeroBladeMoments)            

The 4 functions found in this piece of code (CalculatePositions, CalculateAngularPosVelPAcc, CalculateLinearVelPAcc and CalculateForcesMoments) calculate an object called RtHSdat which seems to contains partial moments, forces and linear and angular accelerations which will be used to calculate the same total quantities in the ED_CalcOutput routine. The RtHSdat object contains the following properties, among others:

REAL(R8Ki) , DIMENSION(1:3)  :: rZ      !< Position vector from inertia frame origin to platform reference (point Z) [m]
REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE  :: PAngVelEX      !< Partial angular velocity (and its 1st time derivative) of the platform (body B) in the inertia frame (body E for earth) [-]
REAL(ReKi) , DIMENSION(1:3)  :: AngVelEX      !< Angular velocity of the platform (body X) in the inertia frame (body E for earth) [-]
REAL(ReKi) , DIMENSION(1:3)  :: LinVelEZ      !< Linear velocity of platform reference (point Z) in the inertia frame [-]

My first idea was to prescribe all the fields in RtHSdat that correspond to the linear and angular positions, velocities and accelerations from the inertia frame origin to platform reference (point Z) but there seem to be some platform motions (including the actual platform DOF like surge and pitch) which relate to the objects x%QT, x%QDT and m%Q2DT as can be seen below (ElastoDyn.f90 line 931):

     ! Platform Motions:

   m%AllOuts( PtfmTDxt) =  DOT_PRODUCT(       m%RtHS%rZ, m%CoordSys%a1 )
   m%AllOuts( PtfmTDyt) = -DOT_PRODUCT(       m%RtHS%rZ, m%CoordSys%a3 )
   m%AllOuts( PtfmTDzt) =  DOT_PRODUCT(       m%RtHS%rZ, m%CoordSys%a2 )
   m%AllOuts( PtfmTDxi) = x%QT  (DOF_Sg )
   m%AllOuts( PtfmTDyi) = x%QT  (DOF_Sw )
   m%AllOuts( PtfmTDzi) = x%QT  (DOF_Hv )
   m%AllOuts( PtfmTVxt) =  DOT_PRODUCT( m%RtHS%LinVelEZ, m%CoordSys%a1 )
   m%AllOuts( PtfmTVyt) = -DOT_PRODUCT( m%RtHS%LinVelEZ, m%CoordSys%a3 )
   m%AllOuts( PtfmTVzt) =  DOT_PRODUCT( m%RtHS%LinVelEZ, m%CoordSys%a2 )
   m%AllOuts( PtfmTVxi) = x%QDT (DOF_Sg )
   m%AllOuts( PtfmTVyi) = x%QDT (DOF_Sw )
   m%AllOuts( PtfmTVzi) = x%QDT (DOF_Hv )
   m%AllOuts( PtfmTAxt) =  DOT_PRODUCT(                 LinAccEZ, m%CoordSys%a1 )
   m%AllOuts( PtfmTAyt) = -DOT_PRODUCT(                 LinAccEZ, m%CoordSys%a3 )
   m%AllOuts( PtfmTAzt) =  DOT_PRODUCT(                 LinAccEZ, m%CoordSys%a2 )
   m%AllOuts( PtfmTAxi) = m%QD2T(DOF_Sg  )
   m%AllOuts( PtfmTAyi) = m%QD2T(DOF_Sw  )
   m%AllOuts( PtfmTAzi) = m%QD2T(DOF_Hv  )
   m%AllOuts( PtfmRDxi) = x%QT  (DOF_R )*R2D
   m%AllOuts( PtfmRDyi) = x%QT  (DOF_P )*R2D
   m%AllOuts( PtfmRDzi) = x%QT  (DOF_Y )*R2D
   m%AllOuts( PtfmRVxt) =  DOT_PRODUCT( m%RtHS%AngVelEX, m%CoordSys%a1 )*R2D
   m%AllOuts( PtfmRVyt) = -DOT_PRODUCT( m%RtHS%AngVelEX, m%CoordSys%a3 )*R2D
   m%AllOuts( PtfmRVzt) =  DOT_PRODUCT( m%RtHS%AngVelEX, m%CoordSys%a2 )*R2D
   m%AllOuts( PtfmRVxi) = x%QDT (DOF_R )*R2D
   m%AllOuts( PtfmRVyi) = x%QDT (DOF_P )*R2D
   m%AllOuts( PtfmRVzi) = x%QDT (DOF_Y )*R2D
   m%AllOuts( PtfmRAxt) =  DOT_PRODUCT(                 AngAccEX, m%CoordSys%a1 )*R2D
   m%AllOuts( PtfmRAyt) = -DOT_PRODUCT(                 AngAccEX, m%CoordSys%a3 )*R2D
   m%AllOuts( PtfmRAzt) =  DOT_PRODUCT(                 AngAccEX, m%CoordSys%a2 )*R2D
   m%AllOuts( PtfmRAxi) = m%QD2T(DOF_R )*R2D
   m%AllOuts( PtfmRAyi) = m%QD2T(DOF_P )*R2D
   m%AllOuts( PtfmRAzi) = m%QD2T(DOF_Y )*R2D

Also, the RtHSdat properties of interest are calculated from x%QT, x%QDT and m%Q2DT and partial angular and linear velocities generally speaking. Now, I think I should prescribed x%QT, x%QDT and m%Q2DT and partial angular and linear velocities that are used to calculate the RtHSdat properties from the inertia frame origin to platform reference (point Z).

Is this reasoning correct or am I completely off of the tracks? If not, I think x%QT, x%QDT and m%Q2DT are straightforward to prescribe since they seem the position, velocity and acceleration of the typical platform degrees of freedom corrected by a constant (which I don't know what it is) but I don't understand what are the partial angular and linear velocities.

Best regards,
Ricardo Amaral

[jjonkman]

Dear Ricardo,

The platform motions are structural DOFs of ElastoDyn. As such, they are solved by calculating F=m*a for the acceleration (a, the unknowns) based on applied loads (F, the knowns) and time integrating. I would generally not recommend trying to prescribe some of the states while calculating others, as, this would likely require a rearrangement of the formulation in terms of which terms are knowns (a for the platform, F for other components) and which are unknowns (F for the platform, a for the other components).

Instead, I would recommend applying loads at the platform reference point that result in the desired motion of the platform. The user-specified external platform called via CompSub = 2 in the OpenFAST primary input file provides a direct means for setting the platform loads.

Similar topics have been discussed a number of times on our forum, e.g. see: https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=2259&p=12800. I suggest review these for further information.

Best regards,

[cc200528]

Dear jjonkman,
Now I cannot directly locate from https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=2259&p=12800 to discussion, could you please provide a new discussion link.

Best regards

[jjonkman]

Dear @cc200528,

Unfortunately, the deep links broke when we switched from the old wind.nrel.gov forum to the new forums.nrel.gov forum, so, I'm not exactly sure which forum topic I linked to previously. But from a quick search on the new forum, I see quite a few related topics: https://forums.nrel.gov/search?q=platform%20motion%20prescribe%20order%3Alatest.

Best regards,

[Amarelio]

Dear jjonkman,

Thank you for the explanation, I understand. On that case, according to what I read about the ExtPtfm module and the Craig-Bampton reduction, I first need to produce the reduced matrices (M,C and K) and reduced loads for the floater using a high-fidelity FEM software.

1. Assuming that I want to prescribe a surge motion in OpenFAST, do I need to then apply loads in the surge direction on the FEM software? If I want to prescribe another type of motion, for instance pitch, do I need to recalculate the reduced matrices?

2. If I still want to prescribed a surge motion but with a different amplitude and since the system is linearized, is it sufficient to multiply the reduced loads in the ExtPtfm input by a given constant (with no change in the reduced matrices) or do I need to recalculate the reduced matrices and loads?

Best regards,
Ricardo Amaral

[jjonkman]

Dear Ricardo,

Actually, I would forget the Craig-Bampton reduction.

Instead, I would set the externally applied force in ExtPtfm equal to

F = K*( DesiredDisplacement - X ) + C*( DesiredVelocity - Xdot )

This is effectively what the Seismic module of FAST v7 does to impose earthquake excitation of the platform. See equations (1)-(3) in the Seismic documentation for more information: https://drive.google.com/file/d/1YbZgeOZ46R5IiMOT-kpHtoElOfyBlwUG/view?usp=sharing.

Best regards,

[Amarelio]

Dear jjonkman,

I think I understood the force that but I'm not sure regarding the M, C and K matrix inputs. Would these be 6x6 matrices considering the 6 platform degrees of freedom where:

• M = Mptfm + Madded(f = f_imposed_motion), where Mptfm comes from ElastoDyn's inertia and Madded from the .1 file from the HydroDyn inputs. I also read that this inertia should be increased so that the applied loads dominate other loads on the platform;
• C(f = f_imposed_motion), coming from the .1 file from HydroDyn inputs;
• K = Khst + Kmooring, where Khst comes from the HydroDyn's .hst file. Kmooring I wouldn't know how to obtain.

At the same time, I have the feeling that the actual matrix values don't matter as long as they are in accordance with the coefficients k and c used to calculate the prescribed forces. At least looking at the equation of motion for pure surge [mxdd + cxd + kx = F = k(xtarget - x) + c(xdtarget - xd)], I can see the terms cancelling on the left and right sides, is this correct?
In this case:

• For a pure degree of freedom, the frequency of the motion is tuned by the m, c and k parameters which are related by equations (1) and (2) of the Seismic module documentation?
• If 2 degrees of freedom are prescribed, then the coupling terms need to be taken into account in the prescribed force?

And finally, on the ExtPtfm input file:

Comment describing the model
---------------------- SIMULATION CONTROL --------------------------------------
False          Echo         - Echo input data to <RootName>.ech (flag)
"default"     DT           - Communication interval for controllers (s) (or "default")
   3          IntMethod    - Integration Method {1:RK4; 2:AB4, 3:ABM4} (switch)
---------------------- REDUCTION INPUTS ----------------------------------------
   1          FileFormat    - File Format {0:Guyan; 1:FlexASCII} (switch)
"ExtPtfm_SE.dat"   Red_FileName    - Path of the file containing Guyan/Craig-Bampton inputs (-)
"NA"         RedCst_FileName - Path of the file containing Guyan/Craig-Bampton constant inputs (-) (currently unused)
-1           NActiveDOFList - Number of active CB mode listed in ActiveDOFList, use -1 for all modes (integer)
1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25     ActiveDOFList  - List of CB modes index that are active, [unused if NActiveDOFList<=0]
0          NInitPosList   - Number of initial positions listed in InitPosList, using 0 implies all DOF initialized to 0  (integer)
0,         InitPosList    - List of initial positions for the CB modes  [unused if NInitPosList<=0 or EquilStart=True]
0          NInitVelList   - Number of initial positions listed in InitVelList, using 0 implies all DOF initialized to 0  (integer)
0,         InitVelList    - List of initial velocities for the CB modes  [unused if NInitVelPosList<=0 or EquilStart=True]
---------------------- OUTPUT --------------------------------------------------
True          SumPrint      - Print summary data to <RootName>.sum (flag)
         1    OutFile      - Switch to determine where output will be placed: {1: in module output file only; 2: in glue code output file only; 3: both} (currently unused)
True          TabDelim     - Use tab delimiters in text tabular output file? (flag) (currently unused)
"ES10.3E2"    OutFmt       - Format used for text tabular output (except time).  Resulting field should be 10 characters. (quoted string) (currently unused)
          0   TStart       - Time to begin tabular output (s) (currently unused)
              OutList      - The next line(s) contains a list of output parameters.  See OutListParameters.xlsx for a listing of available output channels, (-)

Should I put the ActiveDOFList to 0? Is there something else I need to change?

Best regards,
Ricardo Amaral

[jjonkman]

Dear Ricardo,

My understanding is that you want to prescribe the motion of the platform. As such, you don't need to consider any physically correct values for platform mass/inertia, hydrodynamics, or moorings. Instead, as discussed in the Seismic documentation, you should specify an artificially high platform mass/inertia, and choose a stiffness, damping, and force that results in the desired motion.

In your equation of motion, you've double booked the c*xd and k*x terms. Instead, the single DOF equation of motion should be:
m*xdd+c*xd+k*x = c*xdtarget + k*xtarget.

(This is equivalent to m*xdd = c*(xdtarget - xd) + k*(xtarget - x).)

The ExtPtfm module allows you to specify M, C, K, and F for a rigid platform by setting FileFormat = 0 (Guyan). I would choose M such that when it is added to the physical mass/inertia of the wind turbine atop the platform that M dominates. C and K can be derived from M as discussed in the Seismic documentation. Then, F includes the terms that are left on the right-hand side: F = C*xdtarget + K*xtarget. If the mass added via ExtPtfm (M) is large relative to the physical/mass inertia of the wind turbine, then I would expect you could get by with making M, C, and K diagonal matrices (the Seismic module didn't have to worry about off-diagonal terms because Seismic only applied to the platform translation DOFs).

I hope that helps.

Best regards,

[Amarelio]

Dear jjonkman,

I'm coming back to you to report my progress on this. After prescribing the force indicated by the formula F = Cxdtarget + Kxtarget considering: the relation between M, K and C found in the Seismic module documentation, a pure surge motion with Amp = 2 m and f = 0.05 Hz, MoorDyn turned off, the remaining 5 DOF's turned off in ElastoDyn and a prescribed mass in the surge DOF's of 100 times the platform's mass found in the ElastoDyn input, I found out that the motion's amplitude was overpredicted by around 25%. This can be seen below on a prescribed mass sensitivity study I perfomed to see when would the motion converge to the prescribed one. For a prescribed mass of 100 times the platform mass, one can see that the motion has converged by comparison with the purple curve. "PM" is the prescribed mass and "TM" is the turbine's platform mass:


Facing this problem, I changed my approach and decided to prescribe a force based on a single DOF system's response to a constant amplitude force whose formula is depicted below, where F is the prescribed force's amplitude, omega is the force's frequency and c and k are the system's damping and stiffness in the DOF:

(This formula can be found in books on vibrations and system dynamics).
Considering this, one can see that prescribing the force below will result in a motion with the correct amplitude and the same frequency as the force:

This made sense to me because, for such a high prescribed mass and stiffness and damping along with it, the scale of the turbine originated perturbations would be so low compared to the ones prescribed that the system would effectively be in the single DOF approximation which, I think, is why you told me to prescribe such a dominant prescribed mass. After this, I managed to get the correct amplitude of 2 m for prescribed masses which were 10 or 100 times that of the platform's mass:

This force prescription formula was tested for the remaining pure DOF and also for surge with other amplitudes and frequencies and the obtained motion was correct.

I may have missed something in your explanation or made a mistake while launching the first simulations but in any case, what I did seems to have worked so I thought it would be of value to report it back to this repository.

Best regards,
Ricardo Amaral

[jjonkman]

Dear Ricardo,

I'm not following everything you are saying, but from my original description, I was suggesting that the mass in ExtPtfm be set much larger than the full system mass (actual platform mass + tower mass + nacelle mass + rotor mass). The stiffness and damping can be derived from that mass and an actuator frequency that is higher than the highest natural frequency of the turbine model when the platform is fixed. With these values for M, C, and K, you should get the desired motion to match any target motion you specify.

Best regards,