Issues related to the Blade Element Momentum (BEM) theory

[Yingxin.Lv]

Hi everyone
I conducted a dynamic analysis of a 5MW fixed-bottom wind turbine. I applied the Blade Element Momentum Theory to calculate the aerodynamic forces of the wind turbine, with the wind speed set as a steady-state wind profile that accounts for shear variations along the height.
I perform computations for each time step using the following algorithm（No correction models were considered in the calculations.During the simulation, the rotor speed was maintained constant, and servo control systemswere disabled）) :

All parameters in the formulas are sourced from the OpenFAST input files(NRELOffshrBsline5MW_AeroDyn_blade.dat).,“The wind speed V0 was calculated via [v₀ᵢ = HWindSpeed× [(RefHt + (r_hub + r_blade) × cos(Azimuth)) / RefHt]^PLExp， i=1,2,3] and validated against the output variable VUndx.” However, there are certain discrepancies between the calculated Pn​ (normal force) and Pt​ (tangential force) and the OpenFAST output variable ABaN00bFx``ABaN00bFy" ,particularly pronounced within the initial few seconds of the simulation."
I am deeply puzzled about where I made errors in my calculations.

I would appreciate any thoughts or insights.
Thank you and best regards,

------- AERODYN v15 for OpenFAST INPUT FILE -----------------------------------------------
NREL 5.0 MW offshore baseline aerodynamic input properties.
====== General Options ============================================================================
False Echo - Echo the input to “.AD.ech”? (flag)
“default” DTAero - Time interval for aerodynamic calculations {or “default”} (s)
1 WakeMod - Type of wake/induction model (switch) {0=none, 1=BEMT, 2=DBEMT, 3=OLAF} [WakeMod cannot be 2 or 3 when linearizing]
1 AFAeroMod - Type of blade airfoil aerodynamics model (switch) {1=steady model, 2=Beddoes-Leishman unsteady model} [AFAeroMod must be 1 when linearizing]
0 TwrPotent - Type tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}
0 TwrShadow - Calculate tower influence on wind based on downstream tower shadow? (switch) {0=none, 1=Powles model, 2=Eames model}
False TwrAero - Calculate tower aerodynamic loads? (flag)
False FrozenWake - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing]
False CavitCheck - Perform cavitation check? (flag) [AFAeroMod must be 1 when CavitCheck=true]
False Buoyancy - Include buoyancy effects? (flag)
False CompAA - Flag to compute AeroAcoustics calculation [used only when WakeMod = 1 or 2]
“unused” AA_InputFile - AeroAcoustics input file [used only when CompAA=true]
====== Environmental Conditions ===================================================================
“default” AirDens - Air density (kg/m^3)
“default” KinVisc - Kinematic viscosity of working fluid (m^2/s)
“default” SpdSound - Speed of sound in working fluid (m/s)
“default” Patm - Atmospheric pressure ¶ [used only when CavitCheck=True]
“default” Pvap - Vapour pressure of working fluid ¶ [used only when CavitCheck=True]
====== Blade-Element/Momentum Theory Options ====================================================== [unused when WakeMod=0 or 3]
1 SkewMod - Type of skewed-wake correction model (switch) {1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0 or 3]
“default” SkewModFactor - Constant used in Pitt/Peters skewed wake model {or “default” is 15/32*pi} (-) [used only when SkewMod=2; unused when WakeMod=0 or 3]
False TipLoss - Use the Prandtl tip-loss model? (flag) [unused when WakeMod=0 or 3]
False HubLoss - Use the Prandtl hub-loss model? (flag) [unused when WakeMod=0 or 3]
TRUE TanInd - Include tangential induction in BEMT calculations? (flag) [unused when WakeMod=0 or 3]
False AIDrag - Include the drag term in the axial-induction calculation? (flag) [unused when WakeMod=0 or 3]
False TIDrag - Include the drag term in the tangential-induction calculation? (flag) [unused when WakeMod=0,3 or TanInd=FALSE]
“Default” IndToler - Convergence tolerance for BEMT nonlinear solve residual equation {or “default”} (-) [unused when WakeMod=0 or 3]
100 MaxIter - Maximum number of iteration steps (-) [unused when WakeMod=0]
====== Dynamic Blade-Element/Momentum Theory Options ============================================== [used only when WakeMod=2]
1 DBEMT_Mod - Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1, 3=constant tau1 with continuous formulation} (-) [used only when WakeMod=2]
4 tau1_const - Time constant for DBEMT (s) [used only when WakeMod=2 and DBEMT_Mod=1 or 3]
====== OLAF – cOnvecting LAgrangian Filaments (Free Vortex Wake) Theory Options ================== [used only when WakeMod=3]
“unused” OLAFInputFileName - Input file for OLAF [used only when WakeMod=3]
====== Beddoes-Leishman Unsteady Airfoil Aerodynamics Options ===================================== [used only when AFAeroMod=2]
3 UAMod - Unsteady Aero Model Switch (switch) {2=B-L Gonzalez, 3=B-L Minnema/Pierce, 4=B-L HGM 4-states, 5=B-L HGM+vortex 5 states, 6=Oye, 7=Boeing-Vertol} [used only when AFAeroMod=2]
True FLookup - Flag to indicate whether a lookup for f’ will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files (flag) [used only when AFAeroMod=2]
0 UAStartRad - Starting radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAStartRad=0]
1 UAEndRad - Ending radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAEndRad=1]
====== Airfoil Information =========================================================================
1 AFTabMod - Interpolation method for multiple airfoil tables {1=1D interpolation on AoA (first table only); 2=2D interpolation on AoA and Re; 3=2D interpolation on AoA and UserProp} (-)
1 InCol_Alfa - The column in the airfoil tables that contains the angle of attack (-)
2 InCol_Cl - The column in the airfoil tables that contains the lift coefficient (-)
3 InCol_Cd - The column in the airfoil tables that contains the drag coefficient (-)
4 InCol_Cm - The column in the airfoil tables that contains the pitching-moment coefficient; use zero if there is no Cm column (-)
0 InCol_Cpmin - The column in the airfoil tables that contains the Cpmin coefficient; use zero if there is no Cpmin column (-)
8 NumAFfiles - Number of airfoil files used (-)
“../5MW_Baseline/Airfoils/Cylinder1.dat” AFNames - Airfoil file names (NumAFfiles lines) (quoted strings)
“../5MW_Baseline/Airfoils/Cylinder2.dat”
“../5MW_Baseline/Airfoils/DU40_A17.dat”
“../5MW_Baseline/Airfoils/DU35_A17.dat”
“../5MW_Baseline/Airfoils/DU30_A17.dat”
“../5MW_Baseline/Airfoils/DU25_A17.dat”
“../5MW_Baseline/Airfoils/DU21_A17.dat”
“../5MW_Baseline/Airfoils/NACA64_A17.dat”
====== Rotor/Blade Properties =====================================================================
True UseBlCm - Include aerodynamic pitching moment in calculations? (flag)
“../5MW_Baseline/NRELOffshrBsline5MW_AeroDyn_blade.dat” ADBlFile(1) - Name of file containing distributed aerodynamic properties for Blade #1 (-)
“../5MW_Baseline/NRELOffshrBsline5MW_AeroDyn_blade.dat” ADBlFile(2) - Name of file containing distributed aerodynamic properties for Blade #2 (-) [unused if NumBl < 2]
“../5MW_Baseline/NRELOffshrBsline5MW_AeroDyn_blade.dat” ADBlFile(3) - Name of file containing distributed aerodynamic properties for Blade #3 (-) [unused if NumBl < 3]
====== Hub Properties ============================================================================== [used only when Buoyancy=True]
0 VolHub - Hub volume (m^3)
0 HubCenBx - Hub center of buoyancy x direction offset (m)
====== Nacelle Properties ========================================================================== [used only when Buoyancy=True]
0 VolNac - Nacelle volume (m^3)
0, 0, 0 NacCenB - Position of nacelle center of buoyancy from yaw bearing in nacelle coordinates (m)
====== Tail Fin Aerodynamics =======================================================================
False TFinAero - Calculate tail fin aerodynamics model (flag)
“unused” TFinFile - Input file for tail fin aerodynamics [used only when TFinAero=True]
====== Tower Influence and Aerodynamics ============================================================ [used only when TwrPotent/=0, TwrShadow/=0, TwrAero=True, or Buoyancy=True]
12 NumTwrNds - Number of tower nodes used in the analysis (-) [used only when TwrPotent/=0, TwrShadow/=0, TwrAero=True, or Buoyancy=True]
TwrElev TwrDiam TwrCd TwrTI TwrCb ! TwrTI used only when TwrShadow=2; TwrCb used only when Buoyancy=True
(m) (m) (-) (-) (-)
0.0000000E+00 6.0000000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
8.5261000E+00 5.7870000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
1.7053000E+01 5.5740000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
2.5579000E+01 5.3610000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
3.4105000E+01 5.1480000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
4.2633000E+01 4.9350000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
5.1158000E+01 4.7220000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
5.9685000E+01 4.5090000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
6.8211000E+01 4.2960000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
7.6738000E+01 4.0830000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
8.5268000E+01 3.8700000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
8.7600000E+01 3.8700000E+00 1.0000000E+00 1.0000000E-01 0.0000000E+00
====== Outputs ====================================================================================
True SumPrint - Generate a summary file listing input options and interpolated properties to “.AD.sum”? (flag)
0 NBlOuts - Number of blade node outputs [0 - 9] (-)
1, 9, 19 BlOutNd - Blade nodes whose values will be output (-)
0 NTwOuts - Number of tower node outputs [0 - 9] (-)
1, 2, 6 TwOutNd - Tower nodes whose values will be output (-)
OutList - The next line(s) contains a list of output parameters. See OutListParameters.xlsx for a listing of available output channels, (-)
“RtArea”
“B1N3Clrnc, B2N3Clrnc, B3N3Clrnc”
“RtAeroFxh”
“RtAeroFyh”
“RtAeroFzh”
“RtAeroMxh”
“RtAeroMyh”
“RtAeroMzh”
“RtAeroFxi”
“RtAeroFyi”
“RtAeroFzi”
“RtAeroMxi”
“RtAeroMyi”
“RtAeroMzi”
END of input file (the word “END” must appear in the first 3 columns of this last OutList line)
====== Outputs for all blade stations (same ending as above for B1N1… =========================== [optional section]
3 BldNd_BladesOut - Number of blades to output all node information at. Up to number of blades on turbine. (-)
50 BldNd_BlOutNd - Future feature will allow selecting a portion of the nodes to output. Not implemented yet. (-)
OutList - The next line(s) contains a list of output parameters. See OutListParameters.xlsx for a listing of available output channels, (-)
“Fx”
“Fy”
“VUndy”
“VUndx”
“VUndz”
“TnInd”
“AxInd”
“Alpha”
“Theta”
“Phi”
“Cl”
“Cd”
“Cx”
“Cy”
END (the word “END” must appear in the first 3 columns of this last OutList line in the optional nodal output section)

---

[Riad.Elhamoud]

Dear @Yingxin.Lv,
I am Riad, wind turbine enthusiat.
I have a few comments:
1- The algorithm you implement has no guaranteed convergence. Meaning that, this algorithm can spend hours searching for the axial and tangential induction factors without giving you any solution.
I suggest you use another variant of the BEMT. This variant is proposed by Professor Ning. In a glimpse, this method searches for the inflow angle instead of searching for the axial and tangential induction factors. You can find the paper here:
https://onlinelibrary.wiley.com/doi/abs/10.1002/we.1636
You dont need to read the whole paper unless you are interested in the mathematics behind. Just simply grab the algorithm proposed at the end of the paper.

2- Does your model takes into account blade flexibility ? I mean in your model, the blades are considered rigid or flexible ? Whatever the case, you should have same considerations with OpenFAST.

3- I know that in the basic BEMT, only the rotor is modeled. This means that tower, foundation and blades degrees of freedom should be disabled. Do u agree ?

4- How to force OpenFAST to have a constant rotor speed ? I mean that in all textbooks, when implementing BEMT, constant rotor speed is considered. How did you achieve this in OpenFAST in order to have a fair comparison with your model ? What about the unsteady aerodynamics model in OpenFAST? Is it disabled ?

5- I suggest also to compare the axial and tangential induction factors between your model and OpenFAST.

Hope that helps.

Best Regards,

Riad

[Yingxin.Lv]

Dear @ Riad.Elhamoud,
I sincerely appreciate your thoughtful reply. The solution you provided was clear and effective. Thanks again!

1. I’ve already downloaded the paper you recommended. Thank you for pointing me towards this valuable resource!
2. 3 .5
   I fully agree with your insights. After disabling the degrees of freedom for the tower and blades in ElastoDyn, I was pleased to find that my results closely approximated those of OpenFAST. I sincerely appreciate your help in guiding me toward this validation.
   
   May I kindly inquire about how I might modify and enhance the BEM algorithm I am using to account for tower and blade flexibility? Are you familiar with any theoretical frameworks or methodologies in this field, or could you suggest any relevant books or papers that address this integration?
   
   Additionally, in OpenFAST, I was wondering if there is an option to enable tower and blade degrees of freedom in ElastoDyn while excluding the influence of their flexibility during aerodynamic force calculations. Your expertise on this would be invaluable.
   
   Thank you again for your generous contributions to my understanding. I would greatly appreciate any further guidance you might provide.
3. Regarding the method to set the rotational speed as a constant in OpenFAST, you can find the discussion in this link Usage question: fixed pitch & rotor speed · Issue #1035 · OpenFAST/openfast.
   Best Regards,
   yingxin

[Riad.Elhamoud]

Dear @Yingxin.Lv,

I am glad that my comments helped you !

As i know OpenFAST as a user, OpenFAST always takes into account the flexibilities of the tower and blades when computing aerodynamic forces since OpenFAST performs aero-elastic computation. Maybe @Jason.Jonkman could confirm what i am saying.

Regarding the consideration of blade deformation into BEMT, i think (but i am not quite sure) that the blade deformation is considered in the BEMT through relative wind speed. You know that when computing aerodynamic loads, one should always use the relative wind speed, so i think that through it, one can consider the aero-elastic computation. If what i am saying is right, you should pay attention to the different transformations you should perform since the degrres of freedom of the blade are defined in blade coordinate system denoted (b in FAST_v7 user manual). So, i think sevral transformations should be done since there are a lot of frames.
But i repeat, this is my opinion, i am not quite sure about it.

Thank you for the link you have provided to me. I will certainly go through it.

Hope that helps.

Best Regards,

Riad

[Jason.Jonkman]

Dear @Riad.Elhamoud and @Yingxin.Lv,

I generally agree with this discussion and don’t have much to add. I confirm that OpenFAST accounts for the structural flexiblities in the aeordynamic calculation, including taking the displaced position and relative velocity (wind minus structural velocity) into account in the aerodynamics calculations. For the latter, it is important to take into consideration the coordinate systems used for different calculations, e.g., as documented for AeroDyn here: 4.2.2. Coordinate systems — OpenFAST v4.0.4 documentation.

Best regards,

[Yingxin.Lv]

Dear @Jason.Jonkman,

I want to understand the relationship between the in-plane/out-of-plane primary responses of the blade and the tower, and the output variables STVx/STVy/STVz in AeroDyn**.** For example, the description of STVx is "x-component of structural translational velocity at each node (local blade coordinate system).

First, I computed the simplest case scenario: I set both the fore-aft tower bending-mode DOF and side-to-side tower bending-mode DOF to False. In the results, the velocity at any blade node in the x₁-x₂-x₃ coordinate system (shown in the figure below) should be [OoPDefl1_dot, IPDefl1_dot - ω·r, 0] (assuming no axial blade deformation).

Next, to transform this vector into the local blade coordinate system (also shown in the figure below), the x₁-x₂-x₃ system requires two angular transformations. These two angles should represent the torsional rotations of the blade about the x₁ and y₁ axes at point *r*, which can be obtained from the outputs Spn1RDxb1 and Spn1RDyb1 in ElastoDyn**.**.

The transformation matrices are:

TRANSM1 = [1, 0, 0;
0, cos(Spn1RDxb1), -sin(Spn1RDxb1);
0, sin(Spn1RDxb1), cos(Spn1RDxb1)];

TRANSM2 = [cos(Spn1RDyb1), 0, sin(Spn1RDyb1);
0, 1, 0;
-sin(Spn1RDyb1), 0, cos(Spn1RDyb1)];

vb(:,) = (TRANSM2)^⁻¹ * (TRANSM1)^⁻¹ * [OoPDefl1_dot; IPDefl1_dot - ω·r; 0];

However, The results derived from [OoPDefl1_dot, IPDefl1_dot - ω·r, 0] do not match the y-component (STVy) and z-component (STVz) of structural translational velocity at each node in the local blade coordinate system. Could you identify where my error lies?

I would appreciate any thoughts or insights.
Thank you and best regards,

Best Regards,

[Jason.Jonkman]

Dear @Yingxin.Lv,

Just a few comments:

• OoPDefl1_dot and IPDefl1_dot are not outputs of ElastoDyn; how are you getting these?
• Your coordinate transformations involving TRANSM1 and TRANSM2 are not identical to what ElastoDyn is using. The coordinate systems used by ElastoDyn are documented here: https://openfast.readthedocs.io/en/main/\_downloads/c9ae932057ea1e64c79154b36d3ab30a/FASTCoordinateSystems.doc.
• AeroDyn outputs STVx, STVy and STVz use the legacy local output system documented here: 4.2.2. Coordinate systems — OpenFAST v4.1.2 documentation, and the plan is to deprecate them from a future release of AeroDyn.

Best regards,

[Yingxin.Lv]

Dear @Jason.Jonkman ,
I have a question about the coordinate systems in OpenFAST. In the document “coordinate systems used by ElastoDyn”, Does the Blade Element-Fixed Coordinate System refer to the same coordinate system as the local blade coordinate system (which is used for many output options in AeroDyn_nodes)?

In the outputs of AeroDyn_nodes, Fx and Fy are in the x-y coordinate system shown in the figure below, That is, the Airfoil system defined in ‘’Section 4.2.2. Coordinate Systems of the OpenFAST v4.1.2 documentation‘’.
What is the relationship between the x-y coordinate system shown in the figure, the Blade Element-Fixed Coordinate System, and the local blade coordinate system? Are they the same coordinate system, or is a coordinate transformation required?

I would appreciate any thoughts or insights.
Thank you and best regards,

[Jason.Jonkman]

Dear @Yingxin.Lv,

The ElastoDyn documentation you are referring is slightly outdated in the sense that the “m” and “te” coordinate systems are not used by OpenFAST (they were used in FAST v7).

Within OpenFAST, for understanding AeroDyn’s coordinate systems, I suggest reviewing the documentation here: 4.2.2. Coordinate systems — OpenFAST v4.1.2 documentation.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

I am currently developing a flywheel (FW) controller for tower fore-aft damping. The controller is based on the global nacelle IMU fore-aft velocity, while the flywheel-induced Coriolis forces are generated individually within each blade.

To better understand the blade-local kinematics, I examined the AeroDyn nodal output VXs, which appears to represent the local out-of-plane velocity of the blade nodes. When plotting VXs for all three blades, I observe the expected 120° phase shift between the blades.

However, I also noticed that, for each individual blade, the VXs component changes sign as the blade passes approximately 180° azimuth. In other words, while a blade is in the 0°–180° azimuth region, the local VXs is positive, and after crossing the opposite side of the rotor (around 180° azimuth), it becomes negative. I observe a corresponding sign change in the calculated blade-local Coriolis force as well.

My question is:

Is this 180° sign reversal of the blade-local VXs an expected consequence of AeroDyn’s rotating local coordinate system and blade kinematics, or is it related to a modal characteristic of the structural model?

Is there any documentation or reference that explains this behaviour?

I would appreciate any insights or references that could help me better understand this phenomenon.

Regards,

[Jason.Jonkman]

Dear @Abhinay.Goga,

I’m not sure which specific azimuth angle output you referring (different azimuth outputs have different references regarding what constitutes “zero”). Regardless, it sounds like you are seeing a dominate once-per revolution oscillation of the out-of-plane mode, which is not surprising given the presence of gravity, shear, and/or skew of the flow.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

Thank you for your response.

Up to now, I have assumed that the ElastoDyn output Azimuth represents the global azimuth angle of Blade 1, where 0° corresponds to the blade pointing vertically upward and 180° corresponds to the blade pointing vertically downward.

Following your comment, I compared the individual blade azimuth outputs from AeroDyn with the global azimuth output from ElastoDyn. It appears that there is an approximate 70° offset between the AeroDyn Blade 1 azimuth and the ElastoDyn global azimuth.

I have calculated the statistics fro a 600 sec Simulation and it follows:

Offset Mean deg Median deg Std deg Min deg Max deg

B1Azimuth - Azimuth 90.00 90.19 16.19 34.01 149.01
B2Azimuth - Azimuth -149.99 -148.49 16.19 -179.99 179.99
B3Azimuth - Azimuth -29.99 -29.80 16.19 -85.98 29.01

So there is a 90° offset between AeroDyn and ElastoDyn result?

Regards,

[Jason.Jonkman]

Dear @Abhinay.Goga,

The ElastoDyn Azimuth output is defined relative to ElastoDyn input AzimB1Up, which sounds like it is set to 0deg for your model.

In AeroDyn, the local blade azimuth angle (B1Azimuth, etc.) is defined based on how the rotor disk-averaged relative velocity projects onto the disk (based on the skew angle between the rotor axis and inflow) rather than some rigid-body rotation of the rotor. This is discussed more in the following forum topic: Azimuth Angle of the Blade.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

After reading through the thread, I realized that using the ElastoDyn azimuth and adding 120° offsets for the following blades is a reasonable approach when treating the blades as rigid bodies. I am sharing my current issue here because I could not find any existing thread discussing something similar.

First, a short description of my current control strategy may provide some context.

An increase in tower fore-aft (FA) velocity, in either direction, can be damped by applying an opposing force from the rotor plane. If the upper and lower blades produce equal Coriolis forces in opposite directions, they generate a net moment at the hub that counteracts the tower FA moment. Alternatively, one can charge the upper-half blades while discharging the lower-half blades (or vice versa) so that the Coriolis forces act in the same direction, causing the blades to flap against the tower FA motion.

My controller therefore works as follows: when the tower FA speed increases (i.e., the tower is moving further away from its equilibrium position), I apply a damping action. When the tower starts returning toward its original position, I reverse the FW action and inject energy again so that the tower is pushed back in the thrust direction. The intention is not to remove the oscillation completely, but to reduce its oscillation amplitude over time, whle mantaining Fluid charge states.

The issue is that my controller works consistently only when I reverse the FW action for all three blades based on one half-rotation of the rotor, i.e.

if sin_bl1 > 0
    % FW action
else
    % Reverse FW action
end

If not, I see excitations instead of damping

In other words, flipping the sign based on Blade 1 (or equivalently the rotor half revolution) produces stable damping. However, I cannot find any physical explanation for this apparent tower symmetry.

I have also tried using the blade out-of-plane motion as an additional feedback signal. The FW model itself is implemented inside a modified BeamDyn model, meaning that the Coriolis forces are computed directly from the BeamDyn blade motion.

Since my current objective is fore-aft damping, I considered the BeamDyn out-of-plane rotational velocities (RVYg) at different blade nodes. I use the difference between the root and approximately 70% span as a measure of the blade out-of-plane angular deflection speed (deg/s).

My reasoning is the following. If the out-of-plane speed difference increases, the blade is deflecting either in the same direction as the tower FA motion or in the opposite direction. Based on the sign convention, if the blade is already deflecting together with the tower and the tip is moving faster than the root, I can apply a larger fluid speed to generate a stronger Coriolis damping force. If the relative blade deflection is smaller, I reduce the commanded fluid speed. Conversely, if the blade is deflecting opposite to the tower motion, I apply no FW action, because the blade flapping itself is already helping to pull the tower back against the wind.

Even with this approach, where both the tower motion and blade motion are considered simultaneously, I still observe the same behaviour. Reversing the FW action around 0° and 180° rotor azimuth works reliably.

If this behaviour is fundamentally related to a 3P effect, then reversing the FW action individually as each blade passes the tower axis would make physical sense. However, reversing the FW action for all three blades simultaneously based on only one blade (or equivalently one half revolution of the rotor) does not appear to be consistent with any physical mechanism that I can identify.

I would appreciate any thoughts or suggestions regarding what physical effect I may be overlooking.

Kind regards,

[Jason.Jonkman]

Dear @Abhinay.Goga,

I agree with you regarding the Azimuth output.

However, I’m not following what you are describing enough to comment on your damping approach.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

Thank you for the response.

Let me take a step back and ask about the physical relationship first.

From a few basic simulations, I observed that changing the hydraulic flywheel fluid motion (which introduces Coriolis forces inside the blades) modifies the blades’ out-of-plane response, and this in turn changes the tower fore-aft vibration response.

Initially, I focused on directly reducing the tower motion. However, I started wondering whether it would make more sense to act on the cause instead. In other words, if I regulate the fluid movement to intentionally modify each blade’s out-of-plane motion, the combined blade loads should produce a different resultant hub force (or moment), which would then influence the tower fore-aft motion.

Looking at the BeamDyn outputs, the blade-root force in the local X-direction (RootFxr) represents the out-of-plane force for each blade. The aggregated low-speed shaft thrust force (LSSGagFxa) appears to resemble the tower fore-aft response quite well in my simulations.

On the other hand, the blade-root out-of-plane bending moment (RootMyr) and the aggregated low-speed shaft moment (LSSGagMYa) do not seem to correlate with the tower response nearly as well.

Would you therefore expect the blade-root out-of-plane forces to be a better quantity for estimating each blade’s contribution to the resultant hub loading than the blade-root bending moments? Or is there another quantity that you would recommend using for this purpose?

Thank you again for your insights.

Best regards,

[Jason.Jonkman]

Dear @Abhinay.Goga,

I’m not surprised that ElastoDyn output LSSGagMya is not well aligned with TTDspFA because LSSGagMya is expressed in the azimuth coordinate system, which spins with the rotor. ElastoDyn output LSSGagMys in the shaft coordinate system is likely more aligned with TTDspFA when the nacelle is not yawed.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

My FW-based fore–aft tower damping works well during certain time intervals for one fluid-flow polarity, but excites the tower during other intervals. [Orange:withoutFW and Blue :withFW]

I repeated the simulations with the fluid-flow polarity reversed. The previously damped intervals then became excited, while some previously excited intervals became damped.

This suggests that the main difficulty is a phase problem rather than FW control logic.

My current interpretation is that the tower and the individual blades can have different out-of-plane motion phases. While the tower is deflecting in the fore–aft direction, a blade may bend either in the same direction as the tower or in the opposite direction.

Using available OpenFAST outputs, I attempted to estimate the relative phase of each blade as follows:

if towerTopDeflectionError*tipDeflection > 0
    bladeDeflectionSync = 1;
else
    bladeDeflectionSync = -1;
end

if nacelleSpeedFA*tipDeflectionDerivative > 0
    bladeSpeedSync = 1;
else
    bladeSpeedSync = -1;
end

bladeAction = bladeDeflectionSync*bladeSpeedSync;

Here:

• towerTopDeflectionError is the tower-top FA displacement (TTDspFA) relative to a moving mean. I considered 3 seconds average based on Tower frequency.

• tipDeflection is the BeamDyn blade-tip out-of-plane translational deflection (B1TipTDxr).

• tipDeflectionDerivative is calculated numerically from the tip-deflection signal.

• nacelleSpeedFA is the nacelle FA translational velocity (Ncl_sp_FA).

I then defined a collective condition:

if blade1Action == 1 && blade2Action == 1 && blade3Action == 1
    fluidPhase = 1;
else
    fluidPhase = 0;
end

The FW was therefore activated only when all three blades were classified as being in the same synchronized state. However, the simulation still produced similar excitation intervals. In addition, this gating made the controller discontinuous and reduced the available FW action.

I am uncertain whether comparing blade-tip deflection and velocity directly with tower displacement and nacelle velocity is sufficient to determine the correct sign of the phase dependency between blade and tower.

I would greatly appreciate any insights on the following questions:

1. Is there a more appropriate way in OpenFAST to determine the instantaneous phase relationship between individual blade out-of-plane motion and tower fore–aft motion?

2. Could the changing relation between fluid-flow direction and tower response be caused by the rotating coordinate transformation, rather than by the blade–tower motion phase alone?

Thank you again for your time and feedback.

Best regards,

[Jason.Jonkman]

Dear @Abhinay.Goga,

The difficulty I see is that the phasing between the blade and tower motion depends on which full-system mode is being excited, and in reality, you likely have a combination of full-system modes being excited.

FYI: You can see how the blade and tower modes are coupled within each full-system mode by reviewing the corresponding eigenvector, e.g., using visualization functionality available through the Automated Campbell Diagram Code (ADCD): GitHub - OpenFAST/acdc: ACDC: Automated Campbell Diagram Code · GitHub.

Best regards,

[Abhinay.Goga]

Dear Jason and all,

I have now become familiar with the ACDC tool.

I am currently investigating an idling extreme-wind scenario based on DLC 6.1/EWM. The blades are pitched to 90°, the yaw error is zero, and GenDOF is enabled, so the rotor is not mechanically locked. Consequently, the rotor freewheels slowly. In my time-domain simulations, its speed oscillates around zero, with negative rotation direction.

Because there is no unique wind-speed–rotor-speed relationship under these idling conditions, I generated ACDC operating points for rotor speeds between approximately 0.1 and 0.5 rpm and wind speeds from 40 to 70 m/s.

For the first three identified modes, I observe the following general behaviour:

• One mode remains close to 0.326 Hz across nearly all wind speeds and rotor speeds.
• A second nearby mode increases from approximately 0.326 to 0.336 Hz with increasing wind speed.
• A third mode has a natural frequency near 0.51 Hz, but its damping ratio increases strongly toward unity as wind speed increases. Correspondingly, its damped frequency decreases considerably.

I also generated another set of operating points covering wind speeds from 35 to 75 m/s and rotor speeds from 0.01 to 0.8 rpm, while keeping the same configuration (GenDOF = True, pitch angle = 90°, neutral yaw).

To obtain valid linearizations with ACDC under these parked/idling conditions, I had to temporarily modify the settings by setting the cut-in wind speed to 35 m/s and the cut-out wind speed to 75 m/s. Otherwise, ACDC did not generate the operating points. I assume this is simply a workaround to keep the turbine model active for the linearization, rather than a physically meaningful operating configuration. Is this the correct way to handle such parked/idling cases, or is there a better approach?

The mode-shape visualization appears to show coupled tower and blade motion. However, I am still uncertain about identifying the first two closely spaced modes reliably as tower fore–aft, tower side–side, or mixed tower–blade modes. I understand that the numerical labels Mode 0, Mode 1, and Mode 2 merely indicate the ordering returned by the eigenanalysis and may not preserve one physical mode.

Could you please comment on the following points?

1. Near-zero rotor speed: Is ACDC/OpenFAST linearization reliable at rotor speeds as low as 0.01–0.8 rpm, or should the rotor be locked for a representative parked/idling eigenanalysis?
2. Mode identification: Is the recommended approach to identify the modes by examining their eigenvector amplitudes and phases rather than by relying on their frequency ordering?
3. Time-domain interpretation: In my nonlinear simulation, I estimate the frequencies of TTDspFA and TTDspSS, using moving average amplitudes cycles. Can these frequencies be used to identify which ACDC full-system mode is active?

I would appreciate your feedback on whether this interpretation of the ACDC results and mode-shape visualization is correct.

Best regards,