Site Object

For a given position, reference wind speed (WSref) and wind direction (WDref), Site provides the local wind condition in terms of wind speed (WS), wind direction (WD), turbulence intensity (TI) and the probability of each combination of wind direction and wind speed. Furthermore, Site is responsible for calculating the down-wind, cross-wind and vertical distance between wind turbines (which in non-flat terrain is different from the straight-line distances).

Intall 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

Predefined example sites

PyWake contains a few predefined sites of different complexities:

IEA37Site: UniformSite (fix wind speed (9.8m/s), predefined wind sector probability).
Hornsrev1: UniformWeibullSite (Weibull distributed wind speed, predefined wind sector propability, uniform wind a over flat wind area).
ParqueFicticioSite: WaspGridSite (position-dependent Weibull distributed wind speed and sector probability. Terrain following distances over non-flat terrain). Loaded from a set of *.grd files exported from WAsP.

First we import all sites and Python elements for later use

[2]:

import numpy as np
import matplotlib.pyplot as plt

[3]:

from py_wake.examples.data.hornsrev1 import Hornsrev1Site
from py_wake.examples.data.iea37 import IEA37Site
from py_wake.examples.data.ParqueFicticio import ParqueFicticioSite

sites = {"IEA37": IEA37Site(n_wt=16),
         "Hornsrev1": Hornsrev1Site(),
         "ParqueFicticio": ParqueFicticioSite()}

PyWake also allows for user-defined sites

You can define your own site using one of the Site classes:

UniformWeibullSite: Site with uniform sector-dependent Weibull distributed wind speed.
WaspGridSite: Site with gridded non-uniform inflow based on *.grd files exported from WAsP.
XRSite: The flexible general base class behind all Sites.

For more information on these classes, please see the API reference on the Site object.

UniformWeibullSite

[4]:

from py_wake.site import UniformWeibullSite

#specifying the necessary parameters for the UniformWeibullSite object
site = UniformWeibullSite(p_wd = [.20,.25,.35,.25],                         # sector frequencies
                          a = [9.176929,  9.782334,  9.531809,  9.909545],  # Weibull scale parameter
                          k = [2.392578, 2.447266, 2.412109, 2.591797],     # Weibull shape parameter
                          ti = 0.1                                          # turbulence intensity, optional
                         )

WaspGridSite

[5]:

from py_wake.site import WaspGridSite
from py_wake.examples.data.ParqueFicticio import ParqueFicticio_path

site = WaspGridSite.from_wasp_grd(ParqueFicticio_path)

XRSite

The XRSite is the most general and flexible Site. For the input dataset there are some required and optional data variables, such as:

Required data variables:
- P: probability of flow case(s)
or
- Weibull_A: Weibull scale parameter(s)
- Weibull_k: Weibull shape parameter(s)
- Sector_frequency: Probability of each wind direction sector
Optional data variables:
- WS: Wind speed, if not present, the reference wind speed ws is used
- Speedup: Factor multiplied to the wind speed
- Turning: Wind direction turning
- TI: Turbulence intensity
- xxx: Custom variables needed by the wind turbines to compute power, ct or loads
Each data variable may be constant or depend on a combination of the following inputs (Note, the input variables must be ordered according to the list, i.e. P(wd,ws) is ok, while P(ws,wd) is not):
- i: Wind turbine position (one position per wind turbine)
- x,y: Gridded 2d position
- x,y,h: Gridded 3d position
- time: Time
- wd: Refernce wind direction
- ws : Reference wind speed

[6]:

from py_wake.site import XRSite
from py_wake.site.shear import PowerShear
import xarray as xr
import numpy as np
from py_wake.utils import weibull
from numpy import newaxis as na

f = [0.036, 0.039, 0.052, 0.07, 0.084, 0.064, 0.086, 0.118, 0.152, 0.147, 0.1, 0.052]
A = [9.177, 9.782, 9.532, 9.91, 10.043, 9.594, 9.584, 10.515, 11.399, 11.687, 11.637, 10.088]
k = [2.393, 2.447, 2.412, 2.592, 2.756, 2.596, 2.584, 2.549, 2.471, 2.607, 2.627, 2.326]
wd = np.linspace(0, 360, len(f), endpoint=False)
ti = .1

# Site with constant wind speed, sector frequency, constant turbulence intensity and power shear
uniform_site = XRSite(
    ds=xr.Dataset(data_vars={'WS': 10, 'P': ('wd', f), 'TI': ti},
                  coords={'wd': wd}),
    shear=PowerShear(h_ref=100, alpha=.2))

# Site with wind direction dependent weibull distributed wind speed
uniform_weibull_site = XRSite(
    ds=xr.Dataset(data_vars={'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti},
                  coords={'wd': wd}))

# Site with a speedup and a turning value per WT
x_i, y_i = np.arange(5) * 100, np.zeros(5)  # WT positions

complex_fixed_pos_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': ('i', np.arange(.8, 1.3, .1)),
                   'Turning': ('i', np.arange(-2, 3)),
                   'P': ('wd', f)},
        coords={'i': np.arange(5), 'wd': wd}),
    initial_position=np.array([x_i, y_i]).T)

# Site with gridded speedup information
complex_grid_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))),
                   'P': ('wd', f)},
        coords={'x': [0, 500, 1000], 'y': [0, 500], 'wd': wd}))

# Site with ws dependent speedup and wd- and ws distributed probability
P_ws = weibull.cdf(np.array([3, 5, 7, 9, 11, 13]), 10, 2) - weibull.cdf(np.array([0, 3, 5, 7, 9, 11]), 10, 2)
P_wd_ws = P_ws[na, :] * np.array(f)[:, na]

complex_ws_site = XRSite(
    ds=xr.Dataset(
        data_vars={'Speedup': (['ws'], np.arange(.8, 1.4, .1)),
                   'P': (('wd', 'ws'), P_wd_ws), 'TI': ti},
        coords={'ws': [1.5, 4, 6, 8, 10, 12], 'wd': wd}))

Gridded non-uniform wind resources as time series with XRSite

This example creates a site with 2D non-uniform time series resources for wind speed (ws), wind direction (wd), and turbulence intensity (ti). The data is provided as arrays with dimensions corresponding to spatial coordinates (x, y) and time.

[7]:

site_x, site_y = np.meshgrid(np.arange(0.1, 1000, 100), np.arange(0.1, 2000, 100))
site_x, site_y = site_x.flatten(), site_y.flatten()
site_time = np.arange(100)
site_ws = np.random.uniform(3.0, 21.0, (len(site_x), len(site_y), len(site_time)))
site_wd = np.random.uniform(0.0, 360.0, (len(site_x), len(site_y), len(site_time)))
ds = xr.Dataset(
    data_vars=dict(
        WS=(["x", "y", "time"], site_ws),
        WD=(["x", "y", "time"], site_wd),
        TI=(["x", "y", "time"], np.ones_like(site_ws) * 0.1),  # hardcoded TI=0.1
        P=1,  # deterministic wind resource
    ),
    coords=dict(
        x=("x", site_x),
        y=("y", site_y),
        time=("time", site_time),
    ),
)
non_uniform_ts_site = XRSite(ds)
wss_at_mean_loc = non_uniform_ts_site.local_wind(
    site_x.mean(), site_y.mean(), time=site_time
)["WS_ilk"]
print(f"Mean wind speed at the mean location: {wss_at_mean_loc.mean():.2f} m/s")

# check the map of the mean wind speed at the site
mean_resource = ds.WS.mean(dim="time").values
plt.contourf(site_x, site_y, mean_resource)
plt.colorbar()
plt.title("Mean wind speed [m/s]")
plt.xlabel("x [m]")
_ = plt.ylabel("y [m]")

Mean wind speed at the mean location: 12.08 m/s

Wake effects from neighbouring wind farms

In some cases, calculation of wake interaction between the wind farm to optimize and neighbouring wind farms considerably slow down an optimization work flow. To avoid this, a site, which includes wake effects from neighbouring wind farms, can be pre-generated and used for the optimization.

The speed up of this solution depends on the number of turbines in both the current and neighbouring wind farms, as well as the type of sites. If the original site is a uniform site, then a pre-generated site with wake effects from neighbouring wind farms may slow down the workflow as it adds interpolation of inflow characteristics in space.

Note also, that a pre-generated site with wake effects from neighbouring wind farms is only eqivalent to the full simulation if the applied deficit model uses the effective wind speed (some models have an option to switch between effective and free-stream local wind speed).

[8]:

# import and setup site and windTurbines
from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines
from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit
from py_wake.wind_turbines import WindTurbine, WindTurbines
from py_wake.wind_farm_models import PropagateDownwind
from py_wake.superposition_models import LinearSum

site = IEA37Site(16)

# setup current, neighbour and all positions
wt_x, wt_y = site.initial_position.T
neighbour_x, neighbour_y = wt_x-4000, wt_y
all_x, all_y = np.r_[wt_x,neighbour_x], np.r_[wt_y,neighbour_y]

windTurbines = WindTurbines.from_WindTurbine_lst([IEA37_WindTurbines(),IEA37_WindTurbines()])
windTurbines._names = ["Current wind farm","Neighbour wind farm"]
types = [0]*len(wt_x) + [1]*len(neighbour_x)

wf_model = PropagateDownwind(site, windTurbines,
                             wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True),
                             superpositionModel=LinearSum())

# Consider wd=270 +/- 30 deg only
wd_lst = np.arange(240,301)

[9]:

#plotting the wake maps for the desired flow case
wsp = 9.8
wdir = 267

plt.figure(figsize=(16, 6))
wf_model(all_x, all_y, type=types, wd=wdir, ws=wsp, h=110).flow_map().plot_wake_map()
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.title('Wake map for'+ f' {wdir} deg and {wsp} m/s')

[9]:

Text(0.5, 1.0, 'Wake map for 267 deg and 9.8 m/s')

Now, we run the simulation of all wind turbines and calculate AEP of current wind farm

[10]:

print("Total AEP: %f GWh"%wf_model(all_x, all_y, type=types, ws=[wsp], wd=wd_lst).aep().isel(wt=np.arange(len(wt_x))).sum())

Total AEP: 85.187662 GWh

We can also calculate the AEP of the current wind farm by enclosing it in a flow box and setting up a new wind farm model

[11]:

#making a flow box covering the area of interest (i.e the current wind farm + 100m)

ext = 1000
flow_box = wf_model(neighbour_x, neighbour_y, wd=wd_lst).flow_box(
    x=np.linspace(min(wt_x) - ext, max(wt_x) + ext, 101),
    y=np.linspace(min(wt_y) - ext, max(wt_y) + ext, 101),
    h=110)

#creating new site based on the flow box

from py_wake.site.xrsite import XRSite
wake_site = XRSite.from_flow_box(flow_box)

Now, we plot the “free-stream” inflow wind speed of the current wind farm.

[12]:

plt.figure(figsize=(16, 6))
wake_site.ds.WS.sel(wd=267).plot(y='y', cmap = 'Blues_r')
windTurbines.plot(all_x, all_y, types, wd=270)

Then, we setup a new wind farm model with the new pre-generated site and calculate the AEP.

[13]:

wf_model_wake_site = PropagateDownwind(wake_site, windTurbines,
                                       wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True),
                                       superpositionModel=LinearSum())

[14]:

print("Total AEP: %f GWh"%wf_model_wake_site(wt_x, wt_y, ws=[wsp], wd=wd_lst).aep().sum())

Total AEP: 85.187032 GWh

Note that the AEP is not exactly equal due to interpolation errors. The discrepancy can be lowered by increasing the resolution of the flow box.

Lastly, we plot the flow map of the current wind farm with the selected flow box.

[15]:

plt.figure(figsize=(16, 6))
wf_model_wake_site(wt_x, wt_y, wd=wdir, ws=wsp, h=110).flow_map().plot_wake_map()
windTurbines.plot(neighbour_x, neighbour_y, type=1, wd=wdir)
plt.xlabel('x [m]')
plt.ylabel('y [m]')

[15]:

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

Distance

We can also calculate the distance between points of a specific site for either flat or complex terrain.

For the IEA37Site and the Hornsrev1 sites the distances between points are straight line distances, as these sites are characterized by flat terrain.

For the ParqueFicticioSite, on the other hand, the down-wind distance is larger as it follows the non-flat terrain.

[22]:

wd = [0, 30,90] # wind direction at source

for name, site in sites.items():
    print ("------- %s -------"%name)
    wt_x, wt_y = site.initial_position[0]
    site.distance.setup(src_x_ilk=[wt_x, wt_x], src_y_ilk=[wt_y, wt_y-1000], src_h_ilk=[70,90], src_z_ilk=[0,0]) # wt2 1000m to the south
    dw_ijlk, cw_ijlk, dh_ijlk = site.distance(wd_l=wd, src_idx=[0], dst_idx=[[1,1,1]])


    print ('Wind direction: \t\t%d deg\t\t%d deg\t\t%d deg'%tuple(wd))
    print ('Down wind distance [m]: \t%.1f\t\t%.1f\t\t%.1f'%tuple(dw_ijlk[0,0,:,0]))
    print ('Cross wind distance [m]: \t%.1f\t\t%.1f\t\t%.1f'%tuple(cw_ijlk[0,0,:,0]))
    print ('Height difference [m]: \t\t%.1f\t\t%.1f\t\t%.1f'%tuple(dh_ijlk[0,0,:,0]))
    print()

------- IEA37 -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         1000.0          866.0           0.0
Cross wind distance [m]:        0.0             500.0           1000.0
Height difference [m]:          20.0            20.0            20.0

------- Hornsrev1 -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         1000.0          866.0           0.0
Cross wind distance [m]:        0.0             500.0           1000.0
Height difference [m]:          20.0            20.0            20.0

------- ParqueFicticio -------
Wind direction:                 0 deg           30 deg          90 deg
Down wind distance [m]:         1023.6          886.5           -0.0
Cross wind distance [m]:        0.0             500.0           1000.0
Height difference [m]:          20.0            20.0            20.0

Wind resource distribution plots

The Site object has a few plot function to visualize its properties, mainly the wind resource given by the wind rose and the probability functions.

[23]:

import matplotlib.pyplot as plt
site = sites['Hornsrev1']

Plotting wind rose.

[24]:

_ = site.plot_wd_distribution(n_wd=12, ws_bins=[0,5,10,15,20,25])

Plotting probability density function for the four sectors studied.

[25]:

_ = site.plot_ws_distribution(wd=[0,90,180,270])

Plotting probablity density function for the four sector studied \(\pm\) 45 degrees.

If include_wd_distribution=true, the wind speed probability distributions are multiplied by the wind direction probability.

The sector size is set to 360 / len(wd). This only makes sense if the wd array is evenly distributed

[26]:

_ = site.plot_ws_distribution(wd=[0,90,180,270], include_wd_distribution=True)