OpenFAST: Linearization of NREL 5MW Onshore Turbine

[Hadi.Hammoud]

Dear Dr.Jonkman,
I’m writing in hope that you could help with an issue I’m encountering during the linearization of the NREL 5MW wind turbine. I will first outline the details of the procedure I followed for the linearization, then I will present the problems that I encountered.

-I used OpenFAST with the input file attached below to linearize the NREL 5MW Onshore turbine in Region 3 about an 18 m/s wind speed.
-In the ServoDyn input file, the following setting were used: PCMode=0; VSContrl=1; VS_RtTq=43093.55; VS_RtGnSp=VS_Rgn2K=VS_Rgn2K=999e-9.
-3 DOFs were activated: GenDOF, FlapDOF1 (first flapwise mode), and TwFADOF1 (first fore-aft mode); that is a total of 10 continuous states.
-Linearization was performed successfully at 36 azimuth steps using pitch trimming, it resulted in a steady-state pitch angle of 14.77 degrees for a nominal rotor speed of 12.1 rpm.
-The resulting linearization output files were then loaded into MATLAB by defining a FileNames cell array containing the names of the 36 files then running the GetMats_f8 script.
-The MBC transformation was performed afterwards, and the resulting LTI A,B,C,D matrices were azimuth averaged using the function mean(X,3), i.e. averaging over 3rd dimension.
-The generator azimuth state was then eliminated from the state-space model by removing the corresponding rows and columns.
-A Simulink model was then built with the state-space matrices and the response of the system was simulated for a 1 degree collective pitch angle disturbance introduced after 10s (with all other inputs set to zero).
-The response was compared to that of the non-linear OpenFAST model for the same operating conditions and disturbance, the response is attached below.
-For the linear Simulink model, the Falpwise bending moment response (attached below) was obtained by performing the inverse MBC transform on the 0,C, and S outputs of the model.

The encountered issues are the following:

-Instead of the conventional 9 inputs that the linearized OpenFAST model should expect, the output linearization files that I obtained stated that 66 inputs are expected, with input 10 through 66 all called “AD User property on blade i, node i”. The corresponding columns in the state-space matrices are however all zeros.
-After applying the MBC transform and azimuth-averaging, the columns of the A and C matrices corresponding to the generator azimuth state are not all zero as expected.
-Although the linear and non-linear models agree for the rotor speed and tower foreaft responses, they are completely different for the flapwise moment.

Attached below are all the relevant documents.
I’m sorry for the lengthy post, I wanted to make sure that everything is clear. I appreciate the help and am extremely thankful.

[Hadi.Hammoud]

For some reason, the linearization file could not be uploaded, so here is the linearization section of the OpenFAST input file:
True Linearize - Linearization analysis (flag)
True CalcSteady - Calculate a steady-state periodic operating point before linearization? [unused if Linearize=False] (flag)
3 TrimCase - Controller parameter to be trimmed {1:yaw; 2:torque; 3:pitch} [used only if CalcSteady=True] (-)
0.0001 TrimTol - Tolerance for the rotational speed convergence [used only if CalcSteady=True] (-)
0.0001 TrimGain - Proportional gain for the rotational speed error (>0) [used only if CalcSteady=True] (rad/(rad/s) for yaw or pitch; Nm/(rad/s) for torque)
0 Twr_Kdmp - Damping factor for the tower [used only if CalcSteady=True] (N/(m/s))
0 Bld_Kdmp - Damping factor for the blades [used only if CalcSteady=True] (N/(m/s))
36 NLinTimes - Number of times to linearize (-) [>=1] [unused if Linearize=False]
30, 60, 80 LinTimes - List of times at which to linearize (s) [1 to NLinTimes] [used only when Linearize=True and CalcSteady=False]
1 LinInputs - Inputs included in linearization (switch) {0=none; 1=standard; 2=all module inputs (debug)} [unused if Linearize=False]
1 LinOutputs - Outputs included in linearization (switch) {0=none; 1=from OutList(s); 2=all module outputs (debug)} [unused if Linearize=False]
False LinOutJac - Include full Jacobians in linearization output (for debug) (flag) [unused if Linearize=False; used only if LinInputs=LinOutputs=2]
False LinOutMod - Write module-level linearization output files in addition to output for full system? (flag) [unused if Linearize=False]

[Jason.Jonkman]

Dear Hadi,

The approach you describe all sounds correct.

Regarding your questions:

• The AD user property inputs (UserProp) are standard inputs that can be ignored if not needed. See a related issue on the OpenFAST github page: github.com/OpenFAST/openfast/issues/778.
• I would suggest plotting the azimuth variation of the MBC-transformed A matrix (etc.) to see the periodicity in the signal (before azimuth averaging). It could be that the average is not exactly zero due to numerical round off, but they are likely much small than other dominant terms in the MBC-transformed A matrix. If you are not understanding the linearization output, I would always suggest simplifying the model to debug, e.g., does the result make more sense when you only have one DOF enabled (GenDOF = True; all others DOFs False)?
• Perhaps you’ve applied the inverse MBC transformation incorrectly for the flapwise bending moment? Do other outputs make more sense, e.g., blade flapwise deflection (requires inverse MBC transform), rotor thrust, tower-top deflection, and tower-base moment (the latter three do not require the inverse MBC transform)?

Best regards,

[Hadi.Hammoud]

Dear Dr. Jonkman,

After looking at the periodicity of the MBC-transformed matrices, it seems that indeed, the azimuth-averaging did not result in exactly zero values for the terms in the columns corresponding to the generator azimuth state. However, I’m not able to tell whether or not these non-zero terms are have a significant effect on the results of the simulation, and whether they could be causing the problems that I’m facing. Could you please suggest some way to verify that? Below are the entries of the azimuth-averaged MBC-transformed A matrix if this can be helpful. Column 2 corresponds to the azimuth state.

Concerning the inverse MBC transformation, I simply combined the linear model outputs corresponding to the three flapwise bending moments for blade 1 using the following equqtion: q1 = q0 + qccos(theta) + qssin(theta) where q0, qc, qs are the outputs of the linearized model and theta was found by integrating the rotor speed.
I find myself in a difficult situation: I’m seeking to validate the results of my linearization by comparing them to the non-linear model, and, although the two models agree for the rotor speed and tower foreaft moment, the discrepancies in the flapwise moment are causing me a lot of trouble as I cannot seem to find a cause for that error. I appreciate any help from your part.

Best,
Hadi Hammoud

[Jason.Jonkman]

Dear Hadi,

The state-matrix A for a first-order system derived from states that have second time derivatives can be partitioned as:

A = [ 0, I; -M^-1 * K, -M^-1 * C ]

where
0 = matrix of zeros
I = identity matrix
M = mass matrix
K = stiffness matrix
C = damping matrix
^-1 = matrix inverse

So, the bottom-left quadrant of your matrix is -M^1 * K. The diagonal elements of this submatrix are likely the most important. And clearly your results show that the (2,2) element of this submatrix, which is associated with the generator azimuth state, is very close to zero relative to the other diagonal elements (1,1) and (3,3) through (5,5). So, I’m quite sure that the value is not exactly zero due to numerical round off. Of course, by plotting the (2,2) element of this submatrix as a function of azimuth angle, hopefully this will become more clear.

Your inverse MBC equation looks correct. When you time-integrate the rotor speed to get the azimuth angle, does the azimuth angle match well between the linear and nonlinear models, or is there some phase shift due, e.g., to an initial nonzero azimuth?

Best regards,

[Hadi.Hammoud]

Dear Dr. Jonkman,

First of all, thankyou so much for your help, I can now make much more sense of the entries in the A matrix. However, I still do not fully understand why, physically, does the effect of the azimuth state cancel out, could you elaborate more on that?
Concerning the flapwise moment issue:

1. I have verified that the rotor azimuth angle (integral of rotor speed) responses of the linear and non-linear models match (so no phase shift).
2. Another observation is that the C and S components of the flapwise moment output of the linear model are quite large, almost equal to the 0 component, is that possible?
3. Also, in the non-linear model, the flapwise moment is not a perfect sinusoid, is that reasonable or could it be the result of some underlying error?

Thanks again for the cooperation.

[Jason.Jonkman]

Dear Hadi,

The effect of the azimuth state should cancel out because you’ve aizimuth-averaged the MBC-transformed matrices, which eliminates any remaining periodicity (azimuth dependence) in the system.

I’m not sure why the C and S components are so large, but I see in your results that they are large even before your step in blade-pitch, so, perhaps this is a start-up transient that needs to die out before using the linear model (it may be possible to reduce the start-up transient) by using better initial conditions of the states in the linear model).

Best regards,

[Hadi.Hammoud]

Dear Dr.Jonkman,

I do not think that the large C and S components are transients since as can be seen in the figure below, they reach a steady-state value at which they settle. I do not understand how can the responses for the rotor speed and the foreaft moment agree so closely while that of outputs in the rotating frame do not. Could this be an indicator of a more fundamental error?


Also, concerning the flapwise moment response in the non-linear simulation, is it natural to observe a non-perfect sinusoid, or does that indicate an underlying error? Refer to the figure below for reference.

Sorry for annoying you with my repetitive questions, thankyou so much for everything.

[Jason.Jonkman]

Dear Hadi,

The nonlinear response is not surprising to me and likely depends on the shaft tilt, deflections, and gravity loads. I would not expect you’d match this exactly with the linear model, but the large difference in amplitude you are showing is concerning.

Regarding your plots of the C and S components, where is the large nonzero mean coming from? If the output equation is written as y = Cx + Du and the states and inputs are initialized to zero (are they?), why would the value at time zero be nonzero? Where have you added the operating point (OP) values to the perturbations derived from the linear model output?

Best regards,

[Hadi.Hammoud]

Dear Dr.Jonkman,

It turns out that I have made the silly mistake of adding the OP values before applying the inverse MBC transformation, which clearly does not make sense and was causing this discrepancy. After correcting that and adding the OP value for Mflapb1 after the inverse MBC transformation, I obtained the following:

The discrepancy is now much smaller, but the outputs still do no match like for other output variables. What’s your outlook on that? i.e. is it safe to say that the linear model is valid or should I look for further errors?
On another note, is it reasonable to azimuth-average the vectors x_op,y_op, and u_op obtained from each of the 36 FAST linearization files?
Thankyou again, your help was invaluable in pinpointing this mistake.

[Jason.Jonkman]

Dear Hadi,

I now don’t see any oscillation in the flapwise moment from the linear mode before the step change in blade-pitch. Do the operating points (OPs) not have an azimuth variation or are you only using the azimuth-averaged OPs?

You can azimuth-average the OP to get the mean value over one rotor rotation, but that doesn’t mean this is a useful to do. I would expect that you’d want to combine the perturbation from the linear model with the azimuthally varying OP, at least for the blade-related quantities.

Best regards,

[Hadi.Hammoud]

Dear Dr.Jonkman,

Your latest comment opened my eyes to a fundamental aspect that I was missing: the periodicity of the operating point should be accounted for when adding it to the disturbance. With the aid of a look-up table including the OP data retrieved from the FAST linearization file, I have managed to add up the periodic OP value to the disturbance based on the azimuth angle. After doing that, we obtain a very nice correspondence between the linear and nonlinear models as can be seen below.

I now understand the linearization process much better. One final question: do you think there is a better way to add the periodic OP values than using a look-up table? Or should I just increase NLinTimes?
I cannot stress how helpful your comments were, thankyou so much.

[Jason.Jonkman]

Dear Hadi,

I’m glad that your linear model is now working as expected! Hopefully this commentary helps other users of the linearization functionality as well.

Regarding the OP, typically I recommend setting NLinTimes = 36, which is a linearization every 10-degrees of azimuth rotation as a reasonable compromise between accuracy and computational expense of the linearization process. You can interpolate if you want the OP at a finer azimuth revolution. But you could always change NLinTimes too see its effect on the result.

Best regards

[Mohamed.Fekry]

Dear All

I have a question regarding the azimuth of the linear model.
I am trying to find it using the integration of the rotor speed of the linear model. However, it gives me an endless ramp (As the rotor speed is a constant ) which exceeds the 360 deg of the rotation.
Can you guide me towards finding this angle correctly from the linear model?

I tried to make the same comparison but using the azimuth that was calculated from the nonlinear model with the same DOFs as mentioned above

3 DOFs were activated: GenDOF, FlapDOF1 (first flapwise mode), and TwFADOF1 (first fore-aft mode); that is a total of 10 continuous states.

But the results were frustrating. Additionally, the azimuth calculated from the linear model is always zero. The connection diagram of the linear model is shown below.

This linear model is calculated by azimuth averaging after applying the MBC transformation.
The operating point is added to the states after the application of MBC transformation from the rotating frame to the non-rotating frame as follows:

function T_inv = fcn_T_inv_OP(azmth) T = 2/3*[0.5 0.5 0.5; cosd(azmth) cosd(azmth+120) cosd(azmth+240); sind(azmth) sind(azmth+120) sind(azmth+240)]; T_inv =[eye(1) zeros(1,8); zeros(3,1) T zeros(3,5); zeros(2,4) eye(2) zeros(2,3); zeros(3,6) T]; end

The input pitch angles are transformed to the non-rotating frame using the following transformation:

function T_inv = fcn_T_inv(azmth) T_inv = 2/3*[0.5 0.5 0.5; cosd(azmth) cosd(azmth+120) cosd(azmth+240); sind(azmth) sind(azmth+120) sind(azmth+240)]; end
The output of the blade RootMFlps are transformed to the rotating frame using the following transformation:

function T = fcn_T(azmth) T=[1 cosd(azmth) sind(azmth); 1 cosd(azmth+120) sind(azmth+120); 1 cosd(azmth+240) sind(azmth+240)]; end
Again, all the used azimuth angles are calculated from the nonlinear model.
The lookup table is used to reflect the azimuth dependant operating points.
The order of the states are:
ED 1st tower fore-aft bending mode DOF (internal DOF index = DOF_TFA1), m
ED 1st flapwise bending-mode DOF of blade 1 (internal DOF index = DOF_BF(1,1)), m
ED 1st flapwise bending-mode DOF of blade 2 (internal DOF index = DOF_BF(2,1)), m
ED 1st flapwise bending-mode DOF of blade 3 (internal DOF index = DOF_BF(3,1)), m
ED First time derivative of 1st tower fore-aft bending mode DOF (internal DOF index = DOF_TFA1), m/s
ED First time derivative of Variable speed generator DOF (internal DOF index = DOF_GeAz), rad/s
ED First time derivative of 1st flapwise bending-mode DOF of blade 1 (internal DOF index = DOF_BF(1,1)), m/s
ED First time derivative of 1st flapwise bending-mode DOF of blade 2 (internal DOF index = DOF_BF(2,1)), m/s
ED First time derivative of 1st flapwise bending-mode DOF of blade 3 (internal DOF index = DOF_BF(3,1)), m/s

The comparison results are shown below: the blue line represents the linear model and the yellow line represents the non-linear model
The Generator speed VS Simulation time ( rpm VS seconds)

The Tower fore-aft VS Simulation time ( m VS seconds)

Regards
M.Fekry

[Mohamed.Fekry]

Continuing from above

The Blade RootMFlp1 VS Simulation time ( kN.m VS seconds)

The A matrix azimuth averaged after MBC application is shown below (Note that: The second column and raw belong to the generator Azimuth state which should be eliminated) :

Regards
M.Fekry

[Jason.Jonkman]

Dear Mohamed,

I’m not following everything you are doing.

How are you calculating Azimuth in the linear model? (I see it used, but not calculated, in your Simulink model.)

The linear model should only involve perturbations about the operating point, but it seems like you are using the sum of the operating point and perturbations.

Also, it looks like you have two functions called fcn_T_inv().

Best regards,

[Mohamed.Fekry]

Dear Dr. Jason

Thank you for the reply.
The used Azimuth here is calculated from the nonlinear model because the integration of generator speed of the linear model gives me the figure below (ramp) and the results become not logical if it is used instead:
Linear Azimuth (Deg) Vs. Time (sec.) (rotor speed integration)

Regarding fcn_T_inv, I am sorry for this confusion, I agree with you, There are two of them, But, they are in different locations and do not reflect each other. For avoiding and confusion. I renamed each one of them.
fcn_T_inv_OP → Transform the operating point into the non-rotating frame (Dimension depends on the order of the operating points).
fcn_T_inv → Transform the pitch angles into the non-rotating frame (3 by 3 Matrix).

I applied your advice regarding the perturbations and the results are enhanced but not perfectly. The linear model is linearized at a rotor speed of 12.1 rpm with a pitch angle of 14.755 deg and wind speed of 18 m/s. I made perturbation of the pitch angle by 1.45 deg. The generator speed and the Tower fore-aft become reasonable. However, the Blade RootMFlp1 still has a large deviation.

Is the azimuth angle may be the cause of that?

The yellow line for the nonlinear model. The blue line for the linear model.
The Blade RootMFlp1 VS Simulation time ( kN.m VS seconds)

The Generator speed VS Simulation time ( rpm VS seconds)

[Mohamed.Fekry]

Continuing from above
The yellow line for the nonlinear model. The blue line for the linear model.

The Tower fore-aft VS Simulation time ( m VS seconds)

If I calculated the azimuth of the linear model by continuous integration of generator speed divided by gear ratio

It gives me a continuous ramp which makes all the results not logical. So, I always use the azimuth of the non-linear model in the transformation matrices between the rotating and non-rotating frame. The comparison between the linear and the nonlinear model in azimuth angle is shown below
The azimuth VS Simulation time ( deg VS seconds)

I hope everything becomes clear.

Regards
M.Fekry

[Jason.Jonkman]

Dear Mohamed,

I would not expect that the linear and nonlinear model would produce identical azimuth angles because the rotor speeds likely have some small differences. So, you should derive the azimuth angle from the linear model if you need to use it in the linear model. To keep the azimuth angle within the 0-360-deg range, you can subtract 360 deg whenever the integration exceeds 360 degrees.

However, your azimuth comparison shows a bigger difference – the slope of the curves are not the same, suggesting some error in the calculation or the use of different rotor speeds.

Best regards,

[Mohamed.Fekry]

Dear Dr Jason

Thank you very much for your response.
I already derive the azimuth from the linear model and keep its values between 0-360 deg.
Now, the comparison of the azimuth angle become reasonable as follows:
The azimuth VS Simulation time ( deg VS seconds)
The yellow line for the nonlinear model. The blue line for the linear model.

.

However, when I used this linear azimuth angle in the transformations of the linear model, I did not notice any enhancement in the large deviation in the Blade RootMFlp1.
The Blade RootMFlp1 VS Simulation time ( kN.m VS seconds)
The yellow line for the nonlinear model. The blue line for the linear model.

Is there any suggestion regarding this?

Regards
M.Fekry