{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Site Object"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Intall PyWake if needed**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Install PyWake if needed\n",
"try:\n",
" import py_wake\n",
"except ModuleNotFoundError:\n",
" !pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Predefined example sites**\n",
"\n",
"PyWake contains a few predefined sites of different complexities:\n",
"\n",
"- **IEA37Site**: `UniformSite` (fix wind speed (9.8m/s), predefined wind sector probability).\n",
"- **Hornsrev1**: `UniformWeibullSite` (Weibull distributed wind speed, predefined wind sector propability, uniform wind a over flat wind area).\n",
"- **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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**First we import all sites and Python elements for later use**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from py_wake.examples.data.hornsrev1 import Hornsrev1Site\n",
"from py_wake.examples.data.iea37 import IEA37Site\n",
"from py_wake.examples.data.ParqueFicticio import ParqueFicticioSite\n",
"\n",
"sites = {\"IEA37\": IEA37Site(n_wt=16), \n",
" \"Hornsrev1\": Hornsrev1Site(), \n",
" \"ParqueFicticio\": ParqueFicticioSite()}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**PyWake also allows for user-defined sites**\n",
"\n",
"You can define your own site using one of the `Site` classes:\n",
"\n",
"- [UniformWeibullSite](#UniformWeibullSite): Site with uniform sector-dependent Weibull distributed wind speed.\n",
"- [WaspGridSite](#WaspGridSite): Site with gridded non-uniform inflow based on *.grd files exported from WAsP.\n",
"- [XRSite](#XRSite): The flexible general base class behind all Sites.\n",
"\n",
"For more information on these classes, please see the [API reference on the Site object](https://topfarm.pages.windenergy.dtu.dk/PyWake/api/Site.html)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### UniformWeibullSite"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from py_wake.site import UniformWeibullSite\n",
"\n",
"#specifying the necessary parameters for the UniformWeibullSite object\n",
"site = UniformWeibullSite(p_wd = [.20,.25,.35,.25], # sector frequencies\n",
" a = [9.176929, 9.782334, 9.531809, 9.909545], # Weibull scale parameter\n",
" k = [2.392578, 2.447266, 2.412109, 2.591797], # Weibull shape parameter\n",
" ti = 0.1 # turbulence intensity, optional\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### WaspGridSite"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from py_wake.site import WaspGridSite\n",
"from py_wake.examples.data.ParqueFicticio import ParqueFicticio_path\n",
"\n",
"site = WaspGridSite.from_wasp_grd(ParqueFicticio_path)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### XRSite\n",
"\n",
"The `XRSite` is the most general and flexible `Site`. For the input dataset there are some required and optional data variables, such as:\n",
"\n",
"- Required data variables:\n",
" - `P`: probability of flow case(s)\n",
"\n",
" or\n",
"\n",
" - `Weibull_A`: Weibull scale parameter(s)\n",
" - `Weibull_k`: Weibull shape parameter(s)\n",
" - `Sector_frequency`: Probability of each wind direction sector\n",
" \n",
"\n",
"- Optional data variables:\n",
"\n",
" - `WS`: Wind speed, if not present, the reference wind speed `ws` is used\n",
" - `Speedup`: Factor multiplied to the wind speed\n",
" - `Turning`: Wind direction turning\n",
" - `TI`: Turbulence intensity\n",
" - xxx: Custom variables needed by the wind turbines to compute power, ct or loads\n",
" \n",
"\n",
"- 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):\n",
"\n",
" - `i`: Wind turbine position (one position per wind turbine)\n",
" - `x`,`y`: Gridded 2d position\n",
" - `x`,`y`,`h`: Gridded 3d position\n",
" - `time`: Time\n",
" - `wd`: Refernce wind direction\n",
" - `ws` : Reference wind speed"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from py_wake.site import XRSite\n",
"from py_wake.site.shear import PowerShear\n",
"import xarray as xr\n",
"import numpy as np\n",
"from py_wake.utils import weibull\n",
"from numpy import newaxis as na\n",
"\n",
"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]\n",
"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]\n",
"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]\n",
"wd = np.linspace(0, 360, len(f), endpoint=False)\n",
"ti = .1\n",
"\n",
"# Site with constant wind speed, sector frequency, constant turbulence intensity and power shear\n",
"uniform_site = XRSite(\n",
" ds=xr.Dataset(data_vars={'WS': 10, 'P': ('wd', f), 'TI': ti},\n",
" coords={'wd': wd}),\n",
" shear=PowerShear(h_ref=100, alpha=.2))\n",
"\n",
"# Site with wind direction dependent weibull distributed wind speed\n",
"uniform_weibull_site = XRSite(\n",
" ds=xr.Dataset(data_vars={'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti},\n",
" coords={'wd': wd}))\n",
"\n",
"# Site with a speedup and a turning value per WT\n",
"x_i, y_i = np.arange(5) * 100, np.zeros(5) # WT positions\n",
"\n",
"complex_fixed_pos_site = XRSite(\n",
" ds=xr.Dataset(\n",
" data_vars={'Speedup': ('i', np.arange(.8, 1.3, .1)),\n",
" 'Turning': ('i', np.arange(-2, 3)),\n",
" 'P': ('wd', f)},\n",
" coords={'i': np.arange(5), 'wd': wd}),\n",
" initial_position=np.array([x_i, y_i]).T)\n",
"\n",
"# Site with gridded speedup information\n",
"complex_grid_site = XRSite(\n",
" ds=xr.Dataset(\n",
" data_vars={'Speedup': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))),\n",
" 'P': ('wd', f)},\n",
" coords={'x': [0, 500, 1000], 'y': [0, 500], 'wd': wd}))\n",
"\n",
"# Site with ws dependent speedup and wd- and ws distributed probability\n",
"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)\n",
"P_wd_ws = P_ws[na, :] * np.array(f)[:, na]\n",
"\n",
"complex_ws_site = XRSite(\n",
" ds=xr.Dataset(\n",
" data_vars={'Speedup': (['ws'], np.arange(.8, 1.4, .1)),\n",
" 'P': (('wd', 'ws'), P_wd_ws), 'TI': ti},\n",
" coords={'ws': [1.5, 4, 6, 8, 10, 12], 'wd': wd}))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Wake effects from neighbouring wind farms\n",
"\n",
"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.\n",
"\n",
"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.\n",
"\n",
"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). "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# import and setup site and windTurbines\n",
"from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines\n",
"from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit\n",
"from py_wake.wind_turbines import WindTurbine, WindTurbines\n",
"from py_wake.wind_farm_models import PropagateDownwind\n",
"from py_wake.superposition_models import LinearSum\n",
"\n",
"site = IEA37Site(16)\n",
"\n",
"# setup current, neighbour and all positions\n",
"wt_x, wt_y = site.initial_position.T\n",
"neighbour_x, neighbour_y = wt_x-4000, wt_y\n",
"all_x, all_y = np.r_[wt_x,neighbour_x], np.r_[wt_y,neighbour_y]\n",
"\n",
"windTurbines = WindTurbines.from_WindTurbine_lst([IEA37_WindTurbines(),IEA37_WindTurbines()])\n",
"windTurbines._names = [\"Current wind farm\",\"Neighbour wind farm\"]\n",
"types = [0]*len(wt_x) + [1]*len(neighbour_x)\n",
"\n",
"wf_model = PropagateDownwind(site, windTurbines,\n",
" wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True),\n",
" superpositionModel=LinearSum())\n",
"\n",
"# Consider wd=270 +/- 30 deg only\n",
"wd_lst = np.arange(240,301)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#plotting the wake maps for the desired flow case\n",
"wsp = 9.8\n",
"wdir = 267\n",
"\n",
"plt.figure(figsize=(16, 6))\n",
"wf_model(all_x, all_y, type=types, wd=wdir, ws=wsp, h=110).flow_map().plot_wake_map()\n",
"plt.xlabel('x [m]')\n",
"plt.ylabel('y [m]')\n",
"plt.title('Wake map for'+ f' {wdir} deg and {wsp} m/s')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Now, we run the simulation of all wind turbines and calculate AEP of current wind farm**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**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**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#making a flow box covering the area of interest (i.e the current wind farm + 100m)\n",
"\n",
"ext = 1000\n",
"flow_box = wf_model(neighbour_x, neighbour_y, wd=wd_lst).flow_box(\n",
" x=np.linspace(min(wt_x) - ext, max(wt_x) + ext, 101),\n",
" y=np.linspace(min(wt_y) - ext, max(wt_y) + ext, 101),\n",
" h=110)\n",
"\n",
"#creating new site based on the flow box\n",
"\n",
"from py_wake.site.xrsite import XRSite\n",
"wake_site = XRSite.from_flow_box(flow_box)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we plot the \"free-stream\" inflow wind speed of the current wind farm."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(16, 6))\n",
"wake_site.ds.WS.sel(wd=267).plot(y='y', cmap = 'Blues_r')\n",
"windTurbines.plot(all_x, all_y, types, wd=270)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, we setup a new wind farm model with the new pre-generated site and calculate the AEP."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"wf_model_wake_site = PropagateDownwind(wake_site, windTurbines,\n",
" wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True),\n",
" superpositionModel=LinearSum())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"Total AEP: %f GWh\"%wf_model_wake_site(wt_x, wt_y, ws=[wsp], wd=wd_lst).aep().sum())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"**Lastly, we plot the flow map of the current wind farm with the selected flow box.**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(16, 6))\n",
"wf_model_wake_site(wt_x, wt_y, wd=wdir, ws=wsp, h=110).flow_map().plot_wake_map()\n",
"windTurbines.plot(neighbour_x, neighbour_y, type=1, wd=wdir)\n",
"plt.xlabel('x [m]')\n",
"plt.ylabel('y [m]')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local wind\n",
"The method `local_wind` is used to calculate the local wind in a wind farm given certain turbine positions or coordinates. The class returns a `LocalWind`-dictionary."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"localWinds = {name: site.local_wind(x=site.initial_position[:,0], # x position\n",
" y = site.initial_position[:,1], # y position\n",
" h=site.initial_position[:,0]*0+70, # height\n",
" ws=None, # defaults to 3,4,..,25\n",
" wd=None, # defaults to 0,1,...,360\n",
" ) for name, site in sites.items()}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`LocalWind.coords` contains the current coordinates, e.g.:\n",
"\n",
"- i: Point number. Points can be wind turbine position or just points in a flow map\n",
"- wd: Ambient reference wind direction\n",
"- ws: Ambient reference wind speed\n",
"- x,y,h: position and height of points\n",
"\n",
"while the dictionary itself contains some data variables:\n",
"\n",
"- WD: Local wind direction\n",
"- WS: Local wind speed\n",
"- TI: Local turbulence intensity\n",
"- P: Probability of flow case (wind direction and wind speed)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `IEA37` site has 16 wind turbines on a uniform site with a fixed wind speed of 9.8 m/s and the data variables therefore only depend on wind direction."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print (localWinds['IEA37'].coords.keys())\n",
"localWinds['IEA37'].P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `Hornsrev1` site has 80 wind turbines on a uniform site and the data variables therefore depend on wind direction and wind speed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"localWinds['Hornsrev1'].P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the `ParqueFicticio` site has 8 turbines in a complex terrain and the data variables therefore depend on wind direction, wind speed, and position."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"localWinds['ParqueFicticio'].P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Wind speeds at the wind turbines for reference wind speed of 3m/s (k=0):**\n",
"\n",
"- `IEA37`: Constant wind speed of **9.8m/s**\n",
"- `Hornsrev1`: Constant wind speed over the site, **3 m/s**\n",
"- `ParqueFicticio`: Winds speed depends on both wind direction and position"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for name, lw in localWinds.items():\n",
" print (name)\n",
" print (lw.WS.values, 'm/s')\n",
" print (\"=\"*100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The ParqueFicticio site models variations within the site, so the local wind speed varies over the area."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"s = sites[\"ParqueFicticio\"]\n",
"x = np.linspace(262878,264778,300)\n",
"y = np.linspace(6504714,6506614,300)\n",
"X,Y = np.meshgrid(x,y)\n",
"lw = s.local_wind(X.flatten(),Y.flatten(),30, ws=[10],wd=[0])\n",
"Z = lw.WS_ilk.reshape(X.shape)\n",
"c = plt.contourf(X,Y,Z, levels=100)\n",
"plt.colorbar(c,label='Wind speed [m/s]')\n",
"plt.title(\"Local wind speed at 10m/s and 0deg\")\n",
"plt.xlabel('x [m]')\n",
"plt.ylabel('y [m]')\n",
"plt.axis('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Distance\n",
"**We can also calculate the distance between points of a specific site for either flat or complex terrain.**\n",
"\n",
"For the `IEA37Site` and the `Hornsrev1` sites the distances between points are straight line distances, as these sites are characterized by flat terrain.\n",
"\n",
"For the `ParqueFicticioSite`, on the other hand, the down-wind distance is larger as it follows the non-flat terrain."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"wd = [0, 30,90] # wind direction at source\n",
"\n",
"for name, site in sites.items():\n",
" print (\"------- %s -------\"%name)\n",
" wt_x, wt_y = site.initial_position[0]\n",
" site.distance.setup(src_x_ilk=[wt_x, wt_x], src_y_ilk=[wt_y, wt_y-1000], src_h_ilk=[70,90]) # wt2 1000m to the south\n",
" dw_ijlk, cw_ijlk, dh_ijlk = site.distance(wd_l=wd, src_idx=[0], dst_idx=[[1,1,1]])\n",
" \n",
"\n",
" print ('Wind direction: \\t\\t%d deg\\t\\t%d deg\\t\\t%d deg'%tuple(wd))\n",
" print ('Down wind distance [m]: \\t%.1f\\t\\t%.1f\\t\\t%.1f'%tuple(dw_ijlk[0,0,:,0]))\n",
" print ('Cross wind distance [m]: \\t%.1f\\t\\t%.1f\\t\\t%.1f'%tuple(cw_ijlk[0,0,:,0]))\n",
" print ('Height difference [m]: \\t\\t%.1f\\t\\t%.1f\\t\\t%.1f'%tuple(dh_ijlk[0,0,:,0]))\n",
" print()\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Wind resource distribution plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"site = sites['Hornsrev1']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting wind rose."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"_ = site.plot_wd_distribution(n_wd=12, ws_bins=[0,5,10,15,20,25])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting probability density function for the four sectors studied."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"_ = site.plot_ws_distribution(wd=[0,90,180,270])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting probablity density function for the four sector studied $\\pm$ 45 degrees.\n",
"\n",
"If **include_wd_distribution=true**, the wind speed probability distributions are multiplied by the wind direction probability.\n",
"\n",
"The sector size is set to 360 / len(wd). This only makes sense if the wd array is evenly distributed"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"_ = site.plot_ws_distribution(wd=[0,90,180,270], include_wd_distribution=True)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.13"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}