Frequency and time response for linearized model

[QusayHawari]

Hi,

Thanks for your recommendation about increasing TrimTol. The simulations runs well now using a tolerance of 0.015. Under CalcSteady = T, the first steady state solution took place after 4800s. This high time has to do with the low TrimGain of 0.0001 I am using.

I performed frequency response analysis on the linearized state space models, and I got the following:

Frequency Reponse

Uniform wind input of 15m/s. Four linearizations that took place after 4800 seconds, and are very close to each other (approx 10s apart). The response is from the command collective pitch angle to the platform pitch. No wave excitation matrices were included. All 6 platform DOFs were on. My question is that shouldn't all responses be identical or closely related to one another since steady state is reached, and the same wind is applied. I also tried running a 7000s simulation without using CalcSteady, and linearized the model well into it's steady state condition at 5000s,5500s,6500s,and 7000s but the linearized Bodes are still different.

The other question is that the frequency response is not found if I plot the Bode from the commanded collective pitch to the output rotor speed. I have also tried the individual pitch angles, but nothing was shown.

Another thing I have checked was the time response of the platform pitch. The natural frequency for the platform pitch (IEA15MW UMaine ActiveFloat/Semi-sub) is around 0.03Hz. From the .out file, I plotted the pitch, and it looks like the frequency is a little lower (around 0.01Hz).

Thanks in Advance.

[jjonkman]

Dear @QusayHawari,

A few comments/questions:

• What DOFs have you enabled in your OpenFAST model?
• Are you applying MBC3 before computing the Bode diagram?
• Have you confirmed that a reasonable steady-state solution has been reached, e.g., by plotting the time series response of various outputs?
• What do you mean when you say, "the frequency response is not found" and "nothing was shown"?
• Regarding the IEA Wind 15-MW wind turbine atop the UMaine semisubmersible, this is probably a better question to ask here: https://github.com/IEAWindTask37/IEA-15-240-RWT/issues.

Best regards,

[QusayHawari]

Hi,

I asked the question where you recommended, but so far I didn't get a response.

Regarding the other questions you had:

• Those were the DOFs and initial conditions I set:

• I did't use MBC3, what I did was simply import the .lin files to Matlab, and plot the bodes. I was looking for a documentation to integrate simulink with OpenFAST on Mac, but I could't find any, so at the moment I am trying to do that on a windows virtual machine. If there's any documentation for Simulink linking for Mac, please let me know.

• Regarding steady state, the platform motions did't reach steady state even after 5000s run. I thin I need more than 10000s for all platform motions to settle. The Calcsteady just checks the rotor speed tolerance, and based on that the linearization is performed. Any idea on how to speed up the response besides increasing TrimTol?

• What I meant with the response not found is that the bode plots do not show any response (gain and phase). i.e, blank figures.

Finally, I wanted to test the linearization on the virtual machine, so I used the exact same input file and settings, but got different results, that is, the linearization took place approx 200s before what I got using Mac as seen below:

Thanks in Advance, Qusay.

[jjonkman]

Dear @QusayHawari,

Just a few comments:

• You have states (blade flapwise bending) in the rotating frame, so, I would recommend applying MBC3. The response is likely also periodic, so, I would recommend linearizing at a number of azimuth steps around a revolution (typically we recommend linearizing 36 times around the revolution, or once every 10 degrees) and then azimuth averaging. And with the generator DOF enabled, it is important after azimuth averaging to eliminate the generator azimuth state, as has been discussed many times on our forum: https://wind.nrel.gov/forum/wind/. I would fix these issues before post-processing the linearization output. It may also help to simplify the model to just a few key DOFs (e.g., platform-pitch and generator) to debug and build complexity in steps (with only these DOFs enabled, MBC3 is not needed as long as you don't care about blade outputs).
• Regarding the steady state, I would expect that the response will be periodic (due to rotor rotation) rather than a pure steady state (we typically use the term "periodic steady state"). And it is not just the rotor speed that is checked for convergence; instead, it is the periodicity of the outputs that is checked. I would be interested to see the time series of various motion outputs plotted.
• There are likely small differences in the OpenFAST output when comparing across platform and when using different compiler settings.
• I'm not sure what Simulink has to do with the linearization analysis you are using, but on Mac, presumably you are using CMake to compile OpenFAST; is that correct? There is an option for compiling OpenFAST for Simulink using CMake, as documented here: https://openfast.readthedocs.io/en/main/source/install/index.html. Otherwise, the process to use OpenFAST within Simulink is the same on both Mac and Windows.

Best regards,

[QusayHawari]

Hi,

Thank you for the clarification,

I sum up what I did in the following points:

• I turned FlapDOF1, GenDOF, and TwFADOF1 on, and the rest off.
• Set NlinTimes = 36 to compute linearization every 10-degs.
• Used the .m file fx_mbc3.m to generate the averaged state space mats based on 36 linearizations.

I got the following result from MATLAB.

Based on what I understood:

• We need to transform the states from rotating to fixed so that we can visualize the individual blade performance.
• We also need to generate an average state space representation based on the .lin files.
• Finally, we need to eliminate the generator azimuth state DOF from the state space matrices.

All the three points can be done by running the MBC3 file? if so, then my states are the averaged ones, i.e (AvgA...etc)? and my last question is that why should we eliminate the generator azimuth DOF from the model?

Thanks in Advance, Qusay

[jjonkman]

Dear @QusayHawari,

Yes, your approach sounds correct. The _fxmbc3.m script performs the MBC3 transformation and will azimuth average across all 36 *.lin files, but it will not eliminate the generator azimuth state for you; you'll need to eliminate the generator azimuth state manually from the azimuth-averaged matrices output from MBC3 (AvgA, etc.). You should eliminate the generator azimuth state, because once azimuth averaged, the resulting matrices will no longer be periodic and there should be no more azimuth dependence. In effect, the azimuth-averaged state-space system becomes linear time invariant (LTI), which is likely how you want to use in your application.

Best regards,

[QusayHawari]

Hi,

Thank you for the clarification.

Looking at the column related to the GeAz state from AvgA state matrix, the values seem to be quit large relative to other states. From what I understood is that the MBC3 functionality removes periodic dependencies and act as a filter. Shouldn't such a column have insignificant effect that we remove it from the state space representation?

Thanks in advance, Qusay.

[jjonkman]

Dear Qusay,

A similar question was recently discussed on our forum--see: https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=2823&p=18214. My comments and questions particularly from my post dated Oct 18, 2021 (and related posts) apply equal well to your results.

Best regards,

[QusayHawari]

Hi,

Thanks for sharing the link.

Focusing on the lower left quadrant of the A matrix, I can see that entry (2,2) is lower than other diagonal entries (-0.1556), but other entries in the quadrant related to GeAz are large (250.6746, -406.5812, and -3.7712). is that alright?

I am asking because I wanted to make sure that the extracted matrices are in full steady state operation.

I performed 2 other linearizations with the following DOFs on. TrimTol = 0.005.

Case 1: FlapDOF1, GenDOF, TwFADOF1, TwSSDOF1, EdgeDOF, PltfmPDOF.

Case 2 FlapDOF1, GenDOF, TwFADOF1, TwSSDOF1, EdgeDOF, PltfmPDOF, PltfmSwDOF.

For the first case, all time response signals settled at a final value, but for case two, linearization (CalcSteady = True) occurred before the platform sway reached steady state as shown below.

Could this be the reason, and should I rely on setting calcSteady = F while using a long simulation time for the sway (or any slow mode) to settle?

Thanks, Qusay.

[jjonkman]

Dear Qusay,

You may need to change the value of TrimTol to ensure that solution is in steady-state before linearizing. I agree that the sway value you are plotting does not look like it is steady state yet. I would recommend plotting the time series to be confident that the solution is in steady state before linearizing. And by plotting the time series, you can be better informed how to set initial conditions to minimize the length of time it takes to reach a steady-state solution in a subsequent simulation.

I'm not sure why the values of A(8,2), A(9,2), and A(10,2) are so large in your case as well as in the example on the forum (perhaps this is because of ill-conditioning in M^-1 for these terms?), but I would still say the value of A(2,2) is likely the most important. And as I suggest in the forum topic linked above, I suggest plotting the periodic variation of the MBC3-transformed elements of A to verify the periodicity of the solution and to see how how well the azimuth-averaged values capture the mean effect.

Best regards,

[QusayHawari]

Hi,

Thanks as always for the clarification :)

As you recommended I plotted the entries related to the GeAz state from all 36 A mats vs the azimuth angle, and are shown below. It looks similar to the results in https://wind.nrel.gov/forum/wind/viewtopic.php?f=30&t=2730&p=17489&hilit=azimuth#p17489

It was mentioned that such spikes were related to the tower influence, but in the forum, the model didn't include any tower DOF and was based on just the GenDOF, yet the spikes were there. Does the tower influence refer to something other than the tower DOFs?

Finally I have some questions regarding validation between the linear and non-linear models. First is the file .out from the linearization simulation based on the linearized system or are those outputs tabulated for the nonlinear model(basically an open-loop nonlinear simulation)? Finally, if I wanted to apply a disturbance to the non-linear model (say 1 deg pitch step or 1m/s step wind), how can I perform that?

Thanks in Advance, Qusay.

[jjonkman]

Dear @QusayHawari,

By "tower influence", I was referring to on the forum is the aerodynamic influence of the tower on the rotor aerodynamic loads. These can be disabled by setting TwrPotent = TwrShadow = 0 in AeroDyn. I would recommend disabling this tower influence to eliminate the "spikes", which appear to be corrupting the azimuth averaged matrices.

The *.out file output from OpenFAST is always from the nonlinear response. This is the file you should use to track progress on convergence of the steady-state solution.

You can simulate a step change in pitch angle through the override pitch maneuver inputs in ServoDyn or via a user-defined pitch controller.

You can simulate a step change in wind speed using WindType = 2 in InflowWind and an appropriate uniform wind data file.

Best regards,

[QusayHawari]

Hi,

Apologies in advance for the lengthy inquiry.

I applied a wind step from 15-16m/s, and compared the fixed frame outputs, with the nonlinear response. I used the Avg state mats extracted from the MBC3 .m file.

The outputs of the linear closely resemble that of the nonlinear as shown below.

The next step was to validate the rotating frame outputs. For that, I relied on the RootMby for blades 1-3.

I applied the inverse MBC3 on the output fixed RootMby by following : q_b = q0 + qc x cos(phi_b) +qs x sin(phi_b), where b is the number of blades, and phi_b = phi + (120 deg)*(b-1). I obtained phi by integrating the output rotor speed. I did that for all three blades to compare the nonlinear individual RootMby for blades 1-3.

I then compared the simulink output with that of the nonlinear. The results were as follows:

These results where similar to that in https://wind.nrel.gov/forum/wind/viewtopic.php?f=4&t=2823&p=18214 .

From what I understood, the difference in amplitudes can be eliminated by including the OPs of the linear RootMby from all 36 .lin files (before applying MBC3). Those 36 X 3 (for the 3 blades) can then be added to the output of the RootMby1,2,3 from simulink by interpolating the current azimuth angle (phi).

Please refer to the simulink block diagram shown below.(I hope it's clear), y1,2,3 are the interpolated OPs (minus the mean OP for RootMby1,2,3 from the 36x3 OPs) according to the input Az (integral of RtSpeed).

When doing that, I got the following results:

Thanks in advance Qusay

[jjonkman]

Dear @QusayHawari,

Overall your approach makes sense, but I'm guessing one of the details is incorrect. Before the step-change in wind speed, which I assume is the operating point you are linearizing about, the perturbed inputs, states, and outputs of the linear model should all be zero. But to compare the linear and nonlinear solutions, you need to add the periodic operating point to the linear perturbation. In your last plot, you show that the linear output is zero (which I assume is just the perturbation), but your nonlinear solution is periodic, which I assume is the periodic operating point. So, it doesn't appear that you are summing the linear perturbation with the periodic output properly, because these should match before the wind speed is perturbed.

Best regards,

[QusayHawari]

Hi,

I tried another rotating frame output (TipRDyb1,2,3) to double check.

I understand that the output of the linear model is just the perturbation, so I added to that output the mean TipRDyb, which I obtained from calculating the mean of all 36 OPs.

Looking at the output comparisons for blade 1, it looks similar to that of RootMyb shown in the previous post.

If we zoom in to see the difference in amplitudes before accounting for the periodicity of the OPs, it can be seen that a difference of 0.33 degs is there between the linear and nonlinear.

If I plot the Azimuth Vs the TipRDyb1(minus the average OP), we can see the periodicity is there, however, the amplitude is higher than what we want, which is something close to 0.33 degs.

I think this is why the output does not match that of the nonlinear, because of you account for that periodicity shown is the last picture above, then the outputs will be similar to what I showed you yesterday.

Just to make sure, I get the OPs from all the 36 *.lin files for TipRDyb1,2,3?

Thanks in advance, Qusay.

[jjonkman]

Dear @QusayHawari,

I believe your issue is that you are only using the mean value of the OP. The OP is periodic, so, you should add the periodic value of the OP--corresponding to the value of the OP at the given azimuth angle--to the perturbation from the linearized solution. Before the step change in wind, when all perturbations are zero, the solution from the linear model should then match that of the nonlinear model and I would think the results would be close after the perturbation as well.

Best regards,

[QusayHawari]

Hi,

My output azimuth angle(integral of RtSpeed) from the linear model is zero before the step change in wind, hence only one OP is used before the step which are the OPs at 0 azimuth.

Azimuth for each blade: Shouldn't that be the case?

Qusay.

[jjonkman]

Dear @QusayHawari,

Just like other variables in the linear model, the rotor speed output from the linear model is also a perturbed value. This perturbation should be added to the OP rotor speed before integrating to calculate the azimuth angle.

Best regards,

[QusayHawari]

Thank you for your help.

I now moved to a different step, and so I have a question related to getting the s-function to use in Simulink using Mac OS (I created a new issue :) ). I will close this issue now, and if I need any help related to said issue I will open a new one.

Qusay