Open and run in Colab (interactive) Edit on Gitlab

Wind Turbine under Yaw Misalignment

In case the wind turbine rotor is not perpendicular to the inflow, its operation and effects on the flow field will be different. In general, it is a quite complicated process. In PyWake, the effects are divided into four sub-effects that are handled invididually:

  1. Change of operation due to reduced inflow wind speed

  2. Reduced deficit, \(deficit_{normal}\), due to reduced inflow wind speed, (\(C_{T,n} \rightarrow C_{T,x}\))

  3. Reduced deficit, \(deficit_{downwind}\), due to misalignment between thrust and downwind direction

  4. Deflection of wake deficit, \(deficit_{deflected}\), due to transversal thrust component reaction

drawing

Install PyWake if needed

[1]:
# Install PyWake if needed
try:
    import py_wake
except ModuleNotFoundError:
    !pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git

First we import some basic Python elements

[2]:
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
from ipywidgets import IntSlider

We also import the site, wind turbines, and engineering wake models necessary for the simulation

[3]:
from py_wake.examples.data.hornsrev1 import V80
from py_wake.flow_map import HorizontalGrid
from py_wake.tests.test_files import tfp
from py_wake.wind_farm_models import PropagateDownwind, All2AllIterative
from py_wake.deficit_models import FugaYawDeficit
from py_wake.deflection_models import FugaDeflection
from py_wake.examples.data.hornsrev1 import Hornsrev1Site
from py_wake.flow_map import XYGrid

wt = V80()
site = Hornsrev1Site()

Effect 1 - Change of operation due to reduced inflow wind speed

The inflow perpendicular to the rotor is reduced by \(cos\theta\).

This means that the pitch and rotor speed setting of the WT will be different. In PyWake, the WT operation (pitch and rotor speed settings) is modeled in terms of the Power and CT curves. If these curves are specified in terms of the wind speed normal to the rotor, then the wind speed must be multiplied by \(cos\theta\) before looking up the power and CT.

This is one of two effects modeled by the SimpleYawModel, which is an additional model that is applied as default to most PowerCtFunctions, see here for more details.

Note, that below rated rotor speed, the CT curve is typically rather flat and reducing the wind speed has therefore limited effect on the CT value in this region. The power, on the other hand, is limited by the maximum power and therefore only affected in the region below rated wind speed (rated wind speed will increase with yaw misalignment).

Furthermore, misalignment between the thrust force (axial induction) and the inflow leads to higher wind speed at the rotor plane, which also affects the WT operation. This effect, however, is not considered in PyWake.

[4]:
#here we can see an example of the thrust and power curves for different yaw misalignment angles
u = np.arange(3,25,.1)
ax1,ax2 = plt.subplots(1,2, figsize=(15,5))[1]

#looping through different values for yaw angles
for yaw in [0,15,30]:
    ax1.plot(u,wt.ct(u*np.cos(np.deg2rad(yaw))), label='%d deg yaw misalignment'%yaw)
    ax2.plot(u,wt.power(u*np.cos(np.deg2rad(yaw))), label='%d deg yaw misalignment'%yaw)

for ax in (ax1, ax2):
    ax.legend()
    ax.set_xlabel('Wind speed [m/s]')
ax1.set_ylabel('CT')
ax2.set_ylabel('Power [W]')
[4]:
Text(0, 0.5, 'Power [W]')
../_images/notebooks_YawMisalignment_8_1.png

Effect 2 - Reduced deficit due to reduced inflow wind speed (\(C_{T,n} \rightarrow C_{T,x}\))

The wake deficit is caused by a reaction to the thrust force, which slows down the inflow.

The thrust force normal to the rotor plane is

\begin{align} T_n & =\frac{1}{2}\rho C_{T,n} A (U cos\theta)^2 \end{align}

In non-aligned inflow, the thrust force that slows down the flow, i.e. in the mean wind direction is

\begin{align} T_x & =\frac{1}{2}\rho C_ {T,x} \left(A cos\theta \right) U^2 \end{align}

From these two equations we can find the relationship between the thrust coefficient in the rotor-normal direction, \(C_{T,n}\), and the thrust coefficient in the down-wind direction, \(C_{T,x}\).

\begin{align} T_x & = T_n cos\theta \\ \frac{1}{2} \rho C_{T,x} \left(A cos\theta \right) U^2 & = \frac{1}{2}\rho C_{T,n} A (U cos\theta)^2 cos\theta\\ C_ {T,x} & = C_{T,n} \cos^2\theta \end{align}

This is the second effect modeled by the SimpleYawModel, which is an additional model that is applied as default to most PowerCtFunctions and is shown below.

SimpleYawModel

The SimpleYawModel is an additional model, which as default is applied by

  • PowerCtFunction

  • PowerCtTabular

  • PowerCtNDTabular

  • PowerCtXr

  • CubePowerSimpleCt

It handles effects 1. and 2. by:

  1. Compute/look up \(C_{T,n}\) based on wind speed \(U\cos\theta\)

  2. Returns \(C_{T,x} = C_{T,n} \cos^2\theta\)

Effect 3 - Reduced deficit due to misalignment between thrust and downwind direction

Most engineering models calculate the deficit normal to the rotor plane, while the deficit impacting downstream WT is the deficit in the downstream direction.

In case of yaw misalignment, the deficit is therefore mapped to downstream direction by

\(deficit_{downwind} = deficit_{normal} cos \theta\)

To model this effect, PyWake uses the FugaYawDeficit model.

This model requires two sets of look-up tabels, namely UL and UT. These tables describes the normalized deficit in the downwind, U direction caused by a unit forcing in the longitudinal, L, and transveral, T, directions.

[5]:
fuga = FugaYawDeficit(tfp + 'fuga/2MW/Z0=0.00001000Zi=00400Zeta0=0.00E+00.nc')
D = 80
R = D/2
x0,x1,y = 1,10,2 #upstream, downstream and crosswind regions to load [D]
xi0 = int(512-x0/(fuga.dx/D))
xi1 = int(512+x1/(fuga.dx/D))+1
yi = int(y/(fuga.dy/D))+1
UL, UT = [-fuga.mirror(v,anti_symmetric=a) for v, a in zip(fuga.load_luts()[:2,0,:yi, xi0:xi1], (False,True))]
/tmp/ipykernel_7910/1431188295.py:5: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  xi0 = int(512-x0/(fuga.dx/D))
/tmp/ipykernel_7910/1431188295.py:6: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  xi1 = int(512+x1/(fuga.dx/D))+1
/tmp/ipykernel_7910/1431188295.py:7: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  yi = int(y/(fuga.dy/D))+1
[6]:
axes = plt.subplots(2,1, figsize=(15,12))[1]
X, Y = np.meshgrid(fuga.x[xi0:xi1], fuga.mirror(fuga.y[:yi], anti_symmetric=True))
wtL = np.array([[0,0],[-D/2,D/2]])
for ax, v, l in zip(axes.flatten(), [UL, UT], ['UL','UT']):
    c = ax.contourf(X, Y, v)
    plt.colorbar(c, ax=ax, label=r'Downwind deficit $\left[\frac{1}{U^2 CT}\right]$')
    s = (slice(None,None,2),slice(None,None,5))
    ax.quiver(X[s], Y[s], -v[s]*(l[0]=='U'), v[s]*(l[0]=='V'))
    ax.plot([0,0],[-D/2,D/2],'r')
    if l[1]=='L':
        ax.arrow(-10,0,-20,0,color='r', width=10,head_length=30)
        ax.arrow(60,0,-20,0,color='r', width=10,head_length=30)
    else:
        ax.arrow(-30,25,0,-20,color='r', width=10,head_length=30)
        ax.arrow(30,25,0,-20,color='r', width=10,head_length=30)

    ax.axis('equal')
    ax.set_xlim([-D,10*D])
    ax.set_ylim([-D,D])
    ax.set_title(l)
../_images/notebooks_YawMisalignment_15_0.png

As Fuga is a linear model, the deficit in the downstream direction can simply be calculated by the linear superposition of UL and UT:

\(deficit = (UL \cdot cos\theta + UT \cdot sin\theta) \cdot U \cdot CT\)

You can move the slider below to see the deficit fields for different yaw-misalignment angles

[7]:
def plot(yaw):
    plt.figure(figsize=(18,5))
    theta = np.deg2rad(yaw)
    co, si = np.cos(theta), np.sin(theta)
    deficit = co**2*(UL*co + UT*si)
    c = plt.contourf(X, Y, deficit,20)#, levels=np.arange(-.4,.1,.05))
    plt.colorbar(c, label=r'Downwind deficit $\left[\frac{1}{U^2 CT}\right]$')
    s = (slice(None,None,2),slice(None,None,5))
    plt.quiver(X[s], Y[s], -deficit[s], deficit[s]*0, width=.002, scale=2)
    plt.plot([si*R, -si*R],[-co*R,co*R],'w',label='Rotor',lw=3)
    plt.arrow(-0,0,-100*co,0,color='r', width=5,head_length=10,length_includes_head=True)
    plt.arrow(-0,-0,0,-100*si,color='r', width=5,head_length=10,length_includes_head=True, zorder=32)
    plt.axis('equal')
    plt.xlim([-D,10*D])
    plt.ylim([-D/2,D/2])
    ax.set_title(l)
_ = interact(plot, yaw=IntSlider(min=-90, max=90, step=1, value=20, continuous_update=False))

Effect 4 - Wake deflection due to traversal thrust component reaction

Wake deflection is modeled by a DeflectionModel, which modifies the downwind, horizontal crosswind and vertical distance between wind turbines.

The FugaDeflection is capable of modeling the wake deflection if VL and VT is known.

These tables describes the normalized deficit in the crosswind, V direction caused by a unit forcing in the longitudinal, L, and transveral, T, directions.

[8]:
fugaDeflection = FugaDeflection(tfp + 'fuga/2MW/Z0=0.00001000Zi=00400Zeta0=0.00E+00.nc')
D = 80
R = D/2
x0,x1,y = 1,10,2 #upstream, downstream and crosswind regions to load [D]
xi0 = int(512-x0/(fugaDeflection.dx/D))
xi1 = int(512+x1/(fugaDeflection.dx/D))+1
yi = int(y/(fugaDeflection.dy/D))+1

# VL is antisymmetric, VT is symmetric
VL, VT = [fugaDeflection.mirror(v,anti_symmetric=a) for v, a in zip(fugaDeflection.load_luts()[2:4,0,:yi, xi0:xi1], (True, False))]
/tmp/ipykernel_7910/2912431644.py:5: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  xi0 = int(512-x0/(fugaDeflection.dx/D))
/tmp/ipykernel_7910/2912431644.py:6: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  xi1 = int(512+x1/(fugaDeflection.dx/D))+1
/tmp/ipykernel_7910/2912431644.py:7: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  yi = int(y/(fugaDeflection.dy/D))+1
[9]:
axes = plt.subplots(2,1, figsize=(15,12))[1]
X, Y = np.meshgrid(fugaDeflection.x[xi0:xi1], fugaDeflection.mirror(fugaDeflection.y[:yi], anti_symmetric=True))
for ax, v, l in zip(axes.flatten(), [-VL, VT], ['VL','VT']):
    c = ax.contourf(X, Y, v)
    plt.colorbar(c, ax=ax, label=r'$WS_{Crosswind} \left[\frac{1}{U^2 CT}\right]$')
    s = (slice(None,None,10),slice(None,None,1))
    ax.quiver(X[s], Y[s], 0, v[s], scale=(.3,3)[l[1]=='T'])
    ax.plot([0,0],[-D/2,D/2],'r')
    if l[1]=='L':
        ax.arrow(-10,0,-20,0,color='r', width=10,head_length=30)
        ax.arrow(60,0,-20,0,color='r', width=10,head_length=30)
    else:
        ax.arrow(-30,25,0,-20,color='r', width=10,head_length=30)
        ax.arrow(30,25,0,-20,color='r', width=10,head_length=30)

    ax.axis('equal')
    ax.set_xlim([-D,10*D])
    ax.set_ylim([-D,D])
    ax.set_title(l)
../_images/notebooks_YawMisalignment_21_0.png

As Fuga is a linear model, the effects in the crosswind direction can simply be calculated by linear superposition of VL and VT:

\(WS_{crosswind} = (VL \cdot cos\theta + VT \cdot sin\theta) \cdot U \cdot CT\)

You can move the slider below to see the crosswind wind speed field for different yaw-misalignment angles. Furthermore, the integrated deflection is shown with solid lines

[10]:
def plot(yaw):
    plt.figure(figsize=(18,5))
    theta = np.deg2rad(yaw)
    co, si = np.cos(theta), np.sin(theta)
    V = co**2*(-VL*co + VT*si)
    X, Y = np.meshgrid(fugaDeflection.x[xi0:xi1], fugaDeflection.mirror(fugaDeflection.y[:yi], anti_symmetric=True))
    c = plt.contourf(X, Y, V)#, levels=np.arange(-.09,.1,.01))
    plt.colorbar(c, label=r'$WS_{crosswind} \left[\frac{1}{U CT}\right]$')
    s = (slice(None,None,6),slice(None,None,1))
    plt.quiver(X[s], Y[s], V[s]*0, V[s], width=.002, scale=2)

    # plot deflection lines
    fL, fT = fugaDeflection.fLT.V.T[:,256-yi:255+yi, xi0:xi1]
    lambda2p = co**2 * (fL * co - fT * si)
    y = fugaDeflection.mirror(fugaDeflection.y[:yi], anti_symmetric=True)

    lambda2 = np.array([np.interp(y, y + l2p, l2p) for l2p in lambda2p.T])
    Yp = Y + lambda2.T
    plt.plot(X[yi, :], Y[yi, :], 'grey', lw=3)
    for x, y, yp in zip(X[1::4], Y[1::4], Yp[1::4]):
        plt.plot(x, y, 'grey', lw=1, zorder=-32)
        plt.plot(x, yp, 'k', lw=1)
    plt.plot(X[yi, :], Yp[yi, :], 'k', lw=3)

    plt.plot([si*R, -si*R],[-co*R,co*R],'w',label='Rotor',lw=3)
    plt.arrow(-0,0,-50*co,0,color='r', width=5,head_length=10,length_includes_head=True, zorder=32)
    plt.arrow(-0,-0,0,-50*si,color='r', width=5,head_length=10,length_includes_head=True, zorder=32)
    plt.axis('equal')
    plt.xlim([-D,10*D])
    plt.ylim([-D,D])
    ax.set_title(l)
_ = interact(plot, yaw=IntSlider(min=-90, max=90, step=1, value=20, continuous_update=False))

Note, the plot above shows the wind speed in cross-wind direction (up), i.e. yellow colors corresponds to a flow towards the top of the figure, while blue corresponds to a flow towards the bottom.

Finally, all effects are put together. You can move the slider below to see the effects on the deficit behind a WT.

[11]:
path = tfp + 'fuga/2MW/Z0=0.00001000Zi=00400Zeta0=0.00E+00.nc'
fugaDeficit = FugaYawDeficit(path)
wfm = All2AllIterative(site,wt,fugaDeficit,blockage_deficitModel=fugaDeficit,
                       deflectionModel=FugaDeflection(path)
                       )
def plot(yaw):
    plt.figure(figsize=(18,5))
    fm = wfm([0],[0],ws=10,wd=270, yaw=[[[yaw]]]).flow_map(XYGrid(x=np.arange(-100,800,10), y=np.arange(-250,250,10)))
    X,Y = np.meshgrid(fm.x, fm.y)
    c1 = plt.contourf(X,Y,fm.WS_eff.squeeze(),np.arange(6.5,10.1,0.1), cmap='Blues_r')
    c2 = plt.contourf(X,Y,fm.WS_eff.squeeze(),np.arange(10,10.2,.005), cmap='Reds')
    plt.colorbar(c2, label='Wind speed (deficit regions) [m/s]')
    plt.colorbar(c1, label='Wind speed (speed-up regions) [m/s]')
    wt.plot([0],[0],wd=270,yaw=yaw)
    max_deficit_line = fm.min_WS_eff(x=np.arange(-100,800,10))
    max_deficit_line.plot(color='k', label='Max deficit line')
    plt.axhline(0, label='Center line')
    plt.legend()


_ = interact(plot, yaw=IntSlider(min=-90, max=90, step=1, value=20, continuous_update=False))

Note, the scaling of wake deficit(blue) and speedup(red) does not match