This module only:

Aerodynamic loads

1. Lift forces, drag forces, and torsional moments on elements

Airfoil sections which are not a part of a wind turbine may be subjected to aerodynamic loads based on interpolation from the static lift, drag, and moment curves, or Morison-type quadratic drag loads.

1.1. Morison-type aerodynamic drag loads (CHTYPE=MORI)

The drag load calculation for Morison-type aerodynamic drag is based on the relative wind velocity in the local system:

\[ \boldsymbol{u_r}=\boldsymbol{u_{wind}}-\boldsymbol{\dot {v}}=u_{rx}\boldsymbol{i_1}+u_{ry}\boldsymbol{i_2}+u_{rz}\boldsymbol{i_3}\]

where $\mathrm {\boldsymbol{u_r}}$ is the relative velocity at the element center, $\mathrm {\boldsymbol{u_{wind}}}$ is the undisturbed incoming wind velocity at the element center, and $\mathrm {\boldsymbol{\dot {v}}}$ is the structural velocity at the element center. In the local system, the x-axis is along the element.

The magnitude of the wind velocity normal to the elements is used in the calculation of the drag load:

\[ u_{rn}=\sqrt{u_{ry}^2+u_{rz}^2}\]

Using the dimensional drag coefficients ($\mathrm {C_{Dx}}$, $\mathrm {C_{Dy}}$, $\mathrm {C_{Dz}}$), the aerodynamic load per unit length is calculated as:

\[ f=C_{Dx}u_{rx}|u_{rx}|\boldsymbol{i_1}+C_{Dy}u_{ry}u_{rn}\boldsymbol{i_2}+C_{Dz}u_{rz}u_{rn}\boldsymbol{i_3}.\]

The force per unit length is multiplied by the dry length (following the same water level integration as the wave loads) and applied in a lumped manner.

For an axisymmetric cross section, it is possible to give non-dimensional aerodynamic coefficients ($\mathrm {C_{dx}}$, $\mathrm {C_{dy}}$, $\mathrm {C_{dz}}$) as input rather than dimensional coefficients:

\[ C_{Dx}=\frac{1}{2}\rho _{air}S_{2D}C_{dx}\]

\[ C_{Dy}=\frac{1}{2}\rho _{air}B_yC_{dx}\]

\[ C_{Dz}=\frac{1}{2}\rho _{air}B_zC_{dx}\]

where $\mathrm {S_{2D}}$ is the cross-sectional dry surface, $\mathrm {B_y}$ is the projected area per unit length in the local $\mathrm {y}$ direction, and $\mathrm {B_z}$ is the projected area per unit length in the local $\mathrm {z}$ direction. For a circular cross section with diameter $\mathrm {D}$,

\[ S_{2D}=\pi D\]

\[ B_y=B_z=D.\]

1.2. Lift forces, drag forces, and torsional moments on elements (CHTYPE=AIRF)

For airfoil sections which are not a part of a wind turbine, aerodynamic loads are based on interpolation from the static lift, drag, and moment curves.The wind velocity and structural velocity are computed at the center of each airfoil element based on the undisturbed wind speed and structural velocity at the nodes. All computations are carried out in the local airfoil coordinate system. The lift force per unit length is denoted $\mathrm {F_L}$, the drag force per unit length is $\mathrm {F_D}$, and the aerodynamic moment per unit length is $\mathrm {M}$.

Figure 1. Local airfoil coordinate system

The relative velocity in the airfoil coordinate system is computed as:

\[\boldsymbol{V_{r}}=\boldsymbol{V_{\mathrm {wind}}}-\boldsymbol{V_{\mathrm {struc}}}\]

The spanwise ($\mathrm {V_{r_{z}}}$) component of the velocity is assumed to be small, and an error is returned if the spanwise component is large compared to the velocity magnitude. The angle of attack ($\mathrm {\alpha }$) is computed based on the local $\mathrm {V_{r_{x}}}$ and $\mathrm {V_{r_{y}}}$ components of the velocity. The Reynolds number (Re) is computed as

\[\text{Re}=\frac{\rho _{\mathrm {air}}c\sqrt{(V_{r_{x}}^2+V_{r_{y}}^2)}}{\nu_{\mathrm {air}}}\]

where $\mathrm {\rho _{air}}$ is the air density, $\mathrm {c}$ is the chord length, and $\mathrm {\nu_{air}}$ is the viscosity of air.

Lift ($\mathrm {C_L}$), drag ($\mathrm {C_D}$), and moment ($\mathrm {C_M}$) coefficients are then obtained by interpolation in the airfoil tables based on Re and $\mathrm {\alpha }$. The forces and moment (about the airfoil coordinate system) per unit length are:

\[F_L=\frac{1}{2}C_L\rho _{\mathrm {air}}c(V_{r_{x}}^2+V_{r_{y}}^2)\]

\[F_D=\frac{1}{2}C_D\rho _{\mathrm {air}}c(V_{r_{x}}^2+V_{r_{y}}^2)\]

\[M=\frac{1}{2}C_M\rho _{\mathrm {air}}c^2(V_{r_{x}}^2+V_{r_{y}}^2)\]

Then, $\mathrm {F_L}$, $\mathrm {F_D}$, and $\mathrm {M}$ are multiplied by the dry length of the element, transformed into the appropriate coordinate system, and are applied to the nodes in a lumped manner. The dry length of the element is determined consistent with the wave load integration definition.

1.3. Loads due to Vortex-Induced Vibrations

See also Loads due to Vortex-Induced Vibrations for the time domain VIV load formulation, which also applies to hyrodynamic VIV loading.

2. Blade element/momentum (BEM) theory for wind turbines

The aerodynamic load model used for wind turbines in RIFLEX and SIMO is based on the blade element momentum (BEM) theory. This theory combines momentum theory and blade element theory. The general principle of the theory is that forces developed locally at the airfoil, based upon empirical lift and drag coefficients, are balanced with the change in momentum of the air flowing through the rotor disk.

2.1. Outline of BEM calculation

The balance of airfoil forces with the change in momentum of the air is the basic step in the calculation. In the implemented calculation, the dynamic wake method is used. Only one solution sequence is performed for each time step, and time is used a relaxation parameter. In a quasi-static analysis, an iterative calculation would be performed for each time step.

An outline of the basic force balance is given below. Each wind turbine blade is divided into blade elements. For each blade element, aerodynamic airfoil coefficients (lift, drag and moment) are defined.

The calculation begins with a suitable guess for the induced velocity at each blade element. Induced velocity is defined as the difference between the velocity of air at the rotor plane and the velocity of air far upstream of the rotor plane.
The relative air velocity with respect to each blade element is calculated.
The angle-of-attack with respect to the airfoil is calculated.
Using airfoil coefficient data, the aerodynamic lift force, drag force, and moment are calculated using the Øye dynamic stall model.
Airfoil forces are rotated into the rotor plane coordinate system.
Annulus-average forces and velocities are computed, for use in the momentum balance. Rotor plane-average remote and induced velocities are computed, for use in computing the Prandtl factor.
The Prandtl factor is calculated. This factor modifies the basic momentum balance to account for a finite number of blades.
Induced velocity (the average over each annulus) is updated based upon momentum balance, using annulus-average forces.
The calculation is so far quasi-static. The next step is the Øye dynamic wake model, which acts as a filter for the induced velocity. The inputs are the quasi static induced velocity and the induced velocity from the last time step.
In cases where the rotor is yawed with respect to the remote wind, the wake skew angle is computed, and the variation of induced velocity with blade azimuth angle is computed.

2.2. Correction factors

A number of correction factors are applied in the BEM method.

2.2.1. Glauert correction

The relation between thrust and induced velocity from BEM theory is not valid for large induction factors. Instead, an empirical correction as described by Burton et al (2001, p67) is used in case of high thrust loading:

\[a=\frac{(C_T/F-C_{T1})}{C_{T2}-C_{T1}}(a_2-a_1)+a_1\]

where

$\mathrm {a_2=1.0}$, $\mathrm {C_{T2}=1.82}$, $\mathrm {a_1=1.0-0.5\sqrt{C_{T2}}}$, $\mathrm {C_{T1}=4a_1(1-a_1)}$, and $\mathrm {F}$ is the Prandtl factor Equation (15). For $\mathrm {F=1}$, the thrust curve used in the calculations is compared to the momentum theory result shown in the figure Figure 2.

The momentum theory is applied for $\mathrm {C_{T} \leq C_{T1}}$ and the Glauert correction is applied for , $\mathrm {C_{T} > C_{T1}}$.

Figure 2. Glauert correction

2.2.2. Prandtl factor

To account for the tip and hub loss on a blade due to a finite number of blades, the Prandtl factor is implemented. The Prandtl factor $\mathrm {F}$ as a function of the radial location $\mathrm {r}$ is computed as:

\[F=\frac{2}{\pi }\cos^{-1}{(e^{-f})}\]

where

\[f=\frac{B}{2}\frac{R-r}{2r\sin{\phi }}\]

$\mathrm {B}$ is the number of blades, $\mathrm {R}$ is the outer radius of the rotor, and $\mathrm {\phi }$ is the angle that the trailing vorticies make with the rotorplane. $\mathrm {\phi }$ is assumed to be the same as the angle of incoming flow, neglecting tangential induced velocity. The default program behavior is to correct $\mathrm {\phi }$ for yawed inflow.

2.2.3. Dynamic wake

The dynamic wake effect is the time lag in induced velocities due to the shedding and downstream convection of vorticity. Dynamic wake effects are most pronounced for heavily loaded rotors, corresponding to high induction factors (low wind speeds). This effect can be modeled by the Stig Øye dynamic inflow model, which acts as a filter for induced velocities, as shown below (Snel and Schepers, 1995). $\mathrm {\boldsymbol{W_{qs}}}$ is the quasi static induced velocity vector, and $\mathrm {\boldsymbol{W}}$ is the new induced velocity vector. $\mathrm {\tau_1}$ and $\mathrm {\tau_2}$ are time constants.

\[\boldsymbol{W}+\tau_2\frac{d\boldsymbol{W}}{dt}=\boldsymbol{W}_{\mathrm {int}}\]

\[\boldsymbol{W}_{\mathrm {int}}+\tau_1\frac{d\boldsymbol{W}_{\mathrm {int}}}{dt}=\boldsymbol{W}_{qs}+0.6\tau_1\frac{d\boldsymbol{W}_{qs}}{dt}\]

where

\[\tau_1= \frac{1.1}{1-1.3a}\frac{R}{V_o}\]

\[\tau_2= (0.39-0.26 \left(\frac{r}{R}\right)^2)\tau_1\]

2.2.4. Dynamic stall

Dynamic stall refers to time-dependent variations in the lift and drag coefficients of an airfoil due to changes in the angle of attack. The Stig Øye model is implemented and gives unsteady lift by filtering the trailing edge separation point with an empirical time constant (see also Hansen, 2008).

The degree of stall ($\mathrm {f_s}$) affects the lift as

\[C_L=f_sC_{L,\mathrm {inv}}(\alpha )+(1-f_s)C_{L,fs}(\alpha )\]

where $C_{L,\mathrm {inv}}$ is the lift coefficient without separation, and $\mathrm {C_{L,fs}}$ is the fully separated lift coefficient. A static value $\mathrm {f_{s}^{st}}$ can be found to reproduce the static airfoil data. The degree of stall is assumed to chase the static value:

\[\frac{df_s}{dt}=\frac{f_s^{st}-f_s}{\tau}\]

where $\mathrm {\tau}$ is a time constant. By integration:

\[f_s(t+\Delta t)=f_s^{st}+(f_s(t)-f_s^{st})e^{-\Delta t/\tau}.\]

For flow that is not fully stalled, the lift coefficient is computed as:

\[C_L=\frac{1}{4}\frac{dC_L}{d\alpha }(\alpha -\alpha _0)(1+\sqrt{1-|f_s|})^2\]

where $\mathrm {\frac{dC_L}{d\alpha }}$ is computed at the full-stall limit (depending on the sign of $\mathrm {\alpha }$), and $\mathrm {\alpha _0}$ is the angle of attack where the lift is zero.

For fully stalled flow:

\[C_L=C_{L,qs}(1+\sqrt{1-|f_s|})^2\]

where $\mathrm {C_{L,qs}}$ is the quasi-static lift.

In order to initialize the dynamic stall calculation, the program identifies several important points in the lift curve: the angle of attack where there is zero lift ($\mathrm {\alpha _0}$), the maximum slope of the lift curve in the linear region for both positive and negative lift values ($[dC_L/d\alpha ]_{\mathrm {max}}(1)$ and $[dC_L/d\alpha ]_{\mathrm {max}}(2)$), and the angle of attack corresponding to full separation for both positive and negative lift values ($\mathrm {\alpha _{fs}(1)}$ and ($\mathrm {\alpha _{fs}(2)}$)) For typical airfoils, these points are relatively easy to identify. For more complex airfoils, the user may input these values directly.

Figure 3. Illustration of dynamic stall initialization parameters.

2.2.5. Skewed inflow

Skewed inflow may occur due to rotor tilt or a yaw angle between the rotor and the oncoming wind. Glauert developed the basic formulation for correcting the induction factor due to skewed inflow: $\mathrm {a}$ is the induction factor, $\mathrm {\chi}$ is the wake skew angle, $\mathrm {K}$ is a user defined parameter (default=1), and $\mathrm {\Psi}$ is the azimuth angle that is zero at the most downwind position of the rotor.

\[a_{\mathrm {skew}}=a[1+K\tan{(\frac{\chi}{2})}(\frac{r}{R})\cos{(\Psi)}]\]

The skew angle used here is based on Hansen (2008), and is the angle between the wind velocity in the wake and the rotor’s rotational axis. The skew angle is assumed constant with radius and is computed at a radial position close to r/R=0.7.

2.2.6. Upwind tower influence

The tower influence (shadow) effect is implemented using potential flow theory. Since the incoming wind has to travel around the tower, the tower has an effect on the local inflow, even for an upwind turbine. The velocity deficit due to the tower is based on the 2D potential solution for constant flow around a circle, with the possibility to include corrections for the tower drag and the Bak coefficient. The coordinate system for the tower influence calculation is shown in figure Figure 4.

Figure 4. Tower influence coordinate system

For an upwind wind turbine, the non-dimensional influence at a location $\mathrm {(x,y,z)}$ is evaluated as:

\[x_w=\frac{2x}{D_{\mathrm {tow}}}\]

\[y_w=\frac{2y}{D_{\mathrm {tow}}}\]

\[x_{\mathrm {infl}}=[1-\frac{(x_w+b)^2-y_w^2}{((x_w+b)^2+y_w^2)^2}+\frac{C_D}{2\pi }\frac{x_w+b}{(x_w+b)^2+{y_w}^2}]\]

\[y_{\mathrm {infl}}=-2[\frac{(x_w+b)y_w}{((x_w+b)^2+y_w^2)^2}+\frac{C_D}{2\pi }\frac{y_w}{(x_w+b)^2+{y_w}^2}]\]

where $D_{\mathrm {tow}}$ is the local tower diameter for the given $\mathrm {z}$ level, $\mathrm {b}$ is the Bak coefficient, and $\mathrm {C_D}$ is the local tower drag coefficient for the given $\mathrm {z}$ level. If the effect of tower drag is not considered, the Bak factor is set to zero and potential flow solution is applied.

In practice, the horizontal wind direction and speed are computed at the location $\mathrm {(x,y,z)}$. This speed is multiplied by the velocity factors $x_{\mathrm {infl}}$ and $y_{\mathrm {infl}}$, and the modified speed is then transposed back to the initial wind direction (Moriarty and Hansen, 2005).

Note that the $\mathrm {x}$-coordinate here is upwind and the velocity is set to zero for any point inside the tower. The tower is assumed to be approximately vertical and no effect on the vertical wind speed is included.

When the blade segment is above the tower height, the velocity is smoothened to match the free wind. The resulting velocity deficit factors is shown in the figures below. The blade is pointing upwards when the azimuth angle is 0 deg, and is closest to the tower at 180 deg. In these figures the distance from the hub to the center of the segment is 90 m, the tower top is located 6 m below the hub and the tower diameter is ranging from 9 m (azimuth=180 deg) to 6.54 m (tower top). There is no cone on the blade line and no tilt of the shaft.

Figure 5. The velocity factors due to the tower shadow for all azimuth angles of the blade

Figure 6. The velocity factors due to the tower shadow close to the tower

2.2.7. Downwind tower influence

In order to model downwind wind turbines, a different method for computing the velocity deficit due to the presence of the tower is needed. The most common models, which are implemented in AeroDyn and Bladed, for example, are based on work by Powles. Powles’ model, as implemented in RIFLEX, assumes a cosine-squared shape for the tower shadow.

The cosine-squared model for the wake factor $u_{\mathrm {wake}}$ is: $ u_\{}=

$ where $\mathrm {d=\sqrt{x^2+y^2}}$. Together with the potential flow model, the local wind velocity is found as: $ U_\{}=(u-u_\{})U_$ $ V_\{}=(v-u_\{})U_$

where $\mathrm {U_\infty}$ is the total instantaneous horizontal velocity, and $\mathrm {u=1}$ when the cosine-squared model is in effect. In Powles’ model, the flow is assumed to not reverse: for $u-u_{\mathrm {wake}}<0$, we assume $U_{\mathrm {local}}=0$. The computed $U_{\mathrm {local}}$ and $V_{\mathrm {local}}$ are then transformed back to the local global reference frame. Thus, the wake is assumed to align with the instantaneous horizontal wind vector.

The resulting velocity deficit for a downwind turbine is shown in figure Figure 7.

Figure 7. Velocity deficit due to the tower, downwind wind turbine

2.3. BEM Limitations

The BEM aerodynamic load module has the following limitations: 1. The code does not truly account for blade cone or large deflections 2. The vortical wake structure must be preserved. For this to be the case, the ratio of blade tip speed to wind speed must not be too low, and the majority of the blade must not be stalled. 3. The cone angle should be small. If the rotor is coned: for purposes of load calculation, define the aerodynamic analysis as if the cone angle were zero; for a lower bound on power output, define the aerodynamic analysis with the actual cone angle.

3. Loads due to Vortex-Induced Vibrations (VIV)

3.1. Load Formulation

The time domain VIV load per unit length is given by:

\[\begin{aligned} F &= (C_A + 1)\rho\frac{\pi D^2}{4}\dot{u}_n - C_A\rho\frac{\pi D^2}{4}\ddot{x}_n + \frac{1}{2}\rho C_D D v_n|v_n| \\ &+ \frac{1}{2}\rho D C_{v,CF}|v_n|(j_3 \times v_n)\cos\Phi_{v,CF} \\ &+ \frac{1}{2}\rho D C_{v,IL}|v_n|v_n\cos\Phi_{v,IL} \\ &+ \frac{1}{2}\rho D C_{hh}|v_n|(j_3 \times v_n)\cos(\Phi_{v,CF} + \Phi_{v,IL}) \end{aligned}\]

Where:

The first line is equivalent to the ordinary Morison load
The second line is the cross-flow vortex shedding force term
The third line is the in-line vortex shedding force term
The fourth line is the higher harmonic vortex shedding force term

3.2. Synchronization Model

The cross-flow vortex shedding force is given by:

\[F_{v,CF} = \frac{1}{2}\rho D C_{v,CF}|v_n|(j_3 \times v_n)\cos\Phi_{v,CF}\]

where $\Phi_{v,CF}$ is the time varying instantaneous phase of the cross-flow vortex shedding force.

For a fixed force frequency this phase would be:

\[\Phi_{v,CF} = 2\pi f_{v,CF}\, t\]

Instead, the frequency is adjusted to synchronize force and response.

$\Phi_{v_n}$: Instantaneous phase of the relative velocity
$\Phi_{v,CF}$: Instantaneous phase of the CF force
$\theta$: Instantaneous phase difference between relative velocity and force, $\theta = \Phi_{v_n} - \Phi_{v,CF}$

Adjust the phase of the cross-flow load $\Phi_{v,CF}$ to reduce the phase difference.

The IL synchronization model is similar to the CF one, as illustrated in Figure 8.

Figure 8. Synchronization curve illustrating the lock-in (frequency) bandwidth and the phase difference between the external forcing and the structural response under synchronization (phase-locking).

The non-dimensional frequency is defined as:

\[\hat{f} = \frac{f_{osc} \cdot D}{U}\]

The phase is updated at each time step as:

\[\Phi_{v,CF,upd} = \Phi_{v,CF} + \frac{2\pi|v_n|}{D}\,\hat{f}\,\Delta t\]

Figure 9. Illustration of the synchronization concept. The phase difference between vortex shedding force and relative structure velocity is reduced by increasing or decreasing instantaneous frequency of the force at every time step.

3.3. References

Thorsen, M. J. (2016): Time Domain Analysis of Vortex-Induced Vibrations. PhD Thesis, Norwegian University of Science and Technology, Trondheim, Norway.
Ulveseter, J. V. (2018): Advances in semi-empirical time domain modelling of vortex-induced vibrations. PhD Thesis, Norwegian University of Science and Technology.
Kim, S. (2021): Non-Linear Time Domain Analysis of Deepwater Riser Vortex-Induced Vibrations. PhD Thesis, Norwegian University of Science and Technology, Trondheim, Norway.