Engineering Wind Farm Models Object

In this example, all the engineering models that PyWake uses for the wake definition and calculation are explained. All wind farm models take a Site and a WindTurbines object as input.

Install PyWake if needed

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

Model overview

The engineering wind farms models in PyWake are composed of one of two wind farms models in combination with a wake deficit model, a superposition model and optionally a blockage deficit and a turbulence model. In addition, a rotor average model , a deflection model and a ground model are also available.

A brief definition of each can be found in the Overview section.

Predefined WindFarmModels

The deficit models comprise:

Name	WindFarmModel	Wake DeficitModel	Blockage DeficitModel	SuperpositionModel
NOJ	PropagateDownwind	NOJDeficit		SquaredSum
Fuga	PropagateDownwind	FugaDeficit		LinearSum
FugaBlockage	All2AllIterative	FugaDeficit	FugaDeficit	LinearSum
BastankhahGaussian	PropagateDownwind	BastankhahGaussianDeficit		SquaredSum
IEA37SimpleBastankhahGaussian	PropagateDownwind	IEA37SimpleBastankhahGaussianDeficit		SquaredSum

• Default rotor-average model: RotorCenter

• Default turbulence model: None

First we import basic Python and PyWake elements

import numpy as np
import matplotlib.pyplot as plt
import os

# import and setup site and windTurbines
import py_wake
from py_wake.examples.data.hornsrev1 import V80, Hornsrev1Site

site = Hornsrev1Site()
windTurbines = V80()
wt_x, wt_y = site.initial_position.T

Engineering WindFarmModel base classes

PropagateDownwind

The PropagateDownwind wind farm model is very fast as it only performs a minimum of deficit calculations. It iterates over all turbines in downstream order. In each iteration it calculates the effective wind speed at the current wind turbine as the free stream wind speed minus the sum of the deficit from upstream sources. Based on this effective wind speed, it calculates the deficit caused by the current turbine on all downstream destinations. Note, that this procedure neglects upstream blockage effects.

for wt in wind turbines in downstream order:
    ws_eff[wt] = ws[wt] - superposition(deficit[from_upstream_src,to_wt])
    ct = windTurbines.ct(ws_eff[wt])
    deficit[from_wt,to_downstream_dst] = wakeDeficitModel(ct, distances[from_wt,to_downstream_dst], ...)

The proceedure is illustrated in the animation below:

• Iteration 1: WT0 sees the free wind (10m/s). Its deficit on WT1 and WT2 is calculated.

• Iteration 2: WT1 sees the free wind minus the deficit from WT0. Its deficit on WT2 is calculated and the effective wind speed at WT2 is updated.

All2AllIterative

The All2AllIterative wind farm model is slower but is capable of handling blockage effects. It iterates until the effective wind speed converges (i.e. less than or equal to the maximum number of turbines that affect each other in the wind farm. The converge tolerance is an input parameter. In each iteration it sums up the deficit from all wind turbine sources and calculates the deficit caused by the all wind turbines turbine on all wind turbines.

while ws_eff not converged:
    ws_eff[all] = ws[all] - superposition(deficit[from_all,to_all])
    ct[all] = windTurbines.ct(ws_eff[all])
    deficit[from_all,to_all] = wakeDeficitModel(ct[all], distances[from_all,to_all], ...)

The proceedure is illustrated in the animation below:

• Iteration 1: All three WT see the free wind (10m/s) and their CT values and resulting deficits are therefore equal.

• Iteration 2: The local effective wind speeds are updated taking into account the wake and blockage effects of the other WT. Based on these wind speeds the CT and deficits are recalculated.

• Iteration 3: Repeat after which the flow field has converged.

Wake deficit models

The wake deficit models compute the wake deficit caused by a single wind turbine. In PyWake, there are several wake deficit models, which include:

• NOJDeficit

• TurboNOJDeficit

• FugaDeficit

• BastankhahGaussianDeficit

• IEA37SimpleBastankhahGaussianDeficit

• NiayifarGaussianDeficit

• ZongGaussianDeficit

• CarbajofuertesGaussianDeficit

• TurboGaussianDeficit

• GCLDeficit

Note: The implementation of the deficit models is highly vectorized and therefore suffixes are used to indicate the dimension of variables. The suffixes used in this context are:

• i: source wind turbines

• j: destination wind turbines

• k: wind speeds

• l: wind directions

This means that deficit_ijlk[0,1,2,3] holds the deficit caused by the first turbine on the second turbine for the third wind direction and fourth wind speed.

#here we import all wake deficit models available in PyWake
from py_wake import NOJ
from py_wake import Fuga
from py_wake import FugaBlockage
from py_wake import BastankhahGaussian

# Path to Fuga look-up tables
lut_path = os.path.dirname(py_wake.__file__)+'/tests/test_files/fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+00.nc'

models = {'NOJ': NOJ(site,windTurbines),
          'Fuga': Fuga(lut_path,site,windTurbines),
          'FugaBlockage': FugaBlockage(lut_path,site,windTurbines),
          'BGaus': BastankhahGaussian(site,windTurbines)
         }

In addition, these models can easily be combined with other models to better describe the wake behind the turbine, e.g. NOJ with linear sum superposition:

from py_wake.superposition_models import LinearSum

models['NOJLinear'] = NOJ(site,windTurbines,superpositionModel=LinearSum())

Or models can be combined in custom ways, e.g. NOJDeficit for the wake, LinearSum superposition and SelfSimilarityDeficit for the blockage:

from py_wake.wind_farm_models import All2AllIterative
from py_wake.deficit_models import NOJDeficit, SelfSimilarityDeficit

models['NOJ_ss'] = All2AllIterative(site,windTurbines,
                                          wake_deficitModel=NOJDeficit(),
                                          superpositionModel=LinearSum(),
                                          blockage_deficitModel=SelfSimilarityDeficit())

Wake Deficit Models - definition and configuration

First, we create some functions to visualize wake maps for different conditions and wake models.

from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from py_wake.deficit_models.deficit_model import WakeDeficitModel, BlockageDeficitModel
from py_wake.deficit_models.no_wake import NoWakeDeficit
from py_wake.site._site import UniformSite
from py_wake.flow_map import XYGrid
from py_wake.turbulence_models import CrespoHernandez
from py_wake.utils.plotting import setup_plot


#turbine diameter
D = 80

def get_flow_map(model=None, grid=XYGrid(x=np.linspace(-200, 500, 200), y=np.linspace(-200, 200, 200), h=70),
                 turbulenceModel=CrespoHernandez()):
    blockage_deficitModel = [None, model][isinstance(model, BlockageDeficitModel)]
    wake_deficitModel = [NoWakeDeficit(), model][isinstance(model, WakeDeficitModel)]
    wfm = All2AllIterative(UniformSite(), V80(), wake_deficitModel=wake_deficitModel, blockage_deficitModel=blockage_deficitModel,
                           turbulenceModel=turbulenceModel)
    return wfm(x=[0], y=[0], wd=270, ws=10, yaw=0).flow_map(grid)


def plot_deficit_map(model, cmap='Blues', levels=np.linspace(0, 10, 55)):
    fm = get_flow_map(model)
    fm.plot(fm.ws - fm.WS_eff, clabel='Deficit [m/s]', levels=levels, cmap=cmap, normalize_with=D)
    setup_plot(grid=False, ylabel="Crosswind distance [y/D]", xlabel= "Downwind distance [x/D]",
               xlim=[fm.x.min()/D, fm.x.max()/D], ylim=[fm.y.min()/D, fm.y.max()/D], axis='auto')


def plot_wake_deficit_map(model):
    cmap = np.r_[[[1,1,1,1],[1,1,1,1]],cm.Blues(np.linspace(-0,1,128))] # ensure zero deficit is white
    plot_deficit_map(model,cmap=ListedColormap(cmap))

def plot_blockage_deficit_map(model):
    from matplotlib import cm
    from matplotlib.colors import ListedColormap, LinearSegmentedColormap
    cmap = np.r_[cm.Reds_r(np.linspace(-0,1,127)),[[1,1,1,1],[1,1,1,1]],cm.Blues(np.linspace(-0,1,128))] # ensure zero deficit is white
    plot_deficit_map(model,cmap=ListedColormap(cmap), levels=np.linspace(-3.5,3.5,113))

NOJDeficit

The NOJDeficit model is implemented according to Niels Otto Jensen, “A note on wind generator interaction.” (1983), i.e. a top-hat wake.

Only valid in the far wake.

from py_wake.deficit_models import NOJDeficit
plot_wake_deficit_map(NOJDeficit())

TurboNOJDeficit

Modified definition of the wake expansion given by Nygaard [1], which assumes the wake expansion rate to be proportional to the local turbulence intensity in the wake. Here the local turbulence intensity is defined as the combination of ambient and wake added turbulence. Using the added wake turbulence model by Frandsen and integrating, an analytical expression for the wake radius can be obtained.

The definition in [1] of ambient turbulence is the free-stream TI and for this the model constant A has been tuned, however a fully consistent formulation of the model should probably use the locally effective TI, which includes the added turbulence from upstream turbines.

[1] Nygaard 2020 J. Phys.: Conf. Ser. 1618 062072 https://doi.org/10.1088/1742-6596/1618/6/062072

from py_wake.deficit_models import TurboNOJDeficit
plot_wake_deficit_map(TurboNOJDeficit())

FugaDeficit

The FugaDeficit model calculates the wake deficit based on a set of look-up tables computed by a linearized RANS solver. The look-up tables are created in advance using the Fuga GUI.

The FugaDeficit model can model the near wake, far wake and blockage deficit.

import py_wake
from py_wake.deficit_models import FugaDeficit

# Path to Fuga look-up tables
lut_path = os.path.dirname(py_wake.__file__)+'/tests/test_files/fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+00.nc'

plot_wake_deficit_map(FugaDeficit(lut_path))

BastankhahGaussianDeficit

The BastankhahGaussianDeficit model is implemented according to Bastankhah M and Porté-Agel F. “A new analytical model for wind-turbine wakes” J. Renew. Energy. 2014;70:116-23.

The model is valid in the far wake only.

from py_wake.deficit_models import BastankhahGaussianDeficit
plot_wake_deficit_map(BastankhahGaussianDeficit())

IEA37SimpleBastankhahGaussianDeficit

The IEA37SimpleBastankhahGaussianDeficit model is implemented according to the IEA task 37 documentation and is equivalent to BastankhahGaussianDeficit for \(beta=1/\sqrt{8} \sim ct=0.9637188\).

The model is valid in the far wake only.

from py_wake.deficit_models import IEA37SimpleBastankhahGaussianDeficit
plot_wake_deficit_map(IEA37SimpleBastankhahGaussianDeficit())

NiayifarGaussianDeficit

This wake deficit model accounts for the local turbulence intensity when evaluating the wake expansion. The expansion rate ‘k’ varies linearly with local turbluence intensity: k = a1 I + a2. The default constants are set according to publications by Porte-Agel’s group, which are based on LES simulations. Lidar field measurements by Fuertes et al. (2018) indicate that a = [0.35, 0.0] is also a valid selection.

It was implemented according to: Amin Niayifar and Fernando Porté-Agel, Analytical Modeling of Wind Farms: A New Approach for Power Prediction, Energies 2016, 9, 741; doi:10.3390/en9090741.

from py_wake.deficit_models import NiayifarGaussianDeficit
plot_wake_deficit_map(NiayifarGaussianDeficit(a=[0.38, 4e-3], use_effective_ti=False))

ZongGaussianDeficit

The wake expansion is also function of the local turbulence intensity and the wake width expression now follows the approach by Shapiro et al. (2018). It is an extension of the Niayifar et al. (2016) model implementation with the addition of the new wake width expression, which uses the near-wake length estimation by Vermeulen (1980).

It was implemented according to: Haohua Zong and Fernando Porté-Agel, A momentum-conserving wake superposition method for wind farm power prediction, J. Fluid Mech. (2020), vol. 889, A8; doi:10.1017/jfm.2020.77.

from py_wake.deficit_models import ZongGaussianDeficit
plot_wake_deficit_map(ZongGaussianDeficit(a=[0.38, 4e-3]))

CarbajoFuertesGaussianDeficit

Carbajo Fuertes et al. derived Gaussian wake model parameters from nacelle liadar measurements from a 2.5MW turbine and found a variation of epsilon with wake expansion, this in fact identical to the formulation by Zong, only that the near-wake length is fixed for Carbajo Fuertes at xth = 1.91 x/D. We took the relationships found by them and incorporated them into the Zong formulation. This wake deficit model includes an empirical correlation for epsilon as well as new constants for the wake expansion factor equation.

It is a modified version from the Zong deficit model with Gaussian constants, as seen in: Fernando Carbajo Fuertes, Corey D. Markfor and Fernando Porté-Agel, “Wind TurbineWake Characterization with Nacelle-Mounted Wind Lidars for Analytical Wake Model Validation”, Remote Sens. 2018, 10, 668; doi:10.3390/rs10050668.

from py_wake.deficit_models import CarbajofuertesGaussianDeficit
plot_wake_deficit_map(CarbajofuertesGaussianDeficit())

TurboGaussianDeficit

Implemented similar to Ørsted’s TurbOPark model (https://github.com/OrstedRD/TurbOPark)

from py_wake.deficit_models import TurboGaussianDeficit
plot_wake_deficit_map(TurboGaussianDeficit())

GCLDeficit

Implemented according to:

Larsen, G. C. (2009). A simple stationary semi-analytical wake model, Risø National Laboratory for Sustainable Energy, Technical University of Denmark. Denmark. Forskningscenter Risoe. Risoe-R, No. 1713(EN)

Description:

Based on an analytical solution of the thin shear layer approximation of the NS equations. The wake flow fields are assumed rotationally symmetric, and the rotor inflow fields are consistently assumed uniform. The effect of expansion is approximately accounted for by imposing suitable empirical downstream boundary conditions on the wake expansion that depend on the rotor thrust and the ambient turbulence conditions, respectively.

from py_wake.deficit_models import GCLDeficit
plot_wake_deficit_map(GCLDeficit())

Comparing wake deficit models

#printing all available wake deficit models in PyWake
from py_wake.utils.model_utils import get_models
from py_wake.deficit_models.deficit_model import WakeDeficitModel

deficitModels = get_models(WakeDeficitModel)
for deficitModel in deficitModels:
    print (deficitModel.__name__)

NOJDeficit
FugaDeficit
FugaMultiLUTDeficit
FugaYawDeficit
BastankhahGaussianDeficit
CarbajofuertesGaussianDeficit
IEA37SimpleBastankhahGaussianDeficit
NiayifarGaussianDeficit
TurboGaussianDeficit
ZongGaussianDeficit
GCLDeficit
NoWakeDeficit
NOJLocalDeficit
TurboNOJDeficit

1) Deficit along center line

plt.figure(figsize=((16,6)))
for i, deficitModel in enumerate(deficitModels):
    fm = get_flow_map(deficitModel(), XYGrid(x=np.linspace(-200,800,300), y=0))
    plt.plot(fm.x, 10-fm.WS_eff.squeeze(), ('-','--')[i//10], label=deficitModel.__name__)
setup_plot(title="Center line deficit", xlabel='x/D', ylabel='Deficit [m/s]')

2) Deficit profile downstream

for d in 2,8,20:
    plt.figure(figsize=((12,6)))
    for i, deficitModel in enumerate(deficitModels):
        fm = get_flow_map(deficitModel(), XYGrid(x=d*D, y=np.linspace(-200,200,300)))
        plt.plot(fm.y, 10-fm.WS_eff.squeeze(), ('-','--')[i//10], label=deficitModel.__name__)
    plt.title(f'Deficit profile {d}D downstream')
    plt.legend()

You can also implement your own deficit models

Deficit models must subclass DeficitModeland thus must implement the calc_deficit method and a class variable, args4deficit specifying the arguments required by its calc_deficit method

class DeficitModel(ABC):
    args4deficit = ['WS_ilk', 'dw_ijlk']

    @abstractmethod
    def calc_deficit(self):
        """Calculate wake deficit caused by the x'th most upstream wind turbines
        for all wind directions(l) and wind speeds(k) on a set of points(j)

        This method must be overridden by subclass

        Arguments required by this method must be added to the class list
        args4deficit

        See class documentation for examples and available arguments

        Returns
        -------
        deficit_ijlk : array_like
        """

from py_wake.deficit_models.deficit_model import WakeDeficitModel
from numpy import newaxis as na
class MyDeficitModel(WakeDeficitModel):

    def calc_deficit(self, WS_ilk, dw_ijlk, cw_ijlk,**_):
        # 30% deficit in downstream triangle
        ws_10pct_ijlk = 0.3*WS_ilk[:,na]
        triangle_ijlk = (self.wake_radius(dw_ijlk=dw_ijlk)>cw_ijlk)
        return ws_10pct_ijlk *triangle_ijlk

    def wake_radius(self, dw_ijlk, **_):
        return (.2*dw_ijlk)

plot_wake_deficit_map(MyDeficitModel())

Superposition models

The superposition models calculate the effective wind speed given the local wind speed and deficits (typically from multiple sources)

There are three different superposition models in PyWake:

• LinearSum: Deficits sum up linearly.

• SquaredSum: Deficits sum as root-sum-square.

• MaxSum: Only the largest deficit is considered.

LinearSum

In general, LinearSum should be used in combination with use_effective_ws=True which makes downstream wind turbines feel the effective wind speed (including wake effects from upstream turbines) instead of the ambient free-stream wind speed. Otherwise negative wind speeds may occur.

linear_sum = BastankhahGaussian(site, windTurbines, superpositionModel=LinearSum(), use_effective_ws=True)
plt.figure(figsize=(10,4))
linear_sum([0,200],[0,0],wd=270,ws=10).flow_map().plot_wake_map(levels=np.arange(0,10))
plt.xlabel('x [m]')
plt.ylabel('y [m]')

Text(0, 0.5, 'y [m]')

SquaredSum

The SquaredSum model is often used in combination with wake deficit models where the downstream wakes are scaled with the ambient free-stream velocity to avoid negative wind speeds. It is, however, a method to compensate for an inconsistent formulation, see section 2.2 in https://backend.orbit.dtu.dk/ws/portalfiles/portal/151671395/Park2_Documentation_and_Validation.pdf

from py_wake.superposition_models import SquaredSum
squared_sum = BastankhahGaussian(site, windTurbines, superpositionModel=SquaredSum())
plt.figure(figsize=(10,4))
squared_sum([0,200],[0,0],wd=270,ws=10).flow_map().plot_wake_map(levels=np.arange(0,10))
plt.xlabel('x [m]')
plt.ylabel('y [m]')

Text(0, 0.5, 'y [m]')

MaxSum

from py_wake.superposition_models import MaxSum
max_sum = BastankhahGaussian(site, windTurbines, superpositionModel=MaxSum())
plt.figure(figsize=(10,4))
max_sum([0,200],[0,0],wd=270,ws=10).flow_map().plot_wake_map(levels=np.arange(0,10))
plt.xlabel('x [m]')
plt.ylabel('y [m]')

Text(0, 0.5, 'y [m]')

Blockage deficit models

The blockage deficit models compute the blockage effects caused by a single wind turbine.

Their structure are quite similar to the wake deficit models. They model upstream blockage effects (wind speed reduction) and in addition, some models also models downstream speed-up effects. There are several blockage models available, which include:

• SelfSimilarityDeficit

• SelfSimilarityDeficit2020

• FugaDeficit

• VortexCylinder

• VortexDipole

• RankineHalfBody

• HybridInduction

• Rathmann

SelfSimilarityDeficit

Simple induction model, described in N. Troldborg, A.R. Meyer Fortsing, Wind Energy, 2016

from py_wake.deficit_models import SelfSimilarityDeficit
plot_blockage_deficit_map(SelfSimilarityDeficit())

SelfSimilarityDeficit2020

This is an updated version of N. Troldborg, A.R. Meyer Fortsing, Wind Energy, 2016. The new features are found in the radial and axial functions:

1. Radially Eq. (13) is replaced by a linear fit, which ensures the induction half width, r12, to continue to diminish approaching the rotor. This avoids unphysically large lateral induction tails, which could negatively influence wind farm simulations.

2. The value of gamma in Eq. (8) is revisited. Now gamma is a function of CT and axial coordinate to force the axial induction to match the simulated results more closely. The fit is valid over a larger range of thrust coefficients and the results of the constantly loaded rotor are excluded in the fit.

from py_wake.deficit_models import SelfSimilarityDeficit2020
plot_blockage_deficit_map(SelfSimilarityDeficit2020())

FugaDeficit

The FugaDeficit model calculates the wake deficit based on a set op look-up tables computed by a linearized RANS solver. The look-up tables be created in advance using the Fuga GUI.

The FugaDeficit models the near wake, far wake and blockage deficit effects.

Note, the present look-up table generator introduces some unphysical wriggles in the blockage deficit/speed-up.

from py_wake.deficit_models import FugaDeficit
plot_blockage_deficit_map(FugaDeficit())

VortexCylinder

Induced velocity from a semi-infinite cylinder of tangential vorticity, extending along the z axis.

This model is an adapted version of the one published by Emmanuel Branlard at https://github.com/ebranlard/wiz/blob/master/wiz/VortexCylinder.py

References:

• 5. Branlard, M. Gaunaa, Cylindrical vortex wake model: right cylinder, Wind Energy, 2014, https://onlinelibrary.wiley.com/doi/full/10.1002/we.1800

• 5. Branlard, Wind Turbine Aerodynamics and Vorticity Based Method, Springer, 2017

• 5. Branlard, A. Meyer Forsting, Using a cylindrical vortex model to assess the induction zone in front of aligned and yawed rotors, in Proceedings of EWEA Offshore Conference, 2015, https://orbit.dtu.dk/en/publications/using-a-cylindrical-vortex-model-to-assess-the-induction-zone-inf

from py_wake.deficit_models import VortexCylinder
plot_blockage_deficit_map(VortexCylinder())

VortexDipole

The vorticity originating from a wind turbine can be represented by a vortex dipole line (see Appendix B in [2]). The induction estimated by such a representation is very similar to the results given by the more complex vortex cylinder model in the far-field r/R > 6 [1,2]. The implementation follows the relationships given in [1,2]. This model is an adapted version of the one published by Emmanuel Branlard: https://github.com/ebranlard/wiz/blob/master/wiz/VortexDoublet.py

References: - [1] Emmanuel Branlard et al 2020 J. Phys.: Conf. Ser. 1618 062036 - [2] Branlard, E, Meyer Forsting, AR. Wind Energy. 2020; 23: 2068– 2086. https://doi.org/10.1002/we.2546

from py_wake.deficit_models import VortexDipole
plot_blockage_deficit_map(VortexDipole())

RankineHalfBody

A simple induction model using a Rankine Half Body to represent the induction introduced by a wind turbine. The source strength is determined enforcing 1D momentum balance at the rotor disc.

References:

• B Gribben, G Hawkes - A potential flow model for wind turbine induction and wind farm blockage - Technical Paper, Frazer-Nash Consultancy, 2019

from py_wake.deficit_models import RankineHalfBody
plot_blockage_deficit_map(RankineHalfBody())

HybridInduction

The idea behind this model originates from [2,3], which advocates to combine near-rotor and farfield approximations of a rotor’s induced velocities. Whereas in [1,2] the motivation is to reduce the computational effort, here the already very fast self-similar model [1] is combined with the vortex dipole approximation in the far-field, as the self-similar one is optimized for the near-field (r/R > 6, x/R < 1) and misses the acceleration around the wake for x/R > 0. The combination of both allows capturing the redistribution of energy by blockage. Location at which to switch from near-rotor to far-field can be altered though by setting switch_radius.

References: 1. N. Troldborg, A.R. Meyer Fortsing, Wind Energy, 2016 2. Emmanuel Branlard et al 2020 J. Phys.: Conf. Ser. 1618 062036 3. Branlard, E, Meyer Forsting, AR. Wind Energy. 2020; 23: 2068– 2086. https://doi.org/10.1002/we.2546

from py_wake.deficit_models import HybridInduction
plot_blockage_deficit_map(HybridInduction())

Rathmann

Ole Sten Rathmann (DTU) developed in 2020 an approximation to the vortex cylinder solution (E. Branlard and M. Gaunaa, 2014). In speed it is comparable to the vortex dipole method, whilst giving a flow-field nearly identical to the vortex cylinder model for x/R < -1. Its centreline deficit is identical to the vortex cylinder model, whilst using a radial shape function that depends on the opening of the vortex cylinder seen from a point upstream. To simulate the speed-up downstream the deficit is mirrored in the rotor plane with a sign change.

from py_wake.deficit_models import Rathmann
plot_blockage_deficit_map(Rathmann())

Comparing different blockage deficit models

from py_wake.utils.model_utils import get_models
from py_wake.deficit_models.deficit_model import BlockageDeficitModel

blockage_deficitModels = get_models(BlockageDeficitModel, exclude_None=True)
for deficitModel in blockage_deficitModels:
    print (deficitModel.__name__)

FugaDeficit
FugaMultiLUTDeficit
FugaYawDeficit
HybridInduction
SelfSimilarityDeficit2020
VortexDipole
RankineHalfBody
Rathmann
RathmannScaled
SelfSimilarityDeficit
VortexCylinder

1) Deficit along center line

plt.figure(figsize=((10,6)))
for i, deficitModel in enumerate(blockage_deficitModels):
    fm = get_flow_map(deficitModel(), XYGrid(x=np.linspace(-400,100,300), y=0))
    plt.plot(fm.x, 10-fm.WS_eff.squeeze(), ('-','--')[i//10], label=deficitModel.__name__)
setup_plot(title="Center line deficit", xlabel='x/D', ylabel='Deficit [m/s]')

2) Deficit profile 1 up- and downstream

for d in [-1,1]:
    plt.figure(figsize=((10,6)))
    for i, deficitModel in enumerate(blockage_deficitModels):
        fm = get_flow_map(deficitModel(), XYGrid(x=d*D, y=np.linspace(-200,200,300)))
        plt.plot(fm.y, 10-fm.WS_eff.squeeze(), ('-','--')[i//10], label=deficitModel.__name__)
    setup_plot(title="%sD %sstream"%(abs(d),('down','up')[d<0]), xlabel='x/D', ylabel='Deficit [m/s]')

Rotor-average models

In the plots below, it is clearly seen that the wind speed varies over the rotor, and that the the rotor-average wind speed is not well-defined by the wind sped at the rotor center. The rotor average model in PyWake defines one or more points at the turbine rotor to calculate a (weighted) rotor-average deficit. It includes:

• RotorCenter: One point at the center of the rotor

• GridRotorAvg: Custom grid in Cartesian coordinates’

• EqGridRotorAvg: Equidistant N x N Cartesian grid covering the rotor

• GQGridRotorAvg: M x N cartesian grid using Gaussian quadrature coordinates and weights

• PolarGridRotorAvg: Custom grid in polar coordinates

• CGIRotorAVG: Circular Gauss Integration

from py_wake.superposition_models import SquaredSum
from py_wake.flow_map import HorizontalGrid, YZGrid

R = D/2
wfm = BastankhahGaussian(site, windTurbines, superpositionModel=SquaredSum())
sim_res = wfm([0,200],[0,0],wd=270,ws=10)

fig,(ax1,ax2) = plt.subplots(1,2, figsize=(10,4))
ax1.set_xlabel("x [m]"), ax1.set_ylabel("y [m]")
sim_res.flow_map(HorizontalGrid(extend=.1)).plot_wake_map(10, ax=ax1, plot_colorbar=False)
sim_res.flow_map(YZGrid(x=199.99)).plot_wake_map(10, ax=ax2)
ax2.plot([-100,100],[70,70],'-.')
ax2.set_xlabel("y [m]"), ax1.set_ylabel("z [m]")

plt.figure()
flow_map = sim_res.flow_map(HorizontalGrid(x=[199.99], y=np.arange(-80, 80)))

for x in [-.5,.5]:
    plt.gca().axvline(x,color='grey')
plt.plot(flow_map.Y[:, 0]/D, flow_map.WS_eff_xylk[:, 0, 0, 0], '-.', label='Hub-height profile')
plt.plot([-.5,.5],[7.73,7.73],label='Rotor average: 7.73m/s')
rc_ws = flow_map.WS_eff_xylk[80, 0, 0, 0]
plt.plot(flow_map.Y[80, 0]/D, rc_ws,'.', ms=10, label='Rotor center: %.2fm/s'%rc_ws)
plt.legend()
plt.xlabel("y/D")
plt.ylabel('Wind speed [m/s]')

Text(0, 0.5, 'Wind speed [m/s]')

First we create a simple function to model all of the rotor-average models available in PyWake.

from py_wake.rotor_avg_models import RotorCenter, GridRotorAvg, EqGridRotorAvg, GQGridRotorAvg, CGIRotorAvg, PolarGridRotorAvg, PolarRotorAvg, polar_gauss_quadrature, GaussianOverlapAvgModel
from py_wake.flow_map import HorizontalGrid

R=D/2

def plot_rotor_avg_model(rotorAvgModel, name):
    plt.figure()
    m = rotorAvgModel
    wfm = BastankhahGaussian(site,windTurbines,rotorAvgModel=m)
    ws_eff = wfm([0, 200], [0, 0], wd=270, ws=10).WS_eff_ilk[1,0,0]
    plt.title(name)
    c = plt.scatter(m.nodes_x, m.nodes_y,c=m.nodes_weight,label="%.2fm/s"%(ws_eff))
    plt.colorbar(c,label='weight')
    plt.gca().add_artist(plt.Circle((0,0),1,fill=False))
    plt.axis('equal')
    plt.xlabel("y/R"); plt.ylabel('z/R')
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,1.5])
    plt.legend()

RotorCenter

Setting rotorAvgModel=None determines the rotor average wind speed from the rotor center point. Alternatively, you can use the RotorCenter model which is equivalent.

plot_rotor_avg_model(RotorCenter(), 'RotorCenter')

GridRotorAvg

The GridRotorAvg defines a set of points in cartesian coordinates.

x = y = np.array([-.6, 0,.6])
plot_rotor_avg_model(GridRotorAvg(x,y,nodes_weight = [0.25,.5,.25]), 'Grid_4')

EqGridRotorAvg

The EqGridRotorAvg defines a NxN equidistant cartesian grid of points and discards points outside the rotor.

plot_rotor_avg_model(EqGridRotorAvg(4), 'Grid_4')

GQGridRotorAvg

The GQGridRotorAvg defines a grid of M x N cartesian grid points using Gaussian quadrature coordinates and weights.

plot_rotor_avg_model(GQGridRotorAvg(4,3), 'GQGrid_4,3')

PolarGridRotorAvg

The PolarGridRotorAvg defines a grid in polar coordinates.

theta = np.linspace(-np.pi,np.pi,6, endpoint=False)
r = 2/3
plot_rotor_avg_model(PolarGridRotorAvg(r=r, theta=theta, r_weight=None, theta_weight=None), 'PolarGrid_6')

The polar grid can be combined with Gaussian Quadrature.

This is similar to the implementation in FusedWake: https://gitlab.windenergy.dtu.dk/TOPFARM/FUSED-Wake/-/blob/master/fusedwake/gcl/fortran/GCL.f

plot_rotor_avg_model(PolarGridRotorAvg(*polar_gauss_quadrature(4,10)), 'PolarGrid_4,10')

CGIRotorAvg

Circular Gauss integration with 4,7,9 or 21 points as defined in Abramowitz M and Stegun A. Handbook of Mathematical Functions. Dover: New York, 1970.

for n in [4,7,9,21]:
    plot_rotor_avg_model(CGIRotorAvg(n), 'CGIRotorAvg_%d'%n)

AreaOverlapModel

The AreaOverlapModel calculates the fraction of the downstream rotor that is covered by the wake from an upstream wind turbine. This model makes sense for tophat wake deficit models only, e.g. NOJDeficit.

The calculation formula can be found in Eq. (A1) of Feng, J., & Shen, W. Z. (2015). Solving the wind farm layout optimization problem using random search algorithm. Renewable Energy, 78, 182-192. https://doi.org/10.1016/j.renene.2015.01.005

|faee9d70701d40c5a8e196831715dd4c|

GaussianOverlapAvgModel

The GaussianOverlapModel computes the integral of the gaussian wake deficit over the downstream rotor. To speed up the computation, normalized integrals have been precalculated and stored in a look-up table. This model need the gaussian width parameter, \(\sigma\), and therefore only applies to gaussian wake deficit models. See Appendix A in https://github.com/OrstedRD/TurbOPark/blob/main/TurbOPark%20description.pdf

|7ce86e2b832242ea868974ce4f40acce|

Comparing rotor-average models

In general, the compuational cost and the accuracy of the estimate increases with the number of points, but the distribution of the points also has an impact.

The plot below shows the absolute error of the estimated rotor-average wind speed for the wind directions 270 \(\pm\) 30\(^\circ\) (i.e. wind directions with more than 1% deficit) for a number of different rotor-average models.

grid_models = [EqGridRotorAvg(i) for i in range(1,10)]
wd_lst = np.arange(240,301)

def get_ws_eff(rotorAvgModel):
    wfm = BastankhahGaussian(site,windTurbines,rotorAvgModel=rotorAvgModel)
    return wfm([0, 200], [0, 0], wd=wd_lst, ws=10).WS_eff_ilk[1,:,0]

ws_ref = get_ws_eff(EqGridRotorAvg(200)) # Use 200x200 points (31700 inside the rotor) to determine the reference value

def get_n_err(rotorAvgModel):
    ws_mean_err = np.abs(get_ws_eff(rotorAvgModel) - ws_ref).mean()
    return len(rotorAvgModel.nodes_x), ws_mean_err


plt.gca().axhline(0, color='grey',ls='--')
plt.plot(*zip(*[get_n_err(m) for m in grid_models]), label='Grid_x')
model_lst = [('RotorCenter', EqGridRotorAvg(1)),
             ('Grid_4', EqGridRotorAvg(4)),
             ('PolarGrid_4,10', PolarRotorAvg(*polar_gauss_quadrature(4,10))),
             ('GQGrid_4,3', GQGridRotorAvg(4, 3))] + \
            [('CGI_%d'%n, CGIRotorAvg(n)) for n in [4,7,9,21]]
for name, model in model_lst:
    n,err = get_n_err(model)
    plt.plot(n,err,'.',ms=10, label="%s (%.4f)"%(name,err))
goam_err = np.abs(get_ws_eff(GaussianOverlapAvgModel()) - ws_ref).mean()
plt.plot([0],[goam_err],'.', ms=10, label="GaussianOverlapAvgModel (%.4f)"%(goam_err))
plt.xlabel('Number of rotor-average points')
plt.ylabel(r'Mean abs error (270$\pm30^\circ$) [m/s]')
plt.legend()

<matplotlib.legend.Legend at 0x7f78300d7220>

Deflection models

The deflection models calculate the deflection of the wake due to yaw-misalignment, sheared inflow etc.

Note, this is one of the four effects of skew inflow that is handled in PyWake, see here.

The deflection models take as input the downwind and crosswind distances between the source wind turbines and the destination wind turbines/sites and calculate a new set of deflected downwind and crosswind distances. This type of model is important for simulations where the turbine experiences a change in angle between the incoming flow and the rotor, for example in active yaw control or wake steering optimization.

In PyWake, there are three different wake deflection models:

• JimenezWakeDeflection

• FugaDeflection

• GCLHillDeflection

First we create a simple function to plot the different deflection models available in PyWake.

def plot_deflection(deflectionModel):
    from py_wake import BastankhahGaussian
    from py_wake.examples.data.iea37._iea37 import IEA37Site, IEA37_WindTurbines

    site = IEA37Site(16)
    x, y = [0, 400, 800], [0, 0, 0]
    windTurbines = V80()
    wfm = BastankhahGaussian(site, windTurbines, deflectionModel=deflectionModel)

    yaw = [-20,20,0]
    D = windTurbines.diameter()

    plt.figure(figsize=(14,4))
    fm = wfm(x, y, yaw=yaw, tilt=0, wd=270, ws=10).flow_map()
    fm.plot_wake_map(normalize_with=D)
    center_line = fm.min_WS_eff()
    plt.plot(center_line.x/D, center_line/D,'--k')
    plt.grid()

JimenezWakeDeflection

The JimenezWakeDeflection model is implemented according to Jiménez, Á., Crespo, A. and Migoya, E. (2010), Application of a LES technique to characterize the wake deflection of a wind turbine in yaw. Wind Energ., 13: 559-572. doi:10.1002/we.380.

It is the most common wake deflection model used in PyWake and has proved to result in a decent representation of the skewed inflow behind the turbine rotor. In the study, the authors carried out a Large Eddy Simulation (LES) to characterize the turbulence behind a wind turbine given the wake deflection created by different yaw angle and thrust coefficient settings.

from py_wake.deflection_models import JimenezWakeDeflection
plt.figure()
plot_deflection(JimenezWakeDeflection())
plt.xlabel('x/D')
plt.ylabel('y/D')

/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:562: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imin = index.get_loc(minval, method="nearest")
/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:563: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imax = index.get_loc(maxval, method="nearest")

Text(0, 0.5, 'y/D')

<Figure size 640x480 with 0 Axes>

FugaDeflection

from py_wake.deflection_models import FugaDeflection
plt.figure()
plot_deflection(FugaDeflection())
plt.xlabel('x/D')
plt.ylabel('y/D')

/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:562: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imin = index.get_loc(minval, method="nearest")
/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:563: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imax = index.get_loc(maxval, method="nearest")

Text(0, 0.5, 'y/D')

<Figure size 640x480 with 0 Axes>

GCLHillDeflection

Deflection model based on Hill’s ring vortex theory.

Implemented according to: Larsen, G. C., Ott, S., Liew, J., van der Laan, M. P., Simon, E., R.Thorsen, G., & Jacobs, P. (2020). Yaw induced wake deflection - a full-scale validation study. Journal of Physics - Conference Series, 1618, [062047]. https://doi.org/10.1088/1742-6596/1618/6/062047.

Note, this model uses the wake centerline deficit magnitude to calculate the deflection. Hence non-gaussian-shaped wake deficit models as well as deficit models with improper near wake fields, e.g. NOJDeficit, BastankhahGaussianDeficit, NiayifarGaussianDeficit, should be avoided.

As default the model uses the WakeDeficitModel specified for the WindFarmModel to calculate the magnitude of the deficit, but a separate model can be specified.

from py_wake.deflection_models import GCLHillDeflection
plt.figure()
plot_deflection(GCLHillDeflection())
plt.xlabel('x/D')
plt.ylabel('y/D')

/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:562: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imin = index.get_loc(minval, method="nearest")
/opt/conda/lib/python3.9/site-packages/xarray/core/missing.py:563: FutureWarning: Passing method to Float64Index.get_loc is deprecated and will raise in a future version. Use index.get_indexer([item], method=...) instead.
  imax = index.get_loc(maxval, method="nearest")

Text(0, 0.5, 'y/D')

<Figure size 640x480 with 0 Axes>

You can also implement your own deflection models

Deficit models must subclass DeficitModel and thus must implement the calc_deflection method and a class variable, args4deflection specifying the arguments required by its calc_deflection method.

class DeflectionModel(ABC):
    args4deflection = ["ct_ilk"]

    @abstractmethod
    def calc_deflection(self, dw_ijl, hcw_ijl, **kwargs):
        """Calculate deflection

        This method must be overridden by subclass

        Arguments required by this method must be added to the class list
        args4deflection

        See documentation of EngineeringWindFarmModel for a list of available input arguments

        Returns
        -------
        dw_ijlk : array_like
            downwind distance from source wind turbine(i) to destination wind turbine/site (j)
            for all wind direction (l) and wind speed (k)
        hcw_ijlk : array_like
            horizontal crosswind distance from source wind turbine(i) to destination wind turbine/site (j)
            for all wind direction (l) and wind speed (k)
        """

from py_wake.deflection_models import DeflectionModel
from numpy import newaxis as na

class MyDeflectionModel(DeflectionModel):

    def calc_deflection(self, dw_ijlk, hcw_ijlk, **_):
        hcw_ijlk =hcw_ijlk + .1*dw_ijlk       # deflect 1/10 of the downstream distance
        dh_ijlk =np.zeros_like(hcw_ijlk)      # no vertical deflection
        return dw_ijlk, hcw_ijlk, dh_ijlk

iea_my_deflection = BastankhahGaussian(site, windTurbines, deflectionModel=MyDeflectionModel())

plt.figure(figsize=(10,4))
iea_my_deflection(x=[0], y=[0], wd=270, ws=10).flow_map().plot_wake_map(normalize_with=D)
plt.xlabel('x/D')
plt.ylabel('y/D');

Turbulence models

The turbulence models in PyWake are used to calculate the added turbulence in the wake from one wind turbine to downstream turbines or sites in the wind farm. These are important when the flow properties behind the rotor must be accurately represented, for example for calculation of fatigue loading of turbine components.

First we create a simple function to plot the different turbulence models available in PyWake.

def plot_turb_map(model, cmap='Blues'):
    fm = get_flow_map(NOJDeficit(),turbulenceModel=model)
    fm.plot(fm.TI_eff, clabel="Added turbulence intensity [-]", levels=100, cmap=cmap, normalize_with=D)
    setup_plot(grid=False, ylabel="Crosswind distance [y/D]", xlabel="Downwind distance [x/D]",
               xlim=[fm.x.min() / D, fm.x.max() / D], ylim=[fm.y.min() / D, fm.y.max() / D], axis='auto')

STF2005TurbulenceModel

Steen Frandsen model implemented according to IEC61400-1, 2005 and weight according to Steen Frandsen’s thesis.

from py_wake.turbulence_models import STF2005TurbulenceModel
plot_turb_map(STF2005TurbulenceModel())

STF2017TurbulenceModel

Steen Frandsen model implemented according to IEC61400-1, 2017 and weight according to Steen Frandsen’s thesis.

from py_wake.turbulence_models import STF2017TurbulenceModel
plot_turb_map(STF2017TurbulenceModel())

STF20XXTurbulenceModel with IEC-based spread angle

The STF2005TurbulenceModel and STF2017TurbulenceModel take a weight_function input which defaults to the bell-shaped FrandsenWeight defined in Steen Frandsen’s thesis. As an alternative the IECWeight applies the full added turbulence in a 21.6\(^\circ\) spread angle up to 10 diameter downstream.

Note, this is a debatable interpretation of the IEC standard which includes a 6% contribution from neighbouring wind turbines when calculating the omni-directional effective turbulence intensity. These 6% maps to a spread angle of 360\(^\circ\cdot\) 6% = 21.6\(^\circ\).

Note, the IEC standard includes more concepts which is not implemented in PyWake.

from py_wake.turbulence_models import STF2017TurbulenceModel, IECWeight
from py_wake.superposition_models import SqrMaxSum
plot_turb_map(STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=SqrMaxSum(),
                                     weight_function=IECWeight(distance_limit=10)))

GCLTurbulence

Gunner Chr. Larsen model implemented according to:

Pierik, J. T. G., Dekker, J. W. M., Braam, H., Bulder, B. H., Winkelaar, D., Larsen, G. C., Morfiadakis, E., Chaviaropoulos, P., Derrick, A., & Molly, J. P. (1999). European wind turbine standards II (EWTS-II). In E. L. Petersen, P. Hjuler Jensen, K. Rave, P. Helm, & H. Ehmann (Eds.), Wind energy for the next millennium. Proceedings (pp. 568-571). James and James Science Publishers.

from py_wake.turbulence_models import GCLTurbulence
plot_turb_map(GCLTurbulence())

CrespoHernandez

Implemented according to: A. Crespo and J. Hernández, Turbulence characteristics in wind-turbine wakes, J. of Wind Eng. and Industrial Aero. 61 (1996) 71-85.

from py_wake.turbulence_models import CrespoHernandez
plot_turb_map(CrespoHernandez())

Comparing turbulence models

#printing all available wake deficit models in PyWake
from py_wake.utils.model_utils import get_models
from py_wake.turbulence_models.turbulence_model import TurbulenceModel

turbulenceModels = get_models(TurbulenceModel, exclude_None=True)

for model in turbulenceModels:
    print (model.__name__)

CrespoHernandez
GCLTurbulence
STF2005TurbulenceModel
STF2017TurbulenceModel

1) Turbulence intensity along center line

plt.figure(figsize=((16,6)))
for i, model in enumerate(turbulenceModels):
    fm = get_flow_map(NOJDeficit(), turbulenceModel=model(), grid=XYGrid(x=np.linspace(-200,800,300), y=0))
    plt.plot(fm.x, fm.TI_eff.squeeze(), ('-','--')[i//10], label=model.__name__)
setup_plot(title="Center line added turbulence", xlabel='x/D', ylabel='Add turbulence intensity [-]')

2) Deficit profile 2D downstream

d = 2
plt.figure(figsize=((10,6)))
for i, model in enumerate(turbulenceModels):
    fm = get_flow_map(NOJDeficit(), turbulenceModel=model(), grid=XYGrid(x=d*D, y=np.linspace(-200,200,300)))
    plt.plot(fm.y, fm.TI_eff.squeeze(), ('-','--')[i//10], label=model.__name__)
setup_plot(title="%sD %sstream"%(abs(d),('down','up')[d<0]), xlabel='x/D', ylabel='Add turbulence intensity [-]')

Ground models

The ground models in PyWake are used to model the effects that the ground has on the inflow and wake.

Mirror

The Mirror ground model lets the ground mirror the wake deficit. It is implemented by adding wakes from underground mirrored wind turbines

from py_wake.ground_models import Mirror
wfm = NOJ(site, windTurbines, k=.5, groundModel=Mirror())

plt.figure()
wfm([0], [0], wd=0).flow_map(YZGrid(x=0, y=np.arange(-300, 100, 1) + .1, z=np.arange(-100, 200))).plot_wake_map()
plt.xlabel('x [m]')
plt.ylabel('y [m]')

Text(0, 0.5, 'y [m]')