I have updated the issue on github with additional problems that should be resolved.
I replaced this line [CampbellData{i_case}] = campbell_diagram_data(MBC_data{i_case}, BladeLen, TowerLen, , strrep(newFSTName,‘.fst’,‘.MBD.sum’));
to
[CampbellData{i_case}] = campbell_diagram_data(MBC_data{i_case}, BladeLen, TowerLen, strrep(newFSTName,‘.fst’,‘.MBD.sum’));
With this replacement I am now able to write the MBD.sum.xlsx file. I thought the use of was a mistake.
I still have problems generating the campbell diagram. I cannot use the xlswrite and csvwriteCampbell function that you shared with me on github. So I used the writetable function to write the xlsx file. But I get some other new errors (see below) while using the Plot_CampbellData.
-----Error output from matlab--------
Index exceeds matrix dimensions.
Error in Plot_CampbellData (line 93)
CampbellPlotData.lineLabels = txt(2: nLines+1, 1);
Error in runCampbell (line 266)
[num,txt,CampbellPlotData] = Plot_CampbellData(XLSname,wkSheetName);
when I was trying to run mbc file with more than two different .lin files this error message was displaying.
[MBC]=fx_mbc3(‘5MW_Land_DLL_WTurb.1.lin’,‘5MW_Land_DLL_WTurb.2.lin’,‘5MW_Land_DLL_WTurb.3.lin’,‘5MW_Land_DLL_WTurb.4.lin’,‘5MW_Land_DLL_WTurb.5.lin’,‘5MW_Land_DLL_WTurb.6.lin’)
Error using fx_mbc3
Too many input arguments.
I am having trouble reproducing the Campbell diagram of NREL 5MW land based turbine using the test case ‘5MW_Land_BD_Linear’ available within the OpenFast dev branch. I have tried to use the matlab runCampbell script with some modifications to make it usable on Mac and generate the excel file to plot the Campbell diagram. The results are very different from the once reported by you. I am trying to understand if there is a problem with the matlab script or the input parameters used in the fst file. Could you please answer the following queries to help me diagnose the problem?
Does the input parameters, like total simulation time of 5s, timestep of 0.001, NLinTimes=1, LinTimes=0 used in the .fst file reproduce the results? Do I need to increase the total simulation time and change the first linearization step in LinTimes?
I am able to generate the excel file by adapting runCampbell script for different rotor speeds but the worksheet where the modes are mapped to the natural frequencies have some index numbers to be zero for all the rotor speeds. Could you help me understand the algorithm used in IdentifyModes function?
I have also tried using the fx_mbc3 and the CampbellDiagram excel sheet shared on the forum to cross-check the results. But the mode descriptions, in particular for the blades, are not very clear. Also it was difficult to copy the output data (DescStates, MBC_ModeShapeMagnitude, MBC_ModeShapePhaseDeg_MBC_DampingRatio) from fx_mbc3 as they were all in descending order unlike what is mentioned in the excel file comments. Is this data not sorted inside the fx_mbc3 script?
The model 5MW_Land_BD_Linear.fst is not set up correctly to generate a Campbell diagram; this test was simplified for testing the compiled code, not producing useful results. For example, to generate a proper Campbell diagram, you’ll need to find a periodic steady-state solution before linearizing and linearize at multiple azimuth steps to cover the full rotor rotation (requiring increasing TMax and setting NLinTimes = 36 with appropriate LinTimes). And due to the large number of DOFs in BeamDyn, you’ll also need to change the precision of the linearization output (e.g., use OutFmt = “ES17.9E3”) so that the eigenanalysis produces good results.
The IdentifyModes function only partially works. It is probably better to identify the modes manually, at least at first.
You don’t need to use the old CampbellDiagram spreadsheet and manually paste in the data if you are using the runCampbell.m script.
Please note that the OpenFAST linearization process is being improved now. The approach will be more straightforward–including automated calculation of steady-state solutions and mode shape visualizations–once the following OpenFAST pull request is merged in: github.com/OpenFAST/openfast/pull/373.
Thanks for the information. In fact, I did tried several things like a) using TMax and DT values similar to the pull373 branch input parameters used in 5MW land mode shapes case. With TMax=5600 and Dt=0.001 the simulation runtime is long. Except for the results of the 1st tower FA, all other values do not match with previous results (please see the attachment for the figure). Furthermore, simulations for some rotor speeds have aborted (see below for the output from openfast with the fatal error). I tried to rerun the aborted cases to double-check if there is any problem. I think the simulation aborts at the last timestep when the linearization has started as the steady state solution is not found. b) I have seen here in the forum that researchers used previous versions of fast with larger timestep like 0.0125s and smaller simulation time of 240s to get the right results. However, if I use DT>0.002 I see convergence errors and simulation aborts without reaching the linearization step c) I have also tried using the dev branch to see if I get any better results with simulation time like 1000s and Dt = 0.001 but the simulation aborted at after ~620s with a fatal error. In all these cases I was using CompElast and CompServo flags ON. Why is this difference between openfast and previous fast versions? Could you please suggest any reasonable values that can work for generating the Campbell diagram with the OpenFast dev branch?
Yes, I was initially using only the Matlab script to automatically generate the Campbell diagram without copying data from MBC3 output. But since I was not seeing the correct results, I started using the excel sheet manually copying the MBC3 output to cross-check.
Regarding the NLinTimes and Lintimes. In all cases except 0 rpm when NLinTimes=1, Is it right the runcampbell script does some averaging to produce the Campbell diagram?
The use of BeamDyn generally requires a small time step; DT=0.001 s is what NREL decided was appropriate for a numerically converged solution. You won’t be able to increase this without disabling model features (e.g., disabling degrees of freedom or reducing the order of the finite element in BeamDyn). Simulations using DT = 0.0125 s likely used ElastoDyn to model the blades, not BeamDyn. I’m not sure why you want to use BeamDyn instead of ElastoDyn to model the blades, but switching to ElastoDyn will speed up the simulations a lot. From your screenshot, it looks like perhaps you already switched to using BeamDyn?
Your results are looking OK for the tower mode; I also see the 1st blade flapwise collective, regressive, and progressive modes in your figures (although these are labeled incorrectly in your plot), which look correct up to 10 rpm. I would assume that the steady-state solution is not converged at the 12 and 14 rpm speeds, or perhaps the labels are not shown correctly at these speeds (you’ll have to manually identify the modes at each speed).
Correct. The runCampbell.m script calls MBC3, which azimuth-averages the matrices after applying the Coleman transformation at each azimuth step and before performing the eigenanalysis.
Thanks for your clarification on the simulation time and the timestep. There was no particular reason to use BeamDyn. I started using the test case 5MW_Land_ModeShapes from pull373 branch where BeamDyn was used by default. I did not change it assuming the case was already tested based on the discussions in the forum. Yes, I did switch to ElastoDyn with TMax=10000s and DT=0.00625. However, some of my simulations (rotor speeds=8, 10 rpm) are still aborting with fatal error (Unable to find steady-state solution). I think TMax=10000 s is sufficient enough to reach steady state. I am unable to understand why the simulation fails for these two speeds. Do you think TMax should be further increased?
I think runCampbell script has also some issues identifying the modes correctly because it does not produce the same results as the manual approach to generate the Campbell excel sheet from the MBC3 output.
I’m sure TMax = 10000 is more than sufficient to reach a periodic steady-state. What settings have you defiined in the “Linearization” section of the OpenFAST primary input file? You should be able to tell from looking at the time series if the solution is converging towards a steady state or not; what does the time series of the response look like?
The runCampbell.m does not identify the modes for you; it simply runs the many OpenFAST simulations, post-processes the linearization output files using MBC3, and gathers the data into the MS Excel spreadsheet. It is still up to the user to properly identify the modes.
I got the 5MW_Land_ModeShapes case files from links available at github.com/bjonkman/r-test/blob … ModeShapes. However, now I see the link is not working anymore. I think this could be due to some recent updates within OpenFAST github repositories.
I can see 5MW_Land_ModeShapes.fst available at github.com/OpenFAST/r-test/blob … Shapes.fst and using ElastoDyn. Also, I did use this test case from the above link but I still get the error saying steady-solution not found. Maybe the code from pull request 373 has some issues running this test case for rotor speeds 8 and 10 rpm?
Please see attached the fast, elastodyn, and servodyn files that I used for 8 rpm rotor speed.
See below for the my settings in the “Linearization” section:
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.001 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 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]
The “Linearization” section you included in your forum post differs from the “Linearization” section in your attached zip archive. Which one are you using?
TrimCase = 3 doesn’t really function when you disable aerodynamics (because changing the blade-pitch in this case doesn’t effect the torque), as it appears that you have in your attachment. TrimCase = 2 should still work, and settle out on zero torque for all speeds without aerodynamics. If you’ve done this, and it is still not converging, I would suggest plotting the time series to see what is going wrong.
Sorry for the confusion. At the moment I am using the linearization section present in the fast file in the attachment. The fast file in the zip file fails with the fatal error.
Please see attached the rotor speed time series for the first and last 1000s. I see that for the first 1000s the rotor speed is fluctuating around 8 and for time series between 1000-9000s it is fluctuating slightly above 8.005. In the last 1000s the rotor speed drops from 8.05 to 7.975.
It appears that there is not enough damping in the system to reach a steady state and I’m not sure I understanding what is happening at 10000 s. Reducing TrimGain may help to reduce the oscillations. That said, the solution is likely sufficiently converged by 200 s, so you can probably stop it there (increasing TrimTol a bit). Without aerodynamics, it is likely also fine to disable the trim solution altogether (TrimCase = 0) and simply linearize after 200 s or so.
As you suggested I have tried with TMax=200s, TrimTol=0.01, TrimCase=2 (there is no option to set TrimCase=0), and TrimGain=100. This time the simulation for rotor speed=10rpm finished without any problem. But for the rotor speed=8rpm, I still keep getting the same problem. So I reduced the TrimGain to 80 and increased the TrimTol to 0.1. This time the simulation completes without any errors. However, I have to check if the linearization results with these settings are still reliable.
Just to clarify. If CalcSteady = False, then I need to set LinTimes from 200s. Is it okay to choose any 36 (NLinTimes) for LinTimes starting from 200s?
Correct, but you’ll want to set the LinTimes such that you cover a full rotor revolution of the rotor over NLineTimes steps, i.e.:
deltaT = (60 s/min)/( RotSpeed*NLinTimes ) where RotSpeed is in rpm
LinTimes(1) = 200 s (or whatever is needed to be in steady state)
LinTimes(2) = LinTimes(1) + deltaT
LinTimes(3) = LinTimes(2) + deltaT
etc.
Can someone explain to me where does it come from? Is it related to the harmonics of the turbine? Is it normal? And if it isn’t, how to overcome this issue?
I would guess there is some excitation happening at the period of the response, which I estimate at about 3.6 s or 0.278 Hz, e.g., wave excitation or a harmonic of the rotor speed. What incident waves (WaveTp) have you defined? What is the rotor speed you are simulating?
