# 4.2.1.8. Unsteady Aerodynamics

The Unsteady Aerodynamic (UA) models account for flow hysteresis, including unsteady attached
flow, trailing-edge flow separation, dynamic stall, and flow reattachment.
*Dynamic stall* refers to rapid aerodynamic changes that may bring about
or delay stall behavior [ad-Bra17]. Rapid changes in wind speed (for example, when
the blades pass through the tower shadow) cause a sudden detachment and
then reattachment of air flow along the airfoil. Such effects at the
blade surface cannot be predicted with steady state aerodynamics, but
may affect turbine operation, not only when the blades encounter tower
shadow, but also during operation in skewed flows and turbulent wind conditions. Dynamic
stall effects occur on time scales of the order of the time for the
relative wind at the blade to traverse the blade chord, approximately
\(c/\Omega r\). For large wind turbines, this might be on the order
of \(0.5\) seconds at the blade root to \(0.001\) seconds at the
blade tip. Dynamic stall can result in high transient forces as the wind
speed increases, but stall is delayed.

## 4.2.1.8.1. Theory

The different dynamic stall models implemented in AeroDyn are presented below.

### 4.2.1.8.1.1. Notations and Definitions

See Section 4.2.1.3.3 for a comprehensive description of all the inputs present in the profile input file (including some of the ones repeated below).

The airfoil section coordinate system and main variables are presented in Fig. 4.8 and further described below:

Aerodynamic Center (AC): point of the airfoil cross section where the aerodynamic forces and moment are assumed to act. Usually close to the 1/4 chord point for a regular airfoil and at the center for a circular cross section

“3/4” chord point: in the original formulation this point refers to the point on the chord axis located 3/4 chord behind the leading edge. This concept is here generalized to the point located mid-way between the aerodynamic center and the trailing edge, to account for aerodynamic center positions that differ strongly from a 1/4 chord point. The notation \(3/4\) is kept in this document.

\(\omega\): rotational speed of the airfoil section (pitching/torsional rate) positive around z.

\(\boldsymbol{v}_{ac}\): velocity vector at the aerodynamic center \(\boldsymbol{v}_{ac}=[v_{x,ac}, v_{y,ac}]\) (coordinates assumed to be expressed in the airfoil section coordinate system)

\(\boldsymbol{v}_{34}\): velocity vector at the 3/4 chord point \(\boldsymbol{v}_{34}=[v_{x,34}, v_{y,34}]\)(coordinates assumed to be expressed in the airfoil section coordinate system) The velocity is obtained from the velocity at the 1/4 chord point and the rotational speed of the section: \(\boldsymbol{v}_{34}=\boldsymbol{v}_{ac}+\omega d_{34} \hat{\boldsymbol{x}}_s\) where \(d_{34}\) is the distance between the aerodynamic center and the 3/4 chord point.

\(U_{ac}\): velocity norm at the aerodynamic center. \(U_{ac}=\lVert\boldsymbol{v}_{ac}\rVert=\sqrt{v_{x,ac}^2 + v_{y,ac}^2}\)

\(\alpha_{ac}\): angle of attack at the aerodynamic center \(\alpha_{ac}=\operatorname{atan2}(v_{x,ac},v_{y,ac})\)

\(\alpha_{34}\): angle of attack at the 3/4 chord point \(\alpha_{34}=\operatorname{atan2}(v_{x,34},v_{y,34})\)

\(\boldsymbol{x}\): the vector of states used by the continuous formulations

\(c\): airfoil chord

\(C_l^{st}, C_d^{st}, C_m^{st}\): static airfoil coefficients

\(\alpha_0\): angle of attack at zero lift, \(C_l^{st}(\alpha_0)=0\)

\(\alpha_1\): angle of attack close to positive stall.

\(\alpha_2\): angle of attack close to negative stall.

\(C_{l,\alpha}\): slope of the steady lift curve about \(\alpha_0\).

\(f^{st}_s(\alpha)\): is the steady separation function, determined from the lift curve \(C_l^{st}(\alpha)\) (see below, and e.g. [ad-HGAM04])

\(A_1\), \(A_2\), \(b_1\), \(b_2\): are four constants, characteristic of the propagation of the wake vorticity (Wagner constants)

**Time constants:**

\(T_u(t) = \frac{c}{2U_{ac}(t)} \in [0.001, 50]\): Time for the flow to go over half the airfoil section. The value is plateaued to avoid unphysical values.

\(T_{f,0}\): Dimensionless time constant associated with leading edge separation. Default is 3.

\(T_{p,0}\): Dimensionless time constant for the boundary-layer,leading edge pressure gradient. Default is 1.7

**Separation function:**

The steady separation function, \(f_s^{st}\), is defined as the separation point on a flat plate for a potential Kirchhoff flow [ad-HGAM04]:

When \(\alpha=\alpha_0\), \(f_s^{st}(\alpha_0)=1\). Away from \(\alpha_0\), the function drops progressively to \(0\). As soon as the function reaches \(0\) on both sides of \(\alpha_0\), then \(f_s^{st}\) is kept at the constant value \(0\).

**Note that for UAMod=5, a different separation function is formed.**
We define an offset for the \(C_n\) function, `cn_offset`

, where
\(C_{n,offset}=\frac{C_n\left(\alpha^{Lower}\right)+C_n\left(\alpha^{Upper}\right)}{2}\). Then, the separation function
is a value between 0 and 1, given by the following equation:

with the fully-attached \(C_n\) curve defined as \(C_n\) between \(alpha^{Lower}\) and \(alpha^{Upper}\) and linear functions outside of that range:

Note that to avoid numerical issues at the \(\pm180\) degree boundary, this function changes slope when the separation function is 0 above \(alpha^{Upper}\) and below \(alpha^{Lower}\). This allow the fully-attached linear sections to be periodic and avoid numerical issues with large magnitudes of angle of attack.

**Inviscid and fully separated lift coefficient:**

The inviscid lift coefficient is \(C_{l,\text{inv}}= C_{l,\alpha} (\alpha-\alpha_0)\). The fully separated lift coefficient may be modelled in different ways ([ad-Bra17]). In most engineering models, the slope of the fully separated lift coefficient around \(\alpha_0\) is \(C_{l,\alpha}/2\). In the Unsteady AeroDynamics sub-module, the fully separated lift coefficient is derived from the steady separation function as:

### 4.2.1.8.1.2. Beddoes-Leishman type models (UAMod=2,3)

The Beddoes-Leishman model account for attached flows and trailing edge stall [ad-LB89].

Two variants are implemented in the Unsteady Aerodynamic module. These two (compressible) models are currently described in the following reference: [ad-DH19]. The models use \(C_n\) and \(C_c\) as main physical quantities. The models use discrete states and cannot be used with linearization.

### 4.2.1.8.1.3. Beddoes-Leishman 4-states model (UAMod=4)

The 4-states (incompressible) dynamic stall model as implemented in OpenFAST is described in [ad-BJP+22] (the model differs slithgly from the original formulation from Hansen-Gaunaa-Madsen (HGM) [ad-HGAM04]).
The model is enabled using `UAMod=4`

. The model uses \(C_l\) as main physical quantity.
Linearization of the model is available.

NOTE: this model might require smaller time steps until a stiff integrator is implemented in AeroDyn-UA.

**State equation:**
The state equation of the model is:

with

**Output equation:**
The unsteady airfoil coefficients
\(C_{l,\text{dyn}}\), \(C_{d,\text{dyn}}\),
\(C_{m,\text{dyn}}\) are obtained from the states as follows:

with:

### 4.2.1.8.1.4. Beddoes-Leishman 5-states model (UAMod=5)

The 5-states (incompressible) dynamic stall model is similar to the Beddoes-Leishman 4-states model (UAMod=4), but
adds a 5th state to represent vortex generation.
It is enabled using `UAMod=5`

. The model uses \(C_n\) and \(C_c\) as main physical quantities.
Linearization of the model is available.

### 4.2.1.8.1.5. Oye model (UAMod=6)

Oye’s dynamic stall model is a one-state (continuous) model, formulated in [ad-Oye91] and described e.g. in [ad-Bra17]. The model attempts to capture trailing edge stall. Linearization of the model is available.

**State equation:**
Oye’s dynamic stall model uses one state, \(\boldsymbol{x}=[f_s]\)
where \(f_s\) is the unsteady separation function.
The state equation is a first-order differential equation:

where \(T_f=T_{f,0} T_u\) is the time constant of the flow separation and \(f_s^{st}\) is the steady state separation function described in Section 4.2.1.8.1.1. The value \(T_{f,0}\) is usually chosen around 6 (different from the default value). It is readily seen that \(f_s\) reaches the value \(f_s^{st}\) when the system is in a steady state (i.e. when \(\frac{df_s(t)}{dt}=0\)).

**Output equation:**
The unsteady lift coefficient is computed as a linear combination of the inviscid lift
coefficient, \(C_{l, \text{inv}}\), and the fully separated lift
coefficient \(C_{l,\text{fs}}\). Both of these lift coefficients are
determined from the steady lift coefficient, usually provided as
tabulated data, noted \(C_l^{st}(\alpha)\), where the superscript
\(st\) stands for “steady”.
The unsteady lift coefficient is
modelled as:

where \(\alpha_{34}\) is the instantaneous angle of attack at the 3/4 chord. \(f_s\) is seen to act as a relaxation factor between the two flow situations.

### 4.2.1.8.1.6. Boeing-Vertol model (UAMod=7)

The Boeing-Vertol is mentioned in the following paper [ad-MB11]. Details of the model were omitted in this reference, so the documentation presented here is inspired from the implementation done in the vortex code CACTUS, which was reproduced to quasi-identity in AeroDyn. Linearization is not possible with this model.

The model as presented in [ad-MB11] is an output-only model, where the dynamic angle of attack is determined using the quasi steady angle of attack and the rate of change of the angle of attack:

where \(k_1\) and \(\gamma\) are constants of the model. In practice, the implementation is different for the lift and drag coefficients, and for negative and positive stall. The model needs a discrete state to calculate the rate of change of the angle of attack and two discrete states to keep track of whether the model is activated or not.

**Airfoil constants:**

The constants \(k_1\), for positive and negative rates of angle of attack, are set to:

The extent of the transition region is computed as:

where \(\alpha_1\) and \(\alpha_2\) are the angle of attack at positive and negative stall respectively (taken as the values from the airfoil input file). The factor 0.9 is a margin to prevent the effective angle of attack to reach \(\alpha_0\) during stall.

**Intermediate variables:**

The variables \(\gamma\) for the lift and drag are computed as function of the thickness to chord ratio of the airfoil \(t_c\) and the Mach number \(M_a\) (assumed to be 0 in the current implementation):

where \(\delta = 0.06-t_c\).

**Update of discrete states (and intermediate variables):**

The rate of change of the angle of attack is computed as:

An additional state was introduced to avoid sudden jump of \(\dot{\alpha}\), by storing its value. Rates that are beyond a fraction of \(\pi \Delta t\) are replaced with the values at the previous time step. This feature is not present in the CACTUS implementation.

The dynamic angle of attack offsets (lags) for the lift and drag are computed as:

The value of \(k_1\) is taken as \(k_{1,n}\) if \(\dot{\alpha}(\alpha_{34}-\alpha_0)<0\), and taken as \(k_{1,p}\) otherwise. The lagged angle of attacks for the lift and drag are:

The distances to positive and negative stall are computed as follows. If \(\dot{\alpha}(\alpha_{34}-\alpha_0)<0\) and the dynamic stall is active:

If \(\dot{\alpha}(\alpha_{34}-\alpha_0)<0\) and the dynamic stall is not active:

If \(\dot{\alpha}(\alpha_{34}-\alpha_0)\ge0\):

The effective angle of attack for the lift coefficient is taken as the lagged angle of attack:

The effective angle of attack for the drag coefficient is obtained from the lagged angle of attack and the deltas to stall:

where \(T=2\Delta\alpha_\text{max}\) is the extent of the “transition” region.

The lift dynamic stall state is activated if \(\dot{\alpha}(\alpha_{34}-\alpha_0) \ge 0\) and if the angle of attack is above \(\alpha_1\) or below \(\alpha_2\). The state is turned off if \(\dot{\alpha}(\alpha_{34}-\alpha_0) < 0\) and the effective angle of attack is below \(\alpha_1\) and above \(\alpha_2\).

The drag dynamic stall state is activated if any of the condition below are true:

The state is turned off otherwise.

**Calculation of outputs:**
The calculation of the dynamic lift and drag coefficients is done as follows

Recalculation of intermediate variables are necessary to obtain \(\alpha_{e,L}\) and \(\alpha_{e,D}\). The moment coefficient is calculated based on values at the aerodynamic center and mid-chord (“50”):

where \(\alpha_{50}\) is computed the same way as \(\alpha_{34}\) (using the velocity at the aerodynamic center and the rotational rate of the airfoil) but using the distance from the aerodynamic center to the mid-chord (see Section 4.2.1.8.1.1).

## 4.2.1.8.2. Inputs

See Section 4.2.1.3.2.6 for a description of the inputs necessary in the AeroDyn primary file (e.g. `UAMod`

)

See Section 4.2.1.3.3 for a more comprehensive description of all the inputs present in the profile input file. Their default values are described in Section 4.2.1.8.3

See Section 4.2.1.8.1.1 for a list of notations and definitions specific to unsteady aerodynamic inputs.

An example of profile data (containing some of the unsteady aerodynamic parameters) is available here
`(here)`

.

## 4.2.1.8.3. Calculating Default Airfoil Coefficients

The default value for `cd0`

is the minimum value of the \(C_d\) curve between \(\pm20\) degrees angle of attack.
\(\alpha_{c_{d0}}\) is defined to be the angle of attack where `cd0`

occurs.

After computing `cd0`

, the \(C_n\) curve is computed by

The slope of the \(C_n\) curve is computed as follows:

\(C_{n,smooth}^{Slope}\) is a smoothed version of \(C_{n}^{Slope}\), calculated using a triweight kernel with a window of 2 degrees.

Using \(C_{n,smooth}^{Slope}\), `alphaUpper`

and `alphaLower`

are computed:

`alphaUpper`

is the smallest angle of attack value between \(\alpha_{c_{d0}}\) and 20 degrees where the \(C_{n,smooth}^{Slope}\) curve has started to decrease to 90% of its maximum slope.

`alphaLower`

is the largest angle of attack value between -20 degrees and \(\alpha_{c_{d0}}\) where the \(C_{n,smooth}^{Slope}\) curve has started to decrease to 90% of its maximum slope.

`Cn1`

is the value of \(C_n(\alpha)\) at the smallest value of \(\alpha\) where \(\alpha >= \alpha^{Upper}\) and the separation function, \(f_{st}(\alpha)\) = 0.7.

`Cn2`

is the value of \(C_n(\alpha)\) at the largest value of \(\alpha\) where \(\alpha <= \alpha^{Lower}\) and the separation function, \(f_{st}(\alpha)\) = 0.7.

`Cn_offset`

is the average value of the \(C_n\) curve at `alphaUpper`

and `alphaLower`

:

`C_nalpha`

is defined as the maximum slope of the smoothed \(C_n\) curve, \(C_{n,smooth}^{Slope}\) between \(\pm20\) degrees angle of attack.

`C_lalpha`

is defined as the maximum slope of the (un-smoothed) \(C_l\) curve, \(C_{l}^{Slope}\) between \(\pm20\) degrees angle of attack.

The default `alpha0`

is computed as the zero-crossing of a line with a slope equal to `C_lalpha`

that goes through the \(C_l\) curve at \(\alpha = \frac{\alpha^{Upper} + \alpha^{Lower}}{2}\)

`Cm0`

is the value of the \(C_m\) curve at `alpha0`

: \(C_{m,0} = C_m\left(\alpha_0\right)\). If the \(C_m\) polar values have not been included, \(C_{m,0} =0\).

`alpha1`

is the angle of attack above `alphaUpper`

where the separation function, \(f_s^{st}\) is 0.7.

`alpha2`

is the angle of attack below `alphaLower`

where the separation function, \(f_s^{st}\) is 0.7.

`Cn1`

is the value of the \(C_n\) curve at `alpha1`

.

`Cn2`

is the value of the \(C_n\) curve at `alpha2`

.

## 4.2.1.8.4. Outputs

Outputting variables of the dynamic stall models is possible, but requires
to set preprocessor variable `UA_OUTS`

and recompile the program (OpenFAST, AeroDyn Driver, or Unsteady Aero driver).
The outputs are written in output files with extension *.UA.out.
To activate these outputs with cmake, compile using `-DCMAKE_Fortran_FLAGS="-DUA_OUTS=ON"`

## 4.2.1.8.5. Driver

A driver is available to run simulations for a single airfoil, using sinusoidal variation of the angle of attack, or user defined time series of angle of attack, relative wind speed and pitch rate.

Using cmake, the driver is compiled using make unsteadyaero_driver, resulting as an executable in the aerodyn folder.

An example of driver input file is available here: `here`

.
An example of time series input is available here: `here`