Drive train dynamics equations

[alvart-98]

Hello,

I am working in the OpenFast with the NREL 5 MW, which has a semi-submersible platform. I have implemented a state-space system for the rotor and generator rotational speeds basing on the bibliography: [Ir,0;0,Ig][theta_rdot2;theta_gdot2]+[Kl,-Kl/Ng;-Kl/Ng,Kl/Ng^2][theta_r;theta_g]+[Bl,-Bl/Ng;-Bl/Ng,Bl/Ng^2][theta_rdot;theta_gdot]=[Tr;-Tg], where: Ir=3.88310^7kgm^2 is de rotor inertia, Ig=534.116kgm^2 is the generator inertia, theta_r is the angular position of the rotor, theta_g the angular position of the generator (dot and dot2 are the first and second time derivatives of the angular speeds), Kl=867637000 Nm/rad is the spring constant, Bl=6215000 Nm*s/rad is the damper constant and Ng=97 is the gearbox ratio.

The issue consists in the implementation of this model in OpenFast, since the calculation of the rotor inertia is made internally. I have simulated two models: the first one, with the rotor inertia Ir=3.88310^7kgm^2 (I have checked on the forum of NREL, the result is almost the same that Mr. Jason Jonkman provides, so I understand that I can assume this value), and the second one, ONLY with the inertia of the hub. The first model is not even close to the results that OpenFast calculates, however, the second one, gives almost exactly the same results than OpenFast (only small error is found), which does not make sense at all.

Maybe, my set-up configuration is wrong, but I do not know in which file can I change the error. Could you provide me some kind of help to adjust the inertia parameters of the rotor so as to solve the issue?

Thank you very much for your attention.

Mr. Ortiz

Graphs:

1) First case: Ir=3.88310^7 kgm^2

2) Second case: Ir=115926 kg*m^2

OpenFAST Version OpenFAST-v2.0.0

---

OpenFAST

Copyright (C) National Renewable Energy Laboratory Copyright (C) Envision Energy USA LTD

This program is licensed under Apache License Version 2.0 and comes with ABSOLUTELY NO WARRANTY. See the "LICENSE" file distributed with this software for details.

---

OpenFAST-v2.0.0 Compile Info:

• Architecture: 64 bit
• Precision: double
• Date: Nov 27 2018
• Time: 17:19:38 Execution Info:
• Date: 06/04/2021
• Time: 10:00:00-11:00:00

System Information (please complete the following information):

• OS: Windows 10
• Compiler: Intel Fortran oneAPI HPC Toolkit

Additional context Matlab & Simulink: 2018b & 2020b

[ebranlard]

Hi @alvart-98

Can you send a link to the bibliography you mention that shows the equations? I'm wondering if something is off in the conversions from low speed shaft to high speed shaft.

As a small test, you can also try a single degree of freedom, neglecting shaft torsion and transferring everything to the low speed shaft (LSS):

J_LSS theta_r_dot2 = Tr-Tg_LSS = RtAeroMxh - GenTq*1000*Ng

with J_LSS = Ir + Ig*Ng^2 ~= 4.38e7 (I find that I usually have to increase the rotor inertia a bit, J_LSS=4.45e7). Maybe that helps?

[alvart-98]

Hello @ebranlard

The book that I have used is the following one: Wind Energy Systems, Control Engineering Designs. By Mario García-Sanz and Constantine H. Houpis.

I have simplified the model shown in Chapter 12.3, just by assuming that the angular position in the generator shaft and the one in the low speed input of the gearbox are related by the ratio Ng=97, so that 97*theta_L=theta_G.

If I consider the simplification you are mentioning, which is that 97*theta_R=theta_G, I obtain such equation. I have also tried that procedure, but I still have the issue. In my opinion, the problem is totally related to the rotor inertia. I do not know if its calculation is wrongly implemented in OpenFast (what would be a huge mistake, so I do not consider this option) or my set-up is wrong in some aspect.

Thank you for answering!

[ebranlard]

Hi @alvart-98

Unfortunately, I don't have access to this reference . I believe the equation that is solved internally by OpenFAST is the following :

[Jr + Jg N^2, -Jg N^2 ][theta_r]_dd + [0, 0][theta_r] + [0, 0][theta_r]_d =  Q_r - Q_g N
[  -Jg N^2  ,   Jg N^2][ nu    ]_dd + [0, K][ nu    ] + [0, B][ nu    ]_d =  Q_g N

where nu = theta_r - theta_g/N is the torsion angle between the two angles. This is probably the most stable system of equations since it keeps track of the difference between the angle and keeps everything on the low speed shaft side. The equations can be rewritten (still keeping everything on the low speed shaft) as:

[Jr , 0    ][theta_r    ]_dd + [ K,-K][theta_r    ] + [ B,-B][theta_r    ]_d =  Q_r 
[0,  Jg N^2][theta_g_LSS]_dd + [-K, K][theta_g_LSS] + [-B, B][theta_g_LSS]_d =  Q_g N

Or, as you suggested, with each degree of freedom on each side of the gearbox, as follows:

[Jr , 0 ][theta_r]_dd + [ K  ,-K/N  ][theta_r] + [ B,  -B/N  ][theta_r]_d =  Q_r 
[0,  Jg ][theta_g]_dd + [-K/N, K/N^2][theta_g] + [-B/N, B/N^2][theta_g]_d =  Q_g 

I've made a small python script here to integrate these equations based on an OpenFAST input file and output file. The first set of equations match OpenFAST exactly, but the two others show some differences, likely due to numerical noise. I use Q_r = RtAeroMxh and Q_g = GenTq*1000, J_r= 3.85e7, J_g=534.116, N=97. I hope this will help troubleshoot your model, let us know what you find out.

[alvart-98]

Hello,

I have just found the problem. The state-space system you are describing is the most suitable, furthermore, the second and third ones also work almost perfectly. The issue is that we were using the variable 'RotTorq' instead of 'RtAeroMxh', basing on the wrong idea that both were the same.

At this point, I would like to ask why they are not the same. I mean, both should represent the torque that the rotor is receiving from its own dynamics, shouldn't they?

In addition, how would you introduce the term due to the losses of the gearbox in the state-space system mentioned before?

Thank you very much for your attention.

[alvart-98]

Hi,

Could you explain me the difference between the variables 'RotTorq' and 'RtAeroMxh'?

Also, how would you introduce the term due to the losses of the gearbox in the state-space system mentioned before?

[jjonkman]

Closed and replaced with https://github.com/OpenFAST/openfast/issues/772.

[JalalHeidari]

Dear @ebranlard,

I wrote a code to calculate rotor speed with this equation: J_LSS theta_r_dot2 = Tr-Tg_LSS = RtAeroMxh - GenTq1000Ng

However, the error accumulates over time. I tried to decrease time step, use smoothing approaches or apply other higher order time integration methods. The result was the same. Below you can see the code, and attached is OpenFast out files. Could you please let me know what is the problem here? Thanks for your help

This is the code:

import pandas as pd import numpy as np import matplotlib.pyplot as plt

data = pd.read_csv('Data.csv')

aerodynamictorque = data['RtAeroMxh[N-m]'].values # Aerodynamic torque in N-m > Tr generatortorque = data['GenTq[kN-m]'].values 1000 # Generator torque in N-m > Tg actual_rotorspeed = data['RotSpeed[rpm]'].values 2 * np.pi / 60 # Rotor speed in rad/s

N = 97 # Gearbox ratio - Jg = 534.116 # Generator inertia kgm^2 Jr = 38677056 # Rotor inertia kgm^2 time_step = 0.01 # Time step in seconds

calculated_rotor_speed = np.zeros_like(actual_rotor_speed) calculated_rotor_speed[0] = actual_rotor_speed[0]

for k in range(1, len(actual_rotor_speed)): acceleration = (aerodynamic_torque[k-1] - generator_torque[k-1] N ) / (Jr + Jg N*2) calculated_rotor_speed[k] = calculated_rotor_speed[k-1] + acceleration time_step

plt.figure(figsize=(10, 5)) plt.plot(actual_rotor_speed, label='Actual Rotor Speed', linestyle='-', color='blue') plt.plot(calculated_rotor_speed, label='Theoretical Rotor Speed', linestyle='-', color='red') plt.xlabel('Time Steps') plt.ylabel('Rotor Speed (rad/s)') plt.title('Comparison of Actual and Calculated Rotor Speed') plt.legend() plt.show()

Data.zip

[jjonkman]

Dear @JalalHeidari,

Just a couple comments:

• Your equation of motion for rotor acceleration assumes that the wind turbine is rigid other than the variable-speed degree of freedom, but I see from your OpenFAST output that the model you are comparing to and extracting torque output from is not rigid.
• Your integrator is very simple (Euler) that is only first-order accurate.

Best regards,

[JalalHeidari]

Dear @jjonkman

Thank you for your reply and help.

1- Regarding using other integrators, I used RK4 and ABM4, but the error persists.

2- How I can make OpenFast to use a rigid structure and verify the mentioned equation? J_LSS theta_r_dot2 = Tr-Tg_LSS = RtAeroMxh - GenTq1000Ng

I think aerodynamic_torque average in steady state should be 43.090971000 =4,179,730 , but RtAeroMxh average is not. However, RotTorq average is.

As a result, If I bring everything to HSS and use "RotTorq", apply this equation: dwg/dt = (RotTorq/N - GenTq*1000)/(Jg)

Also If I bring everything to LSS, dwr/dt = (RtAeroMxh - RotTorq)/ (Jr)

Are these equations correct?

I also not quiet sure whether aerodynamic_torque is a function of RtAeroMxh, RtAeroMyh, and RtAeroMzh? Why we only considering RtAeroMxh in the equations?

You mentioned in #772 that RotTorq is Knu + Bnu_d I wanted to know how can I get theta_r and theta_g in OpenFast and calculate nu and nu_d. Q_GeAz is theta_g and Q_DrTr is theta_r?

If i am correct about that, then I don't understand why Q_DrTr is not changing between 0 and 2pi (similar to Q_GeAz) in my data set.

Best regards

[jjonkman]

Dear @JalalHeidari,

I haven't reviewed all of your math, but just a few comments:

• The following topics on our forum are likely of interest to you: https://forums.nrel.gov/t/resistant-moment-of-the-rotor-and-of-the-electric-generator/408 and https://forums.nrel.gov/t/drivetrain-dynamic-model/1864.
• To model a rigid turbine (other than the variable-speed generator degree of freedom) in OpenFAST, disable all degrees of freedom in ElastoDyn except GenDOF (and make sure BeamDyn is disabled as well).
• For a torque-only model of the drivetrain, the transverse aerodynamic moments (RtAeroMyh, RtAeroMzh) wont influence the drivetrain acceleration/deceleration or torque-related loads.
• I agree that Q_GeAz is theta_g, but Q_DrTr is nu, such that theta_r = Q_GeAz + Q_DrTr.

Best regards,

[JalalHeidari]

Dear @jjonkman

Thank you for your help and links. For a rigid structure, the problem is solved, even for a stochastic wind profile:

Now I want to explore this equation:

[Jr + Jg N^2, -Jg N^2 ][theta_r]_dd + [0, 0][theta_r] + [0, 0][theta_r]_d = Q_r - Q_g N [ -Jg N^2 , Jg N^2][ nu ]_dd + [0, K][ nu ] + [0, B][ nu ]_d = Q_g N

I have these conditions: ---------------------- DEGREES OF FREEDOM -------------------------------------- True FlapDOF1 - First flapwise blade mode DOF (flag) True FlapDOF2 - Second flapwise blade mode DOF (flag) True EdgeDOF - First edgewise blade mode DOF (flag) False TeetDOF - Rotor-teeter DOF (flag) [unused for 3 blades] True DrTrDOF - Drivetrain rotational-flexibility DOF (flag) True GenDOF - Generator DOF (flag) True YawDOF - Yaw DOF (flag) True TwFADOF1 - First fore-aft tower bending-mode DOF (flag) True TwFADOF2 - Second fore-aft tower bending-mode DOF (flag) True TwSSDOF1 - First side-to-side tower bending-mode DOF (flag) True TwSSDOF2 - Second side-to-side tower bending-mode DOF (flag) False PtfmSgDOF - Platform horizontal surge translation DOF (flag) False PtfmSwDOF - Platform horizontal sway translation DOF (flag) False PtfmHvDOF - Platform vertical heave translation DOF (flag) False PtfmRDOF - Platform roll tilt rotation DOF (flag) False PtfmPDOF - Platform pitch tilt rotation DOF (flag) False PtfmYDOF - Platform yaw rotation DOF (flag)

A fixed 20 m/s wind speed is applied. In steady state dwg/dt = (RotTorq/N - GenTq*1000)/(Jg) should be 0. I wrote a simple code to check the mean values in steady state, this is the result:

Average RotTorq between 300 and 600 seconds: 4,180,039.50 N-m Average (GenTq 1000 N) between 300 and 600 seconds: 4,179,730,000 N-m

This slight difference accumulates over time again. I tried to change other degrees of freedom. I disabled all except Generator DOF and Drivetrain rotational-flexibility DOF, and this solved the problem.

From my understanding this equation:

[Jr + Jg N^2, -Jg N^2 ][theta_r]_dd + [0, 0][theta_r] + [0, 0][theta_r]_d = Q_r - Q_g N [ -Jg N^2 , Jg N^2][ nu ]_dd + [0, K][ nu ] + [0, B][ nu ]_d = Q_g N

is not valid when other degrees of freedom is enabled. There is a slight difference.... Can you please let me know what are the equations inside OpenFast when other degrees of freedom are also active? I couldn't find the equations in ElastoDyn theory manual...

[jjonkman]

Dear @JalalHeidari,

I agree that the equations of motion of ElastoDyn are much more complicated than as you stated when degrees of freedom other than generator rotation and drivetrain torsion are enabled. The full equations of motion cannot be stated in a post. Instead, they are documented in the unoffocial ElastoDyn theory manual: https://openfast.readthedocs.io/en/main/source/user/elastodyn/index.html.

Best regards,