Hello Every Body,
I am modeling a rigid offshore floating wind turbine, and like to find AeroDynamic forces on the rotor while the rigid body has harmonic motions in the surge direction in the form of acos(wt) , where “a” is the amplitude of motion, and “w” is the radial frequency of the motion.
To model a rigid turbine, in the primary.FST file, in the part “FEATURE FLAGS” I deactivated all the degrees of freedom except CompAero.
To give harmonic oscillations to the structure in the surge direction, in the file FAST.f90 in the subroutine solver, I changed the value of Q(displacement), QD(velocity),QD2(acceleration) ,QT,QDT,QD2T for the first degree of freedom (surge DOF) to a harmonic function. In details:
The subroutine solver starts with the 2 lines: “QT = Q (:,IC(1))” " QDT = QD(:,IC(1))" before these 2 lines I added the lines :
QT (1) =amp1cos(w(step+1)DT)
QDT(1) = -amp1wsin(w(step+1)DT), where amp1 is the amplitude of motion in surge, and w is the radial frequency.
The next lines in the subroutine solver are related to RK4 time integration method. In each substep, I used the same expressions as above to redefine the displacement and velocity in the surge direction. In the last substep of RK4,where the final values of Q,QD are calculated I added the lines:
Q (1,IC(NMX)) =amp1cos(w*(step+1)DT)
QD(1,IC(NMX)) = -amp1wsin(w(step+1)*DT)
The next lines in the subroutine solver are related to Adams-Bashforth predictor, in this part I added 2 lines the same as the above lines. After the predictor, there is the line “QD2(:,IC(NMX)) = QD2T” before this line I added:
In the Corrector (Adams-Moulton) part, I added the lines:
Q (1,IC(NMX)) =amp1cos(w*(step+1)DT)
QD(1,IC(NMX)) = -amp1wsin(w(step+1)*DT)
After the corrector part there is the line: “QD2(:,IC(NMX)) = QD2T” before this line I added:
I really appreciate it if anyone could tell me whether these changes in the source code are enough to model a rigid wind turbine with harmonic oscillations in the surge direction or not.
To excite the turbine harmonically you could also apply a harmonic force to for example the hub. You can do that in the routine RtHs.
Thanks a lot for you reply to my post. Do you think the way I have excited the turbine harmonically is correct as well?
I am not sure I can help you with that, I have not looked into the solver part of FAST.
I do not feel confortable in changing the dof displacement, velocities and acellerations manually. They should follow from the equations of motions, under influence of forces.
If k is the surge stiffness and u is your desired harmonic displacement amplitude, I would apply a force ku (though, admittedly, applying ku would not give u because there is also the presense of the aerodynamic forces) on the platform reference point, by changing routine RtHs, though I admit that is not a trivial excerise.
Have a look at how FAST applies the following forces to the equations of motion in RtHs:
FSTipDrag(K, Aerodynamic force on blade tip brake
FSAero(K,J, Aerodynamic force on blade element
FKAero Aerodynamic force on tail fin
FTAero(J, Aerodynamic force on tower
From the similarities in the way these forces appear in the equations you might be able to predict (with eyes shut) how to add an excitiation force to the platform. The platform formulas start with
and end with
More correct would be to derive the equations on paper and then modify the code, but give it a quick try, it might work.
What you could also try is to set the surge dof to zero and give it an initial speed, similar to what you do sometimes with the generator DOF. The turbine would disappear behind the horizon at constant speed, I wonder if the simulation would run, but I think it must. Then you see the influence of turbine motion on loads. But there would be no difference with just using another windspeed. And of course you miss the dynamic wake and dynamic stall effects. For the motion to slow down and reverse (harmonic motion) you would have to apply forces that cause these decellerations. Or indeed use your method, but as said I cannot say if your implementation is correct.
Dear Gerrit ,
Many thanks for your attention.
I agree with Gerrit that the displacement, velocity, and acceleration of a DOF should not be changed directly without considering whether the equations of motion are consistent with the specified motion.
If you simply wish to prescribe a platform surge motion time-history, you can use FAST’s seismic loading module, which was developed for a related purpose. Please see my post dated Aug 29, 2012 in the forum topic found here for more information: http://forums.nrel.gov/t/software-releases/522/1.
Many thanks for your reply to my post. My aim is to find the wind thrust forces on the rotor of a rigid wind turbine. My understanding from AeroDyn and FAST was that AeroDyn reads the velocities of the blade elements from FAST, and calculate the induction factors and angle of attack accordingly. Based on this understanding, as I just wanted to have a periodic velocity in the blade elements, I just wrote them as I wished in subroutine solver, and the resulting thrust forces make sence as they are periodic with the same frequency of the turbine oscillations (under the effect of uniform wind).
I appreciate it if you could let me know whether my understanding of the way that AeroDyn and FAST are linked together is wrong, and if it is not wrong, why the approach I have used is wrong.
Your understanding of how FAST and AeroDyn are linked is correct. AeroDyn will use information from FAST regarding the blade positions and velocities (not accelerations) to calculate the aerodynamic loads. So, as long as you specify proper blade positions and velocities, then the aerodynamic loads calculated by AeroDyn should be “correct.”
However, the problem with your approach is that it does not guarantee that the specified displacement, velocity, and acceleration of a DOF are consistent with the equations of motion, from which the loads output by FAST are calculated from. So, while the applied aerodynamic loads may be fine, the reaction loads output by FAST (including rotor thrust, RotThrust) – that are derived from a combination of applied and inertia loads – may not be correct. See the following forum topic for more information on FAST’s rotor thrust output: http://forums.nrel.gov/t/fast-output-rotct-is-not-constant-for-the-same-tsr/167/1.
Many thanks for your helpful reply.
As I mentioned in the previous posts in this topic, I am modeling a rigid wind turbine oscillating harmonically in surge direction. As I am only interested in the AeroDynamic forces, I used the approach that I mentioned in my first post in the topic (changing the specific lines in subroutine solver) to oscillate the turbine rigidly. I understand that this approach will result in wrong reaction loads by FAST, but correct Aerodynamic forces by AeroDyn.
The motion in surge direction is d=acos(wt) and velocity is v=-awsin(w*t). The platform output parameters “PtfmTDxi,PtfmTVxi” have been checked and are in agreement with the input displacement and velocity. As the structure is rigid, it can be concluded that the blades have the same motion and velocity.
I have used Blade Element Momentum theory to find Aerodynamic forces. To validate the code results, I turned off all wake calculations (no induced velocity calculations, no tip,hub loss modeling, no Dynamic stall). Consequently, I expect that the angle of attack be found simply from equation 2 in AeroDyn Theory manual (Ve-ip is zero, and Ve-op is the velocity in surge direction which is prescribed as v=-awsin(wt) ). knowing the angle of attack, lift and drag coefficients can be found from related airfoil data files,consequently lift and drag forces can be found and Thrust=liftcos(phi)+dragsin(phi). Where the thrust is in fact the parameter “ForcN” which is AeroDyn’s output for each blade element in the course of time.
“ForcN” from the code is compared with the Thrust force calculated from the mentioned formulas, they are in perfect agreement as long as, in the calculations Ve-op=awsin(wt) instead of Ve-op=-awsin(w*t) . Meaning that the blade velocity that AeroDyn reads is the minus sign of the Fast’s blade velocity. I appreciate it if you could let me know the reason of this sign change , and whether I am correct in the procedure of calculations.
To investigate the issue further, I went through AeroDyn source codes. It seems the Aerodynamic forces on the blades are calculated in subroutine “elemfrc” in “Aerosubs.f90”, and the blade velocities are “VNB”. I asked the code to write the time history of “VNB” for each blade element ,exactly the same way “ForcN” is written in the Primary.elm file. Plotting the results shows the velocity is V=awsin(wt) instead of V=-awsin(wt).
Further investigation in the code shows that, “VNB” is the same as “VNElement” which is put equal to “-VES(1)” in subroutine “GetVNVT” in “Aerocalc.f90” and it seems that the minus sign used in this part is the cause of the issue (please correct me if I am wrong in tracking the code).
I appreciate it if you could let me know why the minus sign is used.
Your understanding is correct.
The reason for the minus sign in the expression
Ve-op = -v,
Ve-op = awsin(w*t),
is because Eq. (2) in the AeroDyn theory manual corresponds to Figure 2 in the AeroDyn theory manual, where it shown that positive Ve-op corresponds to an increase in the relative wind speed at the airfoil. However, your expression for v=-awsin(w*t) is an expression for the velocity in global coordinates where, in fact, positive v corresponds to a decrease in the relative wind speed at the airfoil. That is, in global coordinates, positive structural motion and positive wind speed are in the same direction, so the relative motion is
Vrel = Vwind - Vstructural,
unlike what is shown in Figure 2.
Many thanks for your reply.
I’m currently using OpenFAST to prescribe sinusoidal surge and pitch motions to a wind turbine using the ExtPtfm module. So far, I’ve been using a constant time step of 0.01 for all the simulations and the turbine undergoes 50 periods of oscillation per simulation but I think I may be using too small and too large of a time step and number of periods, respectively. The goal is to extract the average Cp, Ct, the spanwise loads and their spectra along the blades and the stall margins.
Regarding the time step)
- The range of surge frequencies and amplitudes I’m using is [0;0.1] Hz and [0;10] m which means that in highest frequency oscillation, there are 1000 time steps per period which seems to much using the rule of thumb of having at least 10 time steps per period. However, drawing from some CFD experience, the time step should be somehow a function of the surge velocity during the simulation so that the information propagates accordingly. Is this question about the timely information propagation relevant for OpenFAST since it is not CFD? If so, is there a criterion I should use to define the time step of each simulation?
Regarding using 50 periods)
- I will perform a convergence study to on the variables I want to extract to see whether 50 periods are too much but do you have an opinion on this?
What is the rotor speed?
Rules of thumb for selecting the aerodynamic and structural time steps are summarized in the following forum topic: http://forums.nrel.gov/t/fast-integrator-step-size-precision/247/2. Granted, these rules are for a flexible turbine operating in turbulence, and it sounds like you are modeling a rigid turbine, presumably in uniform inflow, so, you can probably increase the time step relative to these recommendations.
By the way: If you are truly modeling a a rigid wind turbine, are you only interested in aerodynamic quantities, not including mass/inertial effects. If so, the new driver for AeroDyn under development within OpenFAST may be of interest to you as an alternative to modeling a rigid turbine in OpenFAST with motion prescribed via ExtPtfm: github.com/OpenFAST/openfast/pull/688.
The rotor speed is constant and equal to the rated rotor speed = 7.56 rpm. The turbine is the IEA 15 MW and indeed the first set of simulations will be for a rigid turbine under uniform flow. The control is also turned off. But later, realistic wind conditions, control and elasticity will be considered.
So just to confirm some things:
The aerodynamic time step is exclusively determined by the rotor speed and has nothing to do with the relative wind velocity variation caused by the platform motion?
I can drop the structural time step recommendation because the turbine is rigid and the platform has only 1 active DOF which in practical terms is not quite active so the natural frequencies are 0?
Thank you for the link. I had the idea that it would the .fst file’s “DT” that would be defined by the system’s frequencies but I guess I was wrong. What’s the difference between the “DT”'s in the .fst and ElastoDyn input files?
Thanks for letting me know about that functionality.
Yes, we typically suggest choosing the aerodynamic time step based on the rotor speed. In your case, the rotor speed is only slightly higher than the surge frequency (7.56 rpm = 0.126 Hz), so, not too different anyway.
The natural frequencies are dictated by the modes/degrees of freedom enabled; not the time step. The DT in the OpenFAST driver (glue) code is the time step at which the modules exchange data (module-level inputs-outputs). The DT in ElastoDyn is the time step used to integrate the structural equations of motion. Often these are set identically, but sometimes it is desirable to specify the ElastoDyn time step smaller than the OpenFAST-level time step.
In case the surge frequency were higher than the rotor frequency, let’s say 0.2 Hz, then I should instead consider having 200 time steps in this frequency rather than the rotor frequency (0.126 Hz)?
Indeed, I understand that. What I meant is that since my simulation setup does not have any modes/DOF enabled apart from the platform DOF where the motion is being prescribed, then the recommendation “DT <= 1 / ( 10 * <highest full system natural frequency in Hz )” for ElastoDyn need not be considered because the system is almost infinitely rigid.
The rules of thumb were not derived for your situation of forced surge motion, so, I’m not sure. The 200 steps per azimuth came from turbulence excitation. If you aren’t simulating with turbulence, I would guess that 200 steps per azimuth would be more than sufficient even if the surge natural frequency is a bit higher than the rotor speed.
The structural time step rule of thumb would still apply to your case; simply use the highest natural frequency of the platform DOFs enabled (based on the rigid body mass/inertia and spring stiffness you’ve set).