OpenFAST: Linearization of NREL 5MW Onshore Turbine

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

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,

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.

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,

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?
C Component of Mflapwise.JPG
S Component of Mflapwise.JPG
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.
FlapWise.JPG
Sorry for annoying you with my repetitive questions, thankyou so much for everything.

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,

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:

Mflap1.JPG

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.

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,

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.

Mflap1 Corrected.JPG
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.

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

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

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

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

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,

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)

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

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,

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

Dear Mohamed,

What input perturbation have you applied to the linear model? Does the output from your linear model match that of the nonlinear model when no input perturbations to the linear model are specified, such that all states, state derivatives, and input perturbations are zero? This should be the case because the output perturbations from the linear model are then zero, and when the periodic OP is added the perturbation, the result should match the nonlinear model.

Best regards,

Dear Jason

Thank you for your explanation.
However, I found less deviation between the linearized model and the nonlinear model, although I set all perturbations to zeros. I add the periodic OPs that are retrieved from the linearization txt file (*.lin). The operating conditions of the nonlinear model are the same as the linearization point. Moreover, I made another simulation with accounting the periodicity of A, B, C, and D matrices, but, it gives me the same results.
The Blade RootMFlp1 VS Simulation time ( kN.m VS seconds).
The yellow line for the nonlinear model. The blue line for the linear model.

The connection diagram is shown below

Regards
M.Fekry

Dear Mohamed,

I still have not fully understood your Simulink model, but if the perturbations zero, then the linear and nonlinear model should agree if the OP from the linear model is derived from the conditions of the nonlinear model. Are you adding the periodic OP to perturbation (which are all zeros for now) when comparing the linear and nonlinear models, including inputs, states, and outputs?

Best regards,