Aerodyn dynamic stall model

Good afternoon,
I have never used the dynamic stall model of Aerodyn and I am also still not so used with unsteady airfoils aerodynamics.
In the airfoils data file if InclUAdata is set to true, there are 32 coefficients to be given (example is shown below)


True InclUAdata ! Is unsteady aerodynamics data included in this table? If TRUE, then include 30 UA coefficients below this line
!..
-4.2 alpha0 ! 0-lift angle of attack, depends on airfoil.
8 alpha1 ! Angle of attack at f=0.7, (approximately the stall angle) for AOA>alpha0. (deg)
-8 alpha2 ! Angle of attack at f=0.7, (approximately the stall angle) for AOA<alpha0. (deg)
1 eta_e ! Recovery factor in the range [0.85 - 0.95] used only for UAMOD=1, it is set to 1 in the code when flookup=True. (-)
6.2047 C_nalpha ! Slope of the 2D normal force coefficient curve. (1/rad)
3 T_f0 ! Initial value of the time constant associated with Df in the expression of Df and fā€™ā€˜. [default = 3]
6 T_V0 ! Initial value of the time constant associated with the vortex lift decay process; it is used in the expression of Cvn. It depends on Re,M, and airfoil class. [default = 6]
1.7 T_p ! Boundary-layer,leading edge pressure gradient time constant in the expression of Dp. It should be tuned based on airfoil experimental data. [default = 1.7]
11 T_VL ! Initial value of the time constant associated with the vortex advection process; it represents the non-dimensional time in semi-chords, needed for a vortex to travel from LE to trailing edge (TE); it is used in the expression of Cvn. It depends on Re, M (weakly), and airfoil. [valid range = 6 - 13, default = 11]
0.14 b1 ! Constant in the expression of phi_alpha^c and phi_q^c. This value is relatively insensitive for thin airfoils, but may be different for turbine airfoils. [from experimental results, defaults to 0.14]
0.53 b2 ! Constant in the expression of phi_alpha^c and phi_q^c. This value is relatively insensitive for thin airfoils, but may be different for turbine airfoils. [from experimental results, defaults to 0.53]
5 b5 ! Constant in the expression of Kā€™ā€˜ā€™_q,Cm_q^nc, and k_m,q. [from experimental results, defaults to 5]
0.3 A1 ! Constant in the expression of phi_alpha^c and phi_q^c. This value is relatively insensitive for thin airfoils, but may be different for turbine airfoils. [from experimental results, defaults to 0.3]
0.7 A2 ! Constant in the expression of phi_alpha^c and phi_q^c. This value is relatively insensitive for thin airfoils, but may be different for turbine airfoils. [from experimental results, defaults to 0.7]
1 A5 ! Constant in the expression of Kā€™ā€˜ā€™_q,Cm_q^nc, and k_m,q. [from experimental results, defaults to 1]
0 S1 ! Constant in the f curve best-fit for alpha0<=AOA<=alpha1; by definition it depends on the airfoil. [ignored if UAMod<>1]
0 S2 ! Constant in the f curve best-fit for AOA> alpha1; by definition it depends on the airfoil. [ignored if UAMod<>1]
0 S3 ! Constant in the f curve best-fit for alpha2<=AOA< alpha0; by definition it depends on the airfoil. [ignored if UAMod<>1]
0 S4 ! Constant in the f curve best-fit for AOA< alpha2; by definition it depends on the airfoil. [ignored if UAMod<>1]
1.4144 Cn1 ! Critical value of C0n at leading edge separation. It should be extracted from airfoil data at a given Mach and Reynolds number. It can be calculated from the static value of Cn at either the break in the pitching moment or the loss of chord force at the onset of stall. It is close to the condition of maximum lift of the airfoil at low Mach numbers.
-0.5324 Cn2 ! As Cn1 for negative AOAs.
0.19 St_sh ! Strouhalā€™s shedding frequency constant. [default = 0.19]
0.006 Cd0 ! 2D drag coefficient value at 0-lift.
-0.121 Cm0 ! 2D pitching moment coefficient about 1/4-chord location, at 0-lift, positive if nose up. [If the aerodynamics coefficients table does not include a column for Cm, this needs to be set to 0.0]
0 k0 ! Constant in the \hat(x)_cp curve best-fit; = (\hat(x)_AC-0.25). [ignored if UAMod<>1]
0 k1 ! Constant in the \hat(x)_cp curve best-fit. [ignored if UAMod<>1]
0 k2 ! Constant in the \hat(x)_cp curve best-fit. [ignored if UAMod<>1]
0 k3 ! Constant in the \hat(x)_cp curve best-fit. [ignored if UAMod<>1]
0 k1_hat ! Constant in the expression of Cc due to leading edge vortex effects. [ignored if UAMod<>1]
0.2 x_cp_bar ! Constant in the expression of \hat(x)_cp^v. [ignored if UAMod<>1, default = 0.2]
ā€œDEFAULTā€ UACutout ! Angle of attack above which unsteady aerodynamics are disabled (deg). [Specifying the string ā€œDefaultā€ sets UACutout to 45 degrees]
10 filtCutOff ! Cut-off frequency (-3 dB corner frequency) for low-pass filtering the AoA input to UA, as well as the 1st and 2nd derivatives (Hz) [default = 20]


So my question will be very basic as a beginner. I am wondering how the user determines these coefficients for a given blade. It seems that some of them can be calculated directly from the static polars. Perhaps some other have to be tuned properly by an experimented user.
Many thanks for any tips and any help.
Best regards.
Florence Haudin.

Hi Florence,

Setting up the unsteady aerodynamics parameters can be a bit challenging. Itā€™s important to note that from the 32 parameters requested, only 8 parameters must be calculated by the user. The other parameters can be left as ā€œdefaultā€ or 0. The 8 parameters of interest are: alpha0, alpha1, alpha2, C_nalpha, Cn1, Cn2, Cd0, and Cm0.

These 8 parameters are calculated based on the airfoil static polars.

I have modified one MATLAB function originally written by Srinivas Guntur. The new MATLAB file uses a neutral input format (the original one was using AeroDyn 14 files as input).

Attached you can find the MATLAB function (AeroDyn_UA_coefficients.m) and one example (Example_airfoil.xlsx) from OC6 phase III project. Note that the main MATLAB function has one dependency (intersections.m). Once downloaded and stored, you can simply run the MATLAB function with the example provided as:

blade_polars = xlsread('Example_airfoil.xlsx'); UA_coefficients = AeroDyn_UA_coefficients(blade_polars);

The 8 coefficients are stored in the function output (UA_coefficients) and also printed on the screen automatically:

I hope that helps!

Roger
Unsteady_Aerodynamics_coeffs.zip (15.9 KB)

Info for not MATLAB users:
The MATLAB function posted is fully compatible with GNU Octave. So, in case you donā€™t have a MATLAB license, you can still run it.

Roger

Dear Roger,

I want to confirm if at present Matlab-toolbox is still fully compatible with GUN Octave?

If so, it would be very good for users with no license.

Update:

I try to use GNU Octave to run runlinear.m:

clear all; close all; clc;
restoredefaultpath;
addpath(genpath('D:/Program/openfast/matlab-toolbox-main'));
data=ReadFASTLinear('5MW_Land_DLL_WindSteady.1.lin');
subB=data.B(8:13,343:348);

But it aborted with error:

error: structure has no member 'B'

I checked the date: B, C, D are missing compared to Matlab data structure

Regards,
Ran

Dear @Ran.Tu,

Iā€™m not sure of the compatibility of the MATLAB toolbox with GNU Octave. But please note that Roger was not referring to the MATLAB toolbox in his post; instead, Roger was referring to the MATLAB script that he shared on Oct 08 '21, to compute unsteady airfoil aerodynamic coefficients.

Best regards,

Dear @Ran.Tu,

Jason is right. I was talking about the GNU Octave compatibility for the functions that I posted (i.e., AeroDyn_UA_coefficients.m and intersections.m). I donā€™t think the Matlab toolbox for OpenFAST is fully compatible with GNU Octave.

Kind regards,

Roger

Thank you @Roger.Bergua and @Jason.Jonkman , Indeed, it seems that Matlab toolbox is not fully compatible with GNU Octave.

@Roger.Bergua Thanks for this calculator. However, I got question on effect of Straudal number on fatigue loads. If the blade sections got distiguish vortex sheddings resulting in change in amplitude and frequency of the unsteady loads how would you incorporate that? Could you please suggest?

Iā€™m not sure I understand your question, @Ramesh.Kumar. If you are interested in predicting vortex shedding and how it evolves in time, you could use OLAF within the AeroDyn module in OpenFAST. OLAF is a free vortex wake approach that will let you examine the shedding frequency as well as the time-dependent loading during dynamic stall conditions.

@Roger.Bergua Thanks for reply.

Apologies, Let me explain again.

I want to use airfoil section which has modified polar due to its optimisation or aero upgrade on it which results into breaking down the large vortex sheddings into smaller eddies at high frequency and which results in affecting the unsteady flow and unsteady loads over airfoil (ignore the effect of flex on blade). Could we incorporate this? Is there a way to force to have small vortices or smaller ranges on the time dependent loads?

I am still talking to have airfoil operating well below static stall angle and affecting the sheeding due to airfoil shape. Am I making sense?

I seeā€¦ Iā€™m not aware of modifications that could be performed. But I will let others comment. It may be that you need to run higher fidelity models to capture that behavior.

@Roger.Bergua I tried to run the OLAF model of IEA15MW but I get error in reading OLAF.dat file. Could you please help?

Dear @Ramesh.Kumar,

It seems that the OLAF input file that you are using is not compatible with the OpenFAST version that you are trying to run. What is the OpenFAST version that you are running?

Dear @Ramesh.Kumar,

I agree with @Roger.Bergua. I would just add that as with any input file processing error, I would suggest enabling the Echo option (in this case within AeroDyn) to debug errors associated with input file formatting.

Best regards,

Dear @Roger.Bergua,

I am using 3.5.0, I managed to run it with update format of wake extent and discretisation section, but I didnot find useful information yet, I will deep dive and see if i can make progress.

Best regards

Dear @Jason.Jonkman

Is it possible to have different wake/trailing edge vortex shedding or to account unsteady aerodynamic effectā€¦ I can see the unsteady section of the polar has wagner constants.

A1, A2, b1,b2: are four constants, characteristic of the propagation of the wake vorticity (Wagner constants). How does these coefficient changes with change is vorcity at TE (e.g. aifoil with serrations)?

Best regards
Ramesh

Dear @Ramesh.Kumar,

Iā€™m not sure I fully understand your question, but you can certainly change these airfoil data inputs differently between distinct airfoil data files in AeroDyn. To see how these parameters are used within the unsteady airfoil aerodynamics submodel of AeroDyn, I would review the theory basis here: 4.2.1.8. Unsteady Aerodynamics ā€” OpenFAST v3.5.2 documentation

Best regards,

Dear Jason,

Thanks for the link but I struggle to find which account transient loads on the airfoil due to change in vortex shedding.

The transient load on airfoil due to vortex induced vibration would be different for the vortex shedding of baseline airfoil as compared to the shedding of the airfoil with serrations or any aerodynamic addonsā€™.

So instead of having low frequency large vortex sheet in baseline, the upgrade airfoil has high frequency small eddies.

Where or how we can account that in unsteady aerodynamic input?

Best regards
Ramesh

Dear @Ramesh.Kumar,

I asked @Emmanuel.Branlard and here is his response:

Right now, I wouldnā€™t be able to answer the question directly. I would think one could tune these coefficients; they are linked to the time response after a step change in angle of attack. If they have data for a step change of a given airfoil, they could fit the time response to the Wagner function.

Best regards,

Dear @Jason.Jonkman ,

Thanks for checking. I have looked in the source code Unsteadyaero.f90 ā†’ ComputeKelvinChain.

I could not find the calculation where wagner function is implemented? I tried changing the straudal number for all section on the blade, but do not see any effect of unsteady load history.

Best regards
Ramesh