Turbine theory: Aerodynamic load

Aerodynamic loads

The aerodynamic loads are calculated individually on each blade mid-segment frame, at the aerodynamic centre. The loads are due to the inflow $\vec{w}$ at incidence angle $\alpha$:

Most conventional turbines have a clockwise rotation sense and the above figure applies. If the turbine has an anticlockwise rotation sense, then inflow is determined as shown in the below figure:

In these diagrams, we are looking from the blade tip towards the root so the $\vec{z}$ axis points out of the page towards the viewer.

The incidence angle or angle of attack, $\alpha$, is always in the range $-180\degree \le\alpha\le +180\degree$. The first and last angles in your table of coefficient data must be -180° and +180° respectively, and (clearly, since these two angles represent the same direction) the two sets of coefficient values for these angles must coincide.

The lift force, for conventional turbines with clockwise rotation sense, is in the direction $\vec{w}{\times}\vec{z}$. For anticlockwise rotation sense, it is in the direction $-\vec{w}{\times}\vec{z}$. It is calculated as \begin{equation} f\urm{L} = \frac12 \rho\, A\, \C{l}(\alpha)\, \lvert\vec{w}\rvert^2 \end{equation} The drag force, in the $\vec{w}$ direction, is \begin{equation} f\urm{D} = \frac12 \rho\, A\, \C{d}(\alpha)\, |\vec{w}|^2 \end{equation} The pitching moment, for conventional turbines with clockwise rotation sense, is about the $z$-axis. For anticlockwise rotation sense, it is about the $-z$-axis. It is calculated as \begin{equation} m\urm{z} = \frac12 \rho\, A\, c\, \C{m}(\alpha)\, |\vec{w}|^2 \end{equation} where

$\rho$ is the air density

$c$ is the chord

$A$ is the element area, i.e. the product of the chord and the unstretched segment length

$\C{l}(\alpha)$, $\C{d}(\alpha)$, $\C{m}(\alpha)$ are the lift, drag and moment coefficients respectively, obtained for angle of attack $\alpha$ by linear interpolation of the wing type data . If a UA model is active, the modified unsteady coefficients are used.

Blade element momentum theory

When induction is not included, the inflow $\vec{w}$ corresponds simply to the disturbed wind field at, and relative to, the segment's aerodynamic centre, $\vec{v}$, and we set \begin{align} w\urm{x} &= v\urm{x} \\ w\urm{y} &= v\urm{y} \\ w\urm{z} &= 0 \end{align} In this case, the inflow used does not account for the interaction with the rotor itself.

To account for induction effects from the rotor, we can apply blade element momentum (BEM) theory which combines momentum conservation principles with classical blade sectional theory. In this way we can calculate axial and tangential induction factors, $a$ and $a'$, relating the quasi-steady velocity induced by the rotor, $\vec{v}\urm{q}$, to the instantaneous disturbed relative velocity, $\vec{v}$, as \begin{align} v\urm{q,x} &= - a v\urm{x} \label{v_q} \\ v\urm{q,y} &= + a' v\urm{y} \\ v\urm{q,z} &= 0 \end{align} Where quasi-steady refers to induction being calculated from the instantaneous flow conditions. If dynamic inflow is not enabled, the inflow, $\vec{w}$, is then simply the sum of the relative velocity and the quasi-steady induced velocity \begin{equation} \vec{w} = \vec{v} + \vec{v}\urm{q} \label{w_x} \end{equation} When dynamic inflow is enabled, the Øye model is used to approximate the dynamic induced velocity, $\vec{v}\urm{d}$, from the quasi-steady values. The inflow is then calculated as \begin{equation} \vec{w} = \vec{v} + \vec{v}\urm{d} \label{wd_x} \end{equation} In all cases (whether induction is included or not), the velocities are measured with respect to the nominal rotor plane frame, and the relative flow in the $z$ direction is taken as zero: i.e. we assume that there is no flow along the blade axis.

BEM is strictly speaking a 2D theory, framed in the rotor plane: for a simple 2D conventional clockwise rotor, subject to flow normal to the rotor-plane, the terms axial and tangential refer unambiguously to flow downwind and opposite to the rotational velocity respectively.

In OrcaFlex, however, the blades support pre-cone and deflection, so the blades need not be co-planar and each individual blade is not even constrained to a single plane. To generalise the BEM method to allow for this 3D geometry we define, for each blade segment, a nominal rotor plane frame of reference, and it is in this plane that we apply the BEM method. The figure above is drawn in the nominal rotor plane.

Nominal rotor plane

For a conventional clockwise rotor, the nominal rotor plane frame is best expressed in terms of the orientation of the segment's instantaneous geometry frame. The geometry frame, relative to the instantaneous fitting frame associated with that segment's blade, can be defined as three successive Euler angle rotations: $R\urm{x}$, $R\urm{y}$, and $R\urm{z}$. To recover the orientation of the geometry frame from the fitting frame, one would start from the fitting frame and rotate it by $R\urm{x}$ about its $x$-axis first, followed by $R\urm{y}$ about the new $y$-axis, and then finally by $R\urm{z}$ about the new $z$-axis. This final $R\urm{z}$ rotation is equal and opposite to the effective twist; for an undeformed blade the effective twist is the cumulative aerodynamic twist, plus any blade pitch, at the mid-segment arc length.

The orientation of the nominal rotor plane frame, then, is found from the fitting frame by applying only the first two of these Euler angle rotations and neglecting the final effective twist rotation. The result is a frame that has a $z$-axis along the blade segment and $x$ and $y$ axes in the cross-section of the aerofoil. The $x$-axis and $y$-axis define the axial and tangential directions respectively of the BEM calculation. The nominal rotor plane is essentially the geometry frame without aerodynamic twist or blade pitch.

For a rotor with an anticlockwise rotation sense, the nominal rotor plane is determined by following the same procedure. As shown in the above figure, the nominal rotor plane frame's $y$-axis points in the direction of nominal rotation, and (in the case of no twist) is opposite to the geometry frame's $y$-axis.

Calculating the induction factors

We calculate the induction factors $a$ and $a'$ by reducing the problem to one dimension, the local inflow angle $\phi$, as proposed by Ning et al and generalised by NREL in the implementation of AeroDyn. $\phi$ is the angle between $\vec{w}$ and the nominal rotor plane $y$-axis; below, we derive the residual function $\mathcal{R}(\phi)$ which we solve to obtain the induction factors.

We begin with the angle of attack, $\alpha$, calculated as \begin{equation} \alpha = \phi - \textrm{effective twist} \end{equation} For a conventional clockwise rotor, the lift and drag coefficients can be re-expressed in the nominal rotor plane $x$ and $y$ directions as \begin{align} C\urm{x}(\alpha) &= \C{l}(\alpha) \cos{\phi} + \C{d}(\alpha) \sin{\phi} \\ C\urm{y}(\alpha) &= \C{d}(\alpha) \cos{\phi} - \C{l}(\alpha) \sin{\phi} \end{align} For an anticlockwise rotor, \begin{equation} C\urm{y}(\alpha) = \C{l}(\alpha) \sin{\phi} - \C{d}(\alpha) \cos{\phi} \end{equation} Now, for convenience, we calculate two non-dimensional parameters $\kappa$ and $\kappa'$ \begin{equation} \begin{aligned} \kappa &= \frac{\sigma\ C\urm{x}(\alpha)}{4 F {\sin^2}\phi} \\ \kappa' &= -\frac{\sigma\ C\urm{y}(\alpha)}{4 F \sin{\phi} \cos{\phi}} \end{aligned} \label{kappa} \end{equation} where $\sigma$ is the local blade solidity, defined as \begin{equation} \sigma = \frac{n\urm{B}\ c}{2\pi\ r} \end{equation} and

$n\urm{B}$ is the blade count

$c$ is the local chord

$r$ is instantaneous radius, defined to be the displacement of the mid-segment frame from the turbine origin projected onto the turbine $xy$-plane

The factor $F$ in equation (\ref{kappa}) accounts for Prandtl tip and hub loss and is calculated as follows \begin{align} f\urm{tip} &= \frac12 n\urm{B} \frac{x\urm{b}-x}{(x+r\urm{hub}) \lvert\sin{\phi}\rvert} \\ F\urm{tip} &= \frac{2}{\pi} \cos^{-1}\{\exp(-f\urm{tip})\} \\ f\urm{hub} &= \frac12 n\urm{B} \frac{x}{r\urm{hub} \lvert\sin{\phi}\rvert} \\ F\urm{hub} &= \frac{2}{\pi} \cos^{-1}\{\exp(-f\urm{hub})\} \\ F &= F\urm{tip}\ F\urm{hub} \label{Flossfactor} \end{align} where

$x\urm{b}$ is the total unstretched blade length

$x$ is the unstretched mid-segment arc length

If tip loss is not enabled, substitute $F\urm{tip}{=}1$ in (\ref{Flossfactor}); similarly, if hub loss is not enabled, substitute $F\urm{hub}{=}1$.

The induction factors are then calculated. First the solution region is identified. If $\phi\ v\urm{x} \gt 0$, the turbine is deemed to be operating in the momentum/empirical region and the axial induction factor is given by \begin{equation} a = \begin{cases} \dfrac{\kappa}{1+\kappa} & \text{if $\kappa \leq \frac{2}{3}$} \\ \dfrac{\gamma_1-\sqrt{\gamma_2}}{\gamma_3} & \text{otherwise} \end{cases} \end{equation} where $\gamma_i$ are empirical coefficients given by Glauert/Buhl as \begin{equation} \begin{aligned} \gamma_1 &= 2 F\ \kappa - \left( \frac{10}{9} - F \right) \\ \gamma_2 &= 2 F\ \kappa - F \left( \frac{4}{3} - F \right) \\ \gamma_3 &= 2 F\ \kappa - \left( \frac{25}{9} - 2 F \right) \end{aligned} \end{equation} If the turbine is not operating in the momentum/empirical region, it is deemed to be in the propeller brake region and the axial induction factor is given by \begin{equation} a = \frac{\kappa}{\kappa-1} \end{equation} In all regions the tangential induction factor is given by \begin{equation} a' = \begin{cases} \dfrac{\kappa'}{1-\kappa'} & \text{if $v\urm{x} \gt 0$} \\ \dfrac{-\kappa'}{1+\kappa'} & \text{otherwise} \end{cases} \end{equation} Finally, we have the residual function \begin{equation} \mathcal{R}(\phi) = \frac{\sin \phi}{1-a} - \frac{v\urm{x}}{v\urm{y}} (1-\kappa') \cos \phi \end{equation} We solve the equation $\mathcal{R}(\phi)=0$ iteratively for $\phi$, using a bracketed root finding algorithm subject to the user-specified tolerance, and the resulting induction factors are then used to determine the value of $\vec{w}$ used in the aerodynamic loading calculation.

If an induction factor is not included, its value is simply set to zero. If it is included, we impose a limit on its magnitude: the minimum value we apply in each case is -1; the maximum axial induction factor applied is 1.5; and the maximum tangential induction factor applied is 1. If neither is included, the BEM calculation is omitted completely.

The induction factors and aerodynamic loads are calculated for each blade segment, for every iteration of the solver. If, in the final converged dynamic equilibrium, $\kappa{\lt}1$ and the turbine is in the propeller brake region, or $\kappa{\lt}{-}1$ and it is in the momentum/empirical region, the solution is invalid. In this case, or if the BEM method otherwise failed to converge, the induction factors are set to zero and (provided that the simulation time exceeds the monitoring start time) a warning will be raised.

Dynamic inflow

Calculating the induced velocities from the instantaneous flow is known as the equilibrium wake assumption, i.e the wake is assumed to react instantaneously to a change in the conditions. The induced velocities determined under the equilibrium wake assumption are known as the quasi-steady values, $\vec{v}\urm{q}$. By default, these quasi-steady values are used when calculating the inflow velocity, as in equation (\ref{w_x}).

In reality, the wake takes time to respond to a change in the flow conditions. To capture the physics associated with the delayed wake response, the Øye dynamic inflow model can optionally be enabled. In this model, adapted from AeroDyn, a series of first order filters are applied to the quasi-steady induced velocities to approximate the dynamic values $\vec{v}\urm{d}$. Application of the first filter results in intermediate induced velocities, $\vec{v}\urm{int}$. This filter is written in the form of a state equation, as \begin{equation} \dot{\vec{v}}\urm{int} = \frac{1}{\tau} \vec{v}\urm{q} + 0.6 \dot{\vec{v}}\urm{q} - \frac{1}{\tau} \vec{v}\urm{int} \end{equation} where $\tau$ is either a user-specified constant or calculated as \begin{equation} \tau = \frac{1.1}{1-1.3\bar{a}}\frac{r\urm{x}}{\bar{v}\urm{a}} \end{equation} where

$r\urm{x}$ is the maximum instantaneous radius. This is the largest nodal displacement from the turbine origin (projected onto the turbine $xy$-plane), evaluated over all blade nodes

$\bar{v}\urm{a}$ is the mean disturbed flow velocity, normal to the turbine $xy$-plane, evaluated over all blade elements

$\bar{a}$ is the mean axial induction, evaluated over all blade elements

To determine the dynamic induced velocities, the intermediate values are input into a second filter, which is dependent on the blade element's radial position, $r$. This filter is written in the form of a state equation, as \begin{equation} \dot{\vec{v}}\urm{d} = \frac{1}{\tau'} \left( \vec{v}\urm{int} - \vec{v}\urm{d} \right) \end{equation} where \begin{equation} \tau' = \left[0.39-0.26 \left( \frac{r}{r\urm{x}} \right)^2 \right] \tau \end{equation} The dynamic induced velocity is then used to calculate the inflow, using equation (\ref{wd_x}). When dynamic inflow is enabled, the induction factors are recalculated to relate the dynamic induced velocity to the instantaneous disturbed relative velocity, as \begin{align} a = - \frac{v\urm{d,x}}{v\urm{x}} \qquad a' = + \frac{v\urm{d,y}}{v\urm{y}} \end{align} These effective values are used for reporting the induction factors. The axial value is also input into the skewed wake correction calculation (discussed below), if that feature is enabled.

Skewed wake correction

The Pitt and Peters skewed wake correction, if enabled, is applied to the converged axial induction factor $a$: it is applied retrospectively, and is not coupled to the BEM iterative solution process. The axial induction factor, or its effective value in the case of dynamic inflow, is modified \begin{equation} a\urm{skew} = a \left(1 + F_{s} \frac{r}{r\urm{x}} \tan{\frac{\chi}{2}} \sin{\psi} \right) \end{equation} where

$\psi$ is a segment azimuth angle, defined below

$\chi$ is the wake skew angle

$F_{s}$ is the user-specified skewed wake factor. If it is '~', it is calculated as

\begin{equation} F_{s} = \frac{15 \pi}{32} \end{equation}

$\psi$ and $\chi$ are both functions of the mean disturbed flow direction, $\vec{\bar{v}}$. This is the mean direction of the disturbed relative flow, expressed with respect to the global frame, evaluated over all blade segments.

$\psi$ is found by first projecting the segment $z$-axis direction, $\vec{s\urm{z}}$, and $\vec{\bar{v}}$ onto the turbine $xy$-plane to get $\vec{s}\urm{xy}$ and $\vec{\bar{v}}\urm{xy}$ respectively. $\psi$ is then the angle between $\vec{s}\urm{xy}$ and $\vec{\bar{v}}\urm{xy} \times \vec{t\urm{z}}$, where $\vec{t\urm{z}}$ is the turbine frame $z$-axis direction.

The wake skew angle is calculated as \begin{equation} \chi = (0.6a + 1)\chi_0 \end{equation} where $\chi_0$ is the skew angle, defined to be the angle between $\vec{t\urm{z}}$ and $\vec{\bar{v}}$.

The magnitude of $\chi$ is limited to $\frac{\pi}{2}$ before being used to calculate the skewed wake correction.

Finally, having calculated the corrected induction factor $a\urm{skew}$, we apply this factor in place of the original factor $a$ in equation (\ref{v_q}).

Induction scaling

In OrcaFlex, the tip speed ratio is defined as \begin{equation} TSR = \frac{\omega r\urm{x}}{\bar{v}\urm{a}} \end{equation} where

$r\urm{x}$ is the maximum instantaneous radius. This is the largest nodal displacement from the turbine origin (projected onto the turbine $xy$-plane), evaluated over all blade nodes

$\bar{v}\urm{a}$ is the mean disturbed flow velocity, normal to the turbine $xy$-plane, evaluated over all blade elements

$\omega$ is the main shaft angular velocity

If the tip speed ratio is too low, the BEM could be operating in a regime for which the underlying theory is invalid. To account for this, the induction factors can optionally be smoothly scaled from zero to their full value as a function of TSR. An induction weight is calculated \begin{equation} w\urm{i} = \begin{cases} 0 & \text{if $|TSR| \leq TSR\urm{l}$} \\ 1 & \text{if $|TSR| \geq TSR\urm{u}$} \\ \frac12 \left(1 - \cos \pi\frac{|TSR| - TSR\urm{l}}{TSR\urm{u} - TSR\urm{l}} \right) & \text{otherwise} \end{cases} \end{equation} where

$TSR\urm{l}$ is the user-specified lower TSR bound

$TSR\urm{u}$ is the user-specified upper TSR bound

The induction weight, $w\urm{i}$, then is used to scale the induction factors $a$, $a'$, and $a\urm{skew}$ before they are used to calculate the inflow, equation (\ref{w_x}). The TSR bounds can be set to zero to supress induction scaling.

If the induction weight is less than 1, a warning will be raised (provided that the simulation time exceeds the monitoring start time).

Unsteady aerodynamics

In a steady aerodynamic treatment, the static (i.e. steady) aerofoil coefficients, as specified in the wing type coefficient table, are used directly. In reality, due to skewed/turbulent inflow, controller action, turbine motion etc., the angle of attack and inflow velocity are not steady, and unsteady flow effects can become important. Introducing rapid change, delay, and hysteresis of the aerodynamic load.

Two methods: González and Minnema Pierce, based on the semi-empirical Beddoes-Leishman model, are included to account for unsteady aerodynamic effects, including: unsteady attached flow, trailing-edge flow separation, dynamic stall, and flow reattachment. Both models act to modify the aerofoil coefficients, before they are used to determine the aerodynamic loads, as a function of the inflow conditions \begin{equation} \C{l,UA}, \C{d,UA}, \C{m,UA} = f \left(\vec{w}, \alpha, c, \C{l}, \C{d}, \C{m}, \textrm{UA data} \right) \\ \end{equation} where

$\vec{w}$ is the inflow velocity, equation (\ref{w_x}), accounting for dynamic inflow, skewed wake correction, and induction scaling

$\alpha$ is the inflow incidence angle

$c$ is the user-specified air speed of sound

$\C{l}, \C{d}, \C{m}$ are the steady aerofoil coefficients

UA data, is the user-specified wing type unsteady aerodynamics data used to parameterise the UA models

The UA models are implemented based on those developed by NREL for AeroDyn and are documented in Damiani and Hayman. The modifications are made after the induction factors are determined, i.e. the steady aerofoil coefficients are used in the BEM calculation.

For large angles of attack, or low inflow velocities, the UA models might become invalid. If the magnitude of the angle of attack is within 5 degrees of the wing type UA cutout angle, or the magnitude of the inflow velocity has fallen below 1 m/s, the steady and unsteady aerofoil coefficients will be blended \begin{equation} \begin{aligned} \C{l,blend} &= \gamma\urm{UA} \C{l,UA} + \left(1 - \gamma\urm{UA} \right) \C{l} \\ \C{d,blend} &= \gamma\urm{UA} \C{d,UA} + \left(1 - \gamma\urm{UA} \right) \C{d} \\ \C{m,blend} &= \gamma\urm{UA} \C{m,UA} + \left(1 - \gamma\urm{UA} \right) \C{m} \end{aligned} \end{equation} where $\gamma\urm{UA}$ is the UA weight, calculated via \begin{equation} \gamma\urm{UA} = \gamma\urm{\alpha} \gamma\urm{w} \end{equation} where \begin{equation} \gamma\urm{\alpha} = \begin{cases} 0 & \text{if $|\alpha| \geq \alpha\urm{c}$} \\ 1 & \text{if $|\alpha| \leq \alpha\urm{c} - \delta\urm{\alpha}$} \\ 1 - \frac12 \left(1 - \cos \pi\frac{|\alpha| - (\alpha\urm{c} - \delta\urm{\alpha})}{\delta\urm{\alpha}} \right) & \text{otherwise} \end{cases} \end{equation} and \begin{equation} \gamma\urm{w} = \begin{cases} 0 & \text{if $|\vec{w}| \leq w\urm{min}$} \\ 1 & \text{if $|\vec{w}| \geq w\urm{c}$} \\ \frac12 \left(1 - \cos \pi\frac{|\vec{w}| - w\urm{min}}{w\urm{c} - w\urm{min}} \right) & \text{otherwise} \end{cases} \end{equation} here

$\vec{w}$ is the inflow velocity, equation (\ref{w_x}), accounting for dynamic inflow and skewed wake correction

$\alpha$ is the inflow incidence angle

$\alpha\urm{c}$ is the user-specified UA cutout angle

$\delta\urm{\alpha}$, taken to be 5 degrees, is the size of the angle range over which the coefficients are blended

$w\urm{c}$, taken to be 1 m/s, is the inflow velocity below which the coefficient are blended

$w\urm{min}$ is taken to be 0.01 m/s

If blending, or the Mach number exceeds 0.3, OrcaFlex will display a warning, provided the simulation time exceeds the monitoring start time. The Mach number is calculated \begin{equation} M = \frac{\vec{w}}{c} \end{equation} where

$\vec{w}$ is the inflow velocity, equation (\ref{w_x}), accounting for dynamic inflow and skewed wake correction

$c$ is the user-specified air speed of sound