Natural frequency and damping ratio calculation

Wind & Water Structural Analysis
1

Dear Jason,

Sorry to ask again about the natural frequency calculation, since I find there are quite many posts on this topic, but still I have some confusion.

I would like to calculate the natural frequencies and damping ratios for some modes of OC3-Hywind.

Before model linearization, I activated 6 platform DOFs and 4 tower DOFs, set rotor rotation speed and wind speed as 0, disabled incident wave or radiation damping, disabled all control inputs.
After linearization, since there is no rotating states, I used “eianalysis.m” from MBC package to produce natural frequencies and damping ratios. The result is listed below.

surge: 0.0080 Hz
sway: 0.0080 Hz
heave: 0.0324 Hz
roll: 0.0342 Hz
pitch: 0.0342 Hz
yaw: 0.1211 Hz
TwFA1: 0.4855 Hz
TwFA2: 2.6606 Hz
TwSS1: 0.4790 Hz
TwSS2: 2.3715 Hz

I found The report NREL/SR-500-45891 also listed the analysis result of OC3-Hywind as shown in this figure.

It seems the results for all the platform DOFs agree well, but the natural frequencies of tower DOFs are quite different, especially for 2nd order.

Is it because the absence of some other DOFs? or I should introduce aerodynamics?

Kind regards,
Yulin

2

Dear Yulin,

Your linearization process sounds reasonable.

I believe there is quite a bit of coupling between the tower-bending modes and the rotor blade fexibility for this model. This coupling will have a small impact on the 1st tower-bending frequencies and a large impact on the 2nd tower-bending natural frequenices. I suspect you’ll be able to match the results to NREL/SR-500-45891 much better if you enable the blade-bending DOFs in FAST. More information is available in the forum topic found here: Tower Eigenfrequencies of NREL 5MW Turbine.

Best regards,

3

Dear Jason,

Thank you for the quick and clear reply.

You are right that the blade flexibility does have big influence on the 2nd order tower bending mode. Besides, I found the drivetrain also has big influence on tower side-side bending modes.
I followed your advice by activating 9 blade DOFs and also the drivetrain DOF in addition to previous setup, and it now agrees well with the following result in report NREL/TP-5000-48191.

One other question. In the results of “eianalysis”, the order of states seems different with that defined in .lin file. For example, in the .lin file, the order of states is “surge, sway, heave, pitch, roll, yaw”, but in the natural frequency matrix “NaturalFrequencies”, the order might become “yaw, heave, sway, surge, roll, pitch”. I am wondering how to identify their correspondence.

Thanks again.

Best regards,
Yulin

4

Dear Yulin,

I’m glad you now get the answers you expected.

See the forum topic found here for information on how to interpret (and a tool to aid in interpreting) the results of an eigenanalysis: Eigenanalysis FAST.

Best regards,

5

Dear Jason,

Thank you very much. I can manage now~

Sincerely,
Yulin

6

When the user specifies the percent critical damping using “TwrFADmp” and “TwrSSDmp” in the tower.dat files, how is that damping applied to the structure? - i.e., Does FAST use Rayleigh damping?

7

Dear Wystan,

Good question. I’m surprised this one hasn’t been asked before.

The tower damping in FAST is specified as a fraction of critical damping for each mode of the isolated tower. The damping is implented as a stiffness-proportional modal damping, mathematically as follows. If the i’th generalized elastic force is given as:

Felastic_i = K_ij*q_j,

where K_ij is the (i,j) element of the generalized elastic stiffness matrix of the isolated tower and q_j is the j’th tower-bending degree-of-freedom (DOF), then the i’th generalized damping force is given by:

Fdamping_i = K_ijbeta_jqdot_j,

where qdot_j is the first time-derivative of q_j and beta_j is the stiffness-proportional constant for the j’th DOF. (Einstein notation is used here.) The stiffness-proportional constant is derived from the modal damping ratio of the j’th tower-bending mode (zeta_j) and the j’th natural frequency of the isolated tower (f_j) as follows:

beta_j = zeta_j/(pi*f_j),

where:

f_j = 1/(2*pi)*SQRT(K_jj/M_jj),

where M_ij is the (i,j) element of the generalized mass matrix of the isolated tower.

A few points of clarification for FAST:
*i and j both go from 1 to 2 for both the fore-aft and side-to-side bending modes.
*zeta_j represents TwrFADmp(j)/100 and TwrSSDmp(j)/100 for the fore-aft and side-to-side bending modes, respectively.
*The generalized elastic stiffness matrix of the tower depends on the specified distributed tower-bending stiffness (including any scaling or tuning factors) and the mode shapes. It does not depend on the gravititional destiffening or the platform or rotor-nacelle assembly (RNA) DOFs, so, it really represents the elastic stiffness only of the isolated tower.
*The generalized mass matrix of the tower depends on the specified distributed tower mass (including any scaling factors) and the mode shapes. It does not depend on the mass of the platform or RNA, so, it really represents the mass only of the isolated tower.
*Because the stiffness-proportional constant is calculated only for the isolated tower, the damping of the tower-bending modes when coupled with the platform and RNA may be different than specified (i.e., the damping ratio derived for the tower-bending modes through a FAST linearization analysis of the full wind turbine system may differ from what is specified for the isolated tower in the FAST input file).

I hope that helps.

Best regards,

1 Like
8

Jason,

Yes it does! Thank you.

– Wystan

9

Dear Jason
As I know,the Rayleigh damping have two proportional factor about the mass and stiffness.As you said, Fdamping_i = K_ijbeta_jqdot_j
I think the mass is neglected.If not,
Fdamping_i=M_ijalph_jqdot_j+K_ijbeta_jqdot_j
I am not sure about the question.Hope your reply.
BEST REGARD
Ruiliang.Wang

10

Dear Railung,

I agree that general Rayleigh damping includes both stiffness- and mass-proportional terms. For a given alpha and beta, constant for all modes, the end result is that the alpha (mass-proportional) term has more effect at low frequencies and the beta (stiffness-proportional) term has more effect at high frequencies.

However, the tower damping in FAST is not identical to Rayleigh damping. The tower damping in FAST is only based on the stiffness-proportional term and is specified independently for each mode.

Best regards,

11

Hello,

I have another question about this topic: when I do simulations the damping that I see in the simulation output is about half of the damping that I specify in the FAST input.

This possibly relates to Jason’s comment in his post above:

.

To minimize the influence of the other modes I disabled all DOF’s except the tower FA mode. Surely the output damping in that case must correspond to the input damping? Well, I still had the mass of nacelle and rotor not set to zero. Is that fact alone responsible for this factor 2 difference between input and output?

The way I determined the damping in the output was: Turbine running in vacuum (CompAero=false), all DOFs turned off except the first tower mode FA, turbine is parked (so no dynamic damping), initial tower displacement 1 m at the top, then releasing it and inspecting how the amplitude is decreasing.

I did 5 simulations like this, each with a different input damping ratio zeta. I determined the damping from the output by calculating the damping decrement of the decaying tower oscillation.

The simulation results are in the table below (Sorry for the unclear format: it is a table with 1 row header column, 5 simulation result columns and 1 comment column. I could not figure out how to present it more neatly, any spaces and tabs I put disappear).

Simulation 1 Simulation 2 Simulation 3 Simulation 4 Simulation 5 comment logd 0.01 0.02 0.03 0.04 0.05 Suppose we want these logd’s zeta 0.0016 0.0032 0.0048 0.0064 0.0080 =logd/(2*pi) TwrFADmp1 0.159 0.318 0.477 0.637 0.796 =zeta*100 (FAST input value) x3 9.89E-01 9.79E-01 9.69E-01 9.59E-01 9.49E-01 3rd peak in graph below x12 9.29E-01 8.64E-01 8.03E-01 7.46E-01 6.94E-01 15th peak in graph below logd 0.0053 0.0105 0.0156 0.0209 0.0261 =1/12*ln(x3/x12) zeta 0.0008 0.0017 0.0025 0.0033 0.0041 =logd/(2*pi)
Note the bold numbers, the first bold row is the input logd (which I convert to zeta and % before putting into FAST), the last bold row is the damping determined from the decaying oscillation. There is a factor 2 difference.

Can that be right?

12

Dear Gerrit,

Yes, it is the mass and inertia of the RNA that is resulting in a different level of damping in the response than what you input as a damping ratio. The tower damping ratios specified in FAST apply only to the isolated tower, in the absence of the platform or RNA. Effectively, the RNA changes the generalized mass of the tower-bending mode, dropping the natural frequency and requiring a higher level of damping to achieve an equivalent damping ratio when compared to an isolated tower without the RNA. The RNA also impacts the gravitational destiffening of the tower; this term tends to be small, but has more importance the higher the mass of the RNA.

In your case it sounds like doubling the damping ratio for the isolated tower will give the damping ratio you expect for the combined tower + RNA.

Best regards,

13

Dear Jason,

In my model the frequency of the tower+RNA was indeed a factor two lower than the frequency of the isolated tower.

Thanks

14

Dear @Jason.Jonkman ,

Is there any reason for choosing stiffness proportional modal damping rather than mass proportional modal damping ?

Best

Riad

15

Dear @Jason.Jonkman ,

Also i’m a bit confused between modal damping where we obtain the damping matrix in the eigenvector basis and we have to transform it in the physical coordinates and the damping matrix where we have not to do this transformation.

Please could you help in this point.

16

Dear @Riad.Elhamoud,

In ElastoDyn, the damping ratio is specified separately for every mode rather than as a full mass or stiffness matrix scaling, so, there is no reason to use Rayleigh damping with a mass-proportional term.

I’m not sure I understand your second question regarding different damping approaches.

Best regards,

17

Dear @Jason.Jonkman,

What i meant that in the eigenvectors coordinates, the damping matrix is diagonal.
Each term is obtained using the formula : c(i,i)=2damping_ratiocircular_natural_frequency.
Once, we obtain the diagonal damping in the eigenvectors coordinates, we should transform it to the physical coordinates. Is it correct ?

Also, i have a question concerning the damping matrix which account for radiation damping when using strip theory. This matrix is not calculated. It is obtained from experiment right ? If not, could you explain how to calculate it ?

Best,

Riad

18

Dear @Riad.Elhamoud,

Whether to transform equations between modal and physical coordinates will depend on the use case. I’m not sure what use case you are referring to.

The strip-theory solution in the HydroDyn module has viscous drag, but no radiation damping. Radiation damping is only considered in the potential-flow solution.

Best regards,

1 Like
19

Dear all,

I have developed a simple model of WP1.5MW turbine in MATLAB. In which tower is modeled as a single spring-damper-mass. Simulated and compared results for two different wind profiles i) Step and ii) NTM.

Simple model:

  • M_eq * X_ddot + C_eqTwr * X_dot + K_eq*X = Thrust
  • I_eq * Omega_dot = T_aero – T_gen

M_eq = M_eqTwr + M_twrTop

I see the results matching well between both OpenFAST model and Inhouse (MATLAB model) for the variables “generator speed” and “pitch angle”, but the “tower top displacement (ttd)” from my MATLAB model are higher than those obtained from OpenFAST. The results are clearly showing that the damping in my model seem to be lesser than that considered in OpenFAST.

Tower top displacement corresponding to step wind profile suggests that stiffness and mass are same between the same (because frequency is same between the two models Fig1_Step_comp_omega_ttd).

Note1 : To make OpenFAST model similar to MATLAB model, In OpenFAST model too, I have kept only tower fore-aft 1st mode is activated, all blades are treated rigid. Gravity set to zero; Tower – uniform tower.

Note2: I have used same value for 'zeta" in my model and “TwrFADmp” in OpenFAST.
Note3: I am using “Elastodyn” for structural dynamics

Below is how I defined damping.

Tower mass = M_Twr ( Just the tower mass alone, NOT including mass at the top )

Tower equivalent mass ( M_eqTwr) = 33/140*M_Twr

Tower stiffness (K_Twr) = 3EI/L^3

Damping ratio = Zeta;

Critical damping ( C_eqTwr ) = 2sqrt(K_TwrM_eqTwr)

Damping ( C ) = 2 * Zeta * C_eqTwr

For the case of Fore_aft mode alone, what is the relation between the tower damping used in OpenFAST and damping defined above. Any other insights on why the tower top displacements are not matching between the two models are appreciated. Thanks.

Fig1_Step_comp_omega_ttd
Fig1_Step_comp_omega_ttd1228Ă—462 58.7 KB

Fig2_Step_Comparison
Fig3_NTM_comp_omega_pitch
Fig3_NTM_comp_omega_pitch1132Ă—415 57.9 KB

Fig4_NTM_TTD_Comparison
Fig4_NTM_TTD_Comparison988Ă—521 64.7 KB

Regards,
Kumara

20

Dear @KumaraRaja.Eedara,

I’m not fully sure how you derived M_eqTwr and K_Twr. For example, for K_Twr, what value of EI are you using (the bending stiffness of the WindPACT 1.5-MW wind turbine is a function of elevation along the tower)? Are where does the factor 33/140 come from?

Regardless, I think the biggest term your simplified in-house code is missing is the aerodynamic damping, which is generally large in the fore-aft direction.

You can quantify the effect of aerodynamic damping–at least with open-loop control–through the OpenFAST linearization process. But it sounds like you have implemented some sort of controller as well, which would impact the closed-loop aerodynamic damping. To estimate this, you’d have to include a linear representation of the controller in the linearized system.

Best regards,