Open and run in Colab (interactive) Edit on Gitlab

Ørsted TurbOPark

This notebook reproduces the results and figures of the two examples provided by Ørsted in https://github.com/OrstedRD/TurbOPark/blob/main/TurbOParkExamples.mlx

The PyWake look-up table used for the GaussianOverlap model is slightly different from look-up table used by Ørsted (which has a finer grid resolution). The results are therefore slightly different (relative error < 1e-5).

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

Setting up the wind turbine and site objects

[2]:

import numpy as np
import matplotlib.pyplot as plt

from py_wake.wind_turbines._wind_turbines import WindTurbine, WindTurbines
from py_wake.wind_turbines.power_ct_functions import PowerCtTabular
from py_wake.site._site import UniformSite
from py_wake.site.shear import PowerShear
from py_wake.utils.plotting import setup_plot

Specifying the wind speed and power and CT curve

[3]:

u = np.arange(0, 25.5, .5)
po = [0, 0, 0, 0, 5, 15, 37, 73, 122, 183, 259, 357, 477, 622, 791, 988, 1212, 1469, 1755, 2009, 2176, 2298, 2388, 2447, 2485, 2500, 2500, 2500,
      2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500, 2500,
      2500, 2500, 2500, 2500, 2500, 2500, 2500, 0]
ct = [0, 0, 0, 0, 0.78, 0.77, 0.78, 0.78, 0.77, 0.77, 0.78, 0.78, 0.78, 0.78, 0.78, 0.78, 0.77, 0.77, 0.77, 0.76, 0.73, 0.7, 0.68, 0.52, 0.42,
      0.36, 0.31, 0.27, 0.24, 0.22, 0.19, 0.18, 0.16, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09, 0.08, 0.08, 0.08, 0.07, 0.07, 0.06, 0.06, 0.06,
      0.05, 0.05, 0.05, 0.04, 0]

wt1 = WindTurbine(name="Ørsted1", diameter=120, hub_height=100, powerCtFunction=PowerCtTabular(u, po, 'kw', ct))

u2 = np.arange(0, 27)
pow2 = [0, 0, 0, 0, 54, 144, 289, 474, 730, 1050, 1417, 1780, 2041, 2199, 2260, 2292, 2299, 2300, 2300, 2300, 2300, 2300, 2300, 2300, 2300, 2300, 0]
ct2  = [0, 0, 0, 0, 0.94, 0.82, 0.76, 0.68, 0.86, 0.83, 0.77, 0.68, 0.66, 0.52, 0.47, 0.41, 0.38, 0.34, 0.27, 0.26, 0.23, 0.22, 0.22, 0.2, 0.16, 0.17, 0]
wt2 = WindTurbine(name="Ørsted2", diameter=80, hub_height=70, powerCtFunction=PowerCtTabular(u2, pow2, 'kw', ct2))
wts = WindTurbines.from_WindTurbine_lst([wt1,wt2])

Plotting the power curve

[4]:

u = np.linspace(0,25)

for t in [0,1]:
    plt.plot(u,wts.power(u, type=t)/1000, label=wts.name(t))
setup_plot(xlabel='Wind speed [m/s]', ylabel='Power [kW]', title='Power curve')
../../_images/notebooks_literature_verification_TurbOPark_8_0.png

Plotting the CT curve

[5]:

for t in [0,1]:
    plt.plot(u,wts.ct(u, type=t), label=wts.name(t))
setup_plot(xlabel='Wind speed [m/s]', ylabel='Thrust coefficient [-]', title='Thrust curve')
../../_images/notebooks_literature_verification_TurbOPark_10_0.png

Wind speed, wind direction and turbulence intensity

[6]:

u0 = [6,10,14]           # [m/s]
wd = 270                 # [deg]
ti0 = [0.09,.1,.11]      # [-]

Example 1 - Square farm with identical turbines

[7]:

y, x = [v.flatten() for v in np.meshgrid(np.arange(4) * 120 * 6, np.arange(4) * 120 * 6)]
wt1.plot(x, y)
../../_images/notebooks_literature_verification_TurbOPark_14_0.png 
[8]:

site = UniformSite(shear=PowerShear(h_ref=90, alpha=.1))

[9]:

from py_wake.literature import Nygaard_2022

wfm = Nygaard_2022(site, wt1)
sim_res = wfm(x, y, ws=u0, wd=wd, TI=ti0)
sim_res.WS_eff

[9]:

<xarray.DataArray 'WS_eff' (wt: 16, wd: 1, ws: 3)>
array([[[ 6.06355051, 10.10591751, 14.14828452]],

       [[ 6.06355051, 10.10591751, 14.14828452]],

       [[ 6.06355051, 10.10591751, 14.14828452]],

       [[ 6.06355051, 10.10591751, 14.14828452]],

       [[ 4.35804581,  7.37377436, 12.8635598 ]],

       [[ 4.35804581,  7.37377436, 12.8635598 ]],

       [[ 4.35804581,  7.37377436, 12.8635598 ]],

       [[ 4.35804581,  7.37377436, 12.8635598 ]],

       [[ 3.82923427,  6.45011265, 12.02767091]],

       [[ 3.82923427,  6.45011265, 12.02767091]],

       [[ 3.82923427,  6.45011265, 12.02767091]],

       [[ 3.82923427,  6.45011265, 12.02767091]],

       [[ 3.48367104,  5.87190035, 11.17215029]],

       [[ 3.48367104,  5.87190035, 11.17215029]],

       [[ 3.48367104,  5.87190035, 11.17215029]],

       [[ 3.48367104,  5.87190035, 11.17215029]]])
Coordinates:
  * ws       (ws) int64 6 10 14
  * wt       (wt) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
  * wd       (wd) int64 270
    type     (wt) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
Attributes:
    Description:  Effective local wind speed [m/s]

We now calculate the power of each turbine in the wind farm

[10]:

sim_res.Power /1000

[10]:

<xarray.DataArray 'Power' (wt: 16, wd: 1, ws: 3)>
array([[[ 495.42964709, 2201.84387294, 2500.        ]],

       [[ 495.42964709, 2201.84387294, 2500.        ]],

       [[ 495.42964709, 2201.84387294, 2500.        ]],

       [[ 495.42964709, 2201.84387294, 2500.        ]],

       [[ 165.6815891 ,  938.26709823, 2500.        ]],

       [[ 165.6815891 ,  938.26709823, 2500.        ]],

       [[ 165.6815891 ,  938.26709823, 2500.        ]],

       [[ 165.6815891 ,  938.26709823, 2500.        ]],

       [[ 105.26495844,  607.53266956, 2485.8301274 ]],

       [[ 105.26495844,  607.53266956, 2485.8301274 ]],

       [[ 105.26495844,  607.53266956, 2485.8301274 ]],

       [[ 105.26495844,  607.53266956, 2485.8301274 ]],

       [[  71.82431459,  446.25608366, 2408.31373411]],

       [[  71.82431459,  446.25608366, 2408.31373411]],

       [[  71.82431459,  446.25608366, 2408.31373411]],

       [[  71.82431459,  446.25608366, 2408.31373411]]])
Coordinates:
  * ws       (ws) int64 6 10 14
  * wt       (wt) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
  * wd       (wd) int64 270
    type     (wt) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0

Then, we calculate the wake losses

[11]:

((1-(sim_res.Power.mean('wt') / sim_res.Power.max('wt')))*100).plot()
setup_plot(ylabel='Wake loss [%]', xlabel='Freestream wind speed [m/s]')
../../_images/notebooks_literature_verification_TurbOPark_20_0.png

Example 2 - Two turbine types plus wind speed gradient

First we set up the two types of wind turbines in the farm

[12]:

type = np.array([1,1,0,0]*4)
wts.plot(x, y, type=type)
../../_images/notebooks_literature_verification_TurbOPark_23_0.png 
[13]:

import xarray as xr
from py_wake.site import XRSite

x_pt = np.arange(-100,3100,100)
Y_pt,X_pt = np.meshgrid(x_pt,x_pt)
grad = ((X_pt-5)**2 + (Y_pt)**2)*10**(-8) + 1
speedup= grad/grad[1,1]
ds = xr.Dataset({'Speedup':(('x','y'), speedup), 'P':1}, coords={'x':x_pt,'y':x_pt})
ds.Speedup.plot()
wts.plot(x, y, type=type)
../../_images/notebooks_literature_verification_TurbOPark_24_0.png

We now specify the gradients of the wind speed

[14]:

gradient_site = XRSite(ds=ds, shear=PowerShear(h_ref=90, alpha=.1))
gradient_site.local_wind(x,y,wts.hub_height(type), ws=1)['WS_ilk']

[14]:

array([[[0.97518172]],

       [[0.98025267]],

       [[1.03157163]],

       [[1.05776616]],

       [[0.98018245]],

       [[0.9852534 ]],

       [[1.03675394]],

       [[1.06294848]],

       [[0.99528606]],

       [[1.00035701]],

       [[1.05240599]],

       [[1.07860052]],

       [[1.02049256]],

       [[1.0255635 ]],

       [[1.07852775]],

       [[1.10472229]]])

[15]:

from py_wake.literature import Nygaard_2022

wfm = Nygaard_2022(gradient_site, wts)
sim_res = wfm(x, y, ws=u0, wd=wd, type=type, TI=ti0)
sim_res.WS_eff

[15]:

<xarray.DataArray 'WS_eff' (wt: 16, wd: 1, ws: 3)>
array([[[ 5.85109034,  9.75181723, 13.65254412]],

       [[ 5.881516  ,  9.80252667, 13.72353733]],

       [[ 6.18942978, 10.31571631, 14.44200283]],

       [[ 6.34659697, 10.57766162, 14.80872627]],

       [[ 4.37635315,  7.29654649, 11.30769686]],

       [[ 4.40032994,  7.33705588, 11.38300503]],

       [[ 4.47086687,  7.59526343, 13.26337806]],

       [[ 4.58382739,  7.82538149, 13.71972616]],

       [[ 4.0259611 ,  6.71637847, 10.09987749]],

       [[ 4.04256332,  6.7422665 , 10.162897  ]],

       [[ 3.98767263,  6.73423506, 12.72628796]],

       [[ 4.08583929,  6.92565811, 13.26102546]],

       [[ 3.94419704,  6.48890696,  9.45282588]],

       [[ 3.95541264,  6.51954795,  9.5148913 ]],

       [[ 3.71960601,  6.2782012 , 12.39637986]],

       [[ 3.80940787,  6.44644865, 13.04456244]]])
Coordinates:
  * ws       (ws) int64 6 10 14
  * wt       (wt) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
  * wd       (wd) int64 270
    type     (wt) float64 1.0 1.0 0.0 0.0 1.0 1.0 ... 0.0 0.0 1.0 1.0 0.0 0.0
Attributes:
    Description:  Effective local wind speed [m/s]

Now we can again calculate the power for the new wind farm configuration

[16]:

sim_res.Power / 1000

[16]:

<xarray.DataArray 'Power' (wt: 16, wd: 1, ws: 3)>
array([[[ 267.40809898, 1325.91692326, 2238.80519141]],

       [[ 271.81981999, 1344.52728661, 2243.13577732]],

       [[ 531.93463744, 2253.03477894, 2500.        ]],

       [[ 577.5131225 , 2311.97909224, 2500.        ]],

       [[  87.87178378,  549.9159009 , 1860.30888155]],

       [[  90.02969426,  560.28630549, 1879.96431321]],

       [[ 179.44575863, 1030.67801758, 2500.        ]],

       [[ 195.74176377, 1133.77090798, 2500.        ]],

       [[  56.33649862,  421.5300164 , 1453.25552871]],

       [[  57.8306992 ,  426.31930299, 1476.13161259]],

       [[ 120.79191747,  701.17144884, 2500.        ]],

       [[ 132.47239357,  765.8724426 , 2500.        ]],

       [[  50.98664016,  379.44778803, 1216.18709616]],

       [[  51.59228249,  385.11637136, 1238.96510883]],

       [[  94.52138919,  557.67834759, 2496.89139586]],

       [[ 103.3219714 ,  606.47010743, 2500.        ]]])
Coordinates:
  * ws       (ws) int64 6 10 14
  * wt       (wt) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
  * wd       (wd) int64 270
    type     (wt) float64 1.0 1.0 0.0 0.0 1.0 1.0 ... 0.0 0.0 1.0 1.0 0.0 0.0
[17]:

#turbulence intensity
sim_res.TI

[17]:

<xarray.DataArray 'TI' (ws: 3)>
array([0.09, 0.1 , 0.11])
Coordinates:
  * ws       (ws) int64 6 10 14
Attributes:
    Description:  
[18]:

for ws in u0:
    plt.figure()
    for row, p in enumerate((sim_res.Power.sel(ws=ws, wd=270).values.reshape((4,4))).T/1000,1):
        plt.plot(sim_res.x[::4], p, label=f'Row {row}')
    plt.plot(sim_res.x.sel(wt=15), sim_res.Power.sel(ws=ws, wd=270,wt=15)/1000,'.r', 'WT15')
    setup_plot(title=f'Wind speed {ws}m/s', ylabel='Power [kW]', xlabel='X coordinate [m]')
../../_images/notebooks_literature_verification_TurbOPark_31_0.png
../../_images/notebooks_literature_verification_TurbOPark_31_1.png
../../_images/notebooks_literature_verification_TurbOPark_31_2.png