Stability analysis of operating wind turbines (NREL 5MW Onshore)

Dear Dr. Jason,

We are currently working on a research project on Wind Turbine stability analysis (CENER/Public University of Navarre) in which we want to apply Floquet Theory to analyze the stability of land based wind turbines using OpenFAST. Our approach requires to obtain linearized state-space (x,xdot,u,A,B,C,D) at several azimuth angles for a full revolution and, based on it, do the required numerical integration.

In the conventional approach, as discussed several times in this forum, the system is linearized around the rotor at several azimuth angles, linearizing around an operating point which needs to be periodically stable (when the transient is gone). Then, Multi Blade Coordinate change is performed to switch the rotating frame DOFs to the non-rotating frame, using MATLAB MBC3 code, for instance. Finally, averaging is performed azimuthally, and the eigensolution of the azimutally-averaged matrix computed. Real part of the eigenvalues will give information about damping ratios and the imaginary part of the eigenvalues of the damped frequencies. To identify which DOF affects most each machine mode, eigenvectors can be used.

On this context, we are currently working with the NREL Onshore 5MW model. We want to study the stability of the operating wind turbine at 13 mps incoming steady wind, rotating at rated velocity (12.1 rpm) (region III). Now we are studying in detail how the eigenvalues change around 1 rotor revolution. We compare our results of damped frequencies and damping ratios with results given in figure 4 shown in the paper:
Multi-blade coordinate transformation and its application to wind turbine analysis (G.Bir 2008).

In our model, DOFs are enabled for blades (3 DOF for each blade), tower (4 DOF), drivetrain and yaw.

We have chosen NLinTimes to be 360 (1º each linearization, being linearization number 1 at 0 deg of azimuth) in order to capture the dynamics of the blades passing in front of the tower. The time step was set to 0,0046 sec, being a multiple of (60 s/min)/(NLinTimes*RotSpeed) to ensure that we get a linearization at each azimuth degree.

We have applied 3 different approaches to obtain the OpenFAST linearization data:

A) Using CalcSteady=True, TrimCase=3 , VSControl=1 VS_RtTq=43094 Nm (rated), TimGenOn=0 and VS_RtGnSp = VS_Rgn2K = VS_SlPc = very small number.
This method has been discussed several times in the forum and is explained in FAST User Manual, so trimming the pitch in region 3 of the machine lead to convergence (0.0001 = TrimTol & 0.0001 = TrimGain) for pitch=6.01 deg (u=13 mps). After this, 360 linearizations are computed by OpenFAST. In the following figure, the real/imaginary part of the eigenvalues of A state-matrix are shown, after applying MBC3 transformation:

The variation of pitch and rotor speed over the linealizations is the following:

After applying the azimuth averaging of the MBC3-transformed A matrices, the damped frequencies given by the eigensolution are, in Hz:
6.1045 3.9450 2.9079 2.9562 2.0990 1.9901 1.7037 1.6940 1.2954 0.8880 0.7444 0.5659 0.3426 0.3359 0.3211

As GenDOF is enabled, a mode with pure real eigenvalues 0.1185 + 0.0000i and -0.2049 + 0.0000i appear. Does the real positive part of this mode reveal any real instability? how is so? Or it is it simply an nonphysical damping ratio?

There are some results that we can not totally understand:
• A 3P influence of the blades passing in front of the tower is observed, however, we expected to observe perfectly 3P periodical functions. As it can be seen, when 1st and 3rd blade pass in front of the tower, similar eigenvalue variation is observed, whereas different peaks are seen with the second blade (lin. Number around 60).
• The results of the damped frequencies do not match very well with the ones reported in the paper mentioned above (Fig.3, G.Bir 2008 ). i.e., in the figures of the article a damped frequency of 3.94 Hz is not shown, or 2 frequencies close to each other at around 1.7Hz.

  So, we have tested another linerization setting:

B) Using CalcSteady=False, , GenDOF=False, VSControl=0, GenModel=1, TimGenOn=9999.9 BldPitch=6.01 deg (for each blade) and RotSpeed=12.1 rpm
In this case we run a non-linear OpenFAST simulation with fixed pitch equal to the one found with CalcSteady functionality in A) and the rotor speed is set to the rated value. Visually a time where the transient is gone is identified, (550sec is chosen) and with the aid of a MATLAB script 360 LinTimes corresponding to Azimuth values of 0º,1º,2º…359º are calculated (and written in a OpenFast input file), choosing the closest values to these integers found in the non-linear time series.

Following an analogous approach, the variation of eigenvalues found now is:

The variation of pitch and rotor speed through the linealizations is the following:

The damped frequencies given by the eigensolution are, in Hz:
6.1040 3.7068 2.9079 2.9560 2.0988 1.9850 1.7037 1.2964 0.8889 0.7447 0.5654 0.3430 0.6303 0.3362 0.3141

Some comments about the results:

• With this approach, once again the expected periodic 3P variation of the eigenvalues is not seen, as it is not completely periodic. See how, around linearization number 300 there are no aggressive peaks on the real part of some eigenvalues.
• This time, the damped frequencies obtained from the eigensolution match well with the results shown in the paper by G.Bir.

Finally, a last linearization setting in OpenFAST was tested:

C) Using CalcSteady=False, , GenDOF=True, VSControl=1, TimeGenOn=0 s , VS_RtTq=43094 Nm, BldPitch=6.01 deg (for each blade) and RotSpeed=12.1 rpm
The procedure was exactly as in the previous approach, this time enabling the GenDOF and limiting the torque to the rated value of region 3 (as in A) approach )
The variation of eigenvalues obtained is:

The variation of pitch and rotor speed over the linealizations is the following:

The damped frequencies given by the eigensolution are, in Hz:

6.1045 3.9453 2.9080 2.9563 2.0987 1.9899 1.7033 1.6943 1.2975 0.8895 0.7447 0.5694 0.3428 0.3359 0.3210

Some comments about the results:
• This time, the variation of the eigenvalues when passing in front of the tower seems to be closer for each of the 3 perturbations (linearizations 60,180 and 360).
• The damped frequencies are very close to the ones with the A approach, but different from those obtained from the B approach, that matched the results of the NREL 5MW Onshore Turbine on the mentioned paper.

To sum up, we would be very pleased in you gave us your vision about the following issues:
1. Which of the 3 approaches (A, B or C) seem better for you to study the stability of a operating wind turbine by eigenanalysis?. Why, if we are linearizing a system with equal pitch and very similar rotor speed (subtle variations about 12.1 rpm in three approaches) we are obtaining different damped frequencies by the eigensolution?. Why approaches A and C do not match with the damped frequencies presented by G.Bir in the above mentioned paper?

2. When GenDOF=True (A and C), the variation of rotor speed versus linearization number is step-like (Rotor speed resolution in lin files is only 4 decimals) and lower in amplitude that when GenDOF is disabled (B approach). Why does this happen, and how do you thing it can affect the eigensolution?

3. Given that the conditions used for C) are extracted form A), and given the negligible variations in rotor sped an pitch in both cases, the linearizations are dynamically equivalent (same constraints and external conditions in the degrees of freedom of the system). So, how is that the eigenvalues at the position in which the blades go over the tower are not coincident?. Are we missing something?

4. What is causing the eigenvalues to be so sensible when blade is at the vicinity of the tower (strong peaks in the real part of eigenvalues), difficulting the obtention of a periodic 3P pattern with the rotor in steady-state and 3 equal blades. Is there a way in which we can obtain the expected perfect 3P periodic behavior using the different linearization schemes A) B) and C).

   Finally a suggestion for a future feature. It takes a lot of simulation time to reach periodic equilibrium before linearization can be done properly. It would be nice to be able to save a given state along a simulation, and later to be able to load it as initial state for other simulation. This will certainly save us and possibly others a lot of time.

Here I add some piece of information that could be useful for debbuging our issues:

• OpenFAST v 2.5 is used in the simulations, with ElastoDyn and AeroDyn v15 modules enabled.
• Model files were taken from Github repository (18-Feb 2021), and the input files, as well as the results of the nonlinear time series and linealizations (.lin) are available in the link (too large file to attach): [url][/url]
ServoDyn and ElastoDyn files were changed accordingly with the explanations above.
• High precision of "ES20.12E3"= OutFmt  is used in the simulations, in order to avoid numerical problems with the eigensolution.

We really appreaciate your help, time, opinion and guidance. Please, if you need more information or details about our procedure or models, do not hesitate to ask for it.

Thank you very much in advance for your time,

Best regards,

Javier Ros, Alvaro Olcoz and Aitor Plaza
Applied and Computational Mechanical Engineering research group
Public University of Navarre (Spain)

Dear Javier, Alvaro, and Aitor,

Overall, it sounds like you understand the linearization process and are following the correct approach. The big difference between method B and methods A and C is that the generator DOF has been disabled in B. This will eliminate the rigid-body mode and have some influence on the natural (and damped) frequencies of the drivetrain torsion, tower side-side and blade edgewise-bending modes. I was not involved in the paper you reference, so, I’m not sure whether the generator DOF was enabled or not, but seeing that your results from B match the paper well, I would guess the generator DOF was diasabled in the paper as well.

Regarding the rigid-body mode, the damping of a rigid-body mode is undefined, so, I’m not sure the real part of the eigenvalue has any physical relevance.

To answer your numbered questions:

  1. I would generally recommend approach A when linearizing with wind above rated. Again, B produces different results from A and B because of the absence of the rigid-body mode of the generator.
  2. I would guess the stepping is related to the precision of the output. It sounds like you have already chosen to output a high precision (as recommend for linearization analyses), but you may want to check the precision of the value you are plotting specifically.
  3. I’m not sure, as I’ve not really plotted the eigenvalues as a function of azimuth angle before. Considering that you are comparing “spikes”, I would guess the solution is extra sensitive to the exact conditions (perhaps a small change in azimuth has a large effect).
  4. I would guess this is related to the influence of the tower on the local wind speeds, again, with perhaps a small change in azimuth having a large effect.

Regarding your feature request, I agree. We are actually working on developing this capability in a collaborative project with Envision Energy. The capability under development will allow both initialization of a model based on a precomputed state and the direct calculation of a steady-state solution at model initialization within OpenFAST.

I hope that helps.

Best regards,