Update IEA37SimpleBastankhahGaussian.ipynb

d136d6d8 · Mads M. Pedersen · 9906e710 · d136d6d8
Commit d136d6d8 authored 6 years ago by Mads M. Pedersen
--- a/_notebooks/elements/IEA37SimpleBastankhahGaussian.ipynb
+++ b/_notebooks/elements/IEA37SimpleBastankhahGaussian.ipynb
@@ -4,7 +4,7 @@
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "# Bastankhah Gaussian WakeModel"
+    "# IEA37 Simple Bastankhah Gaussian WakeModel"
   ]
  },
  {

 %% Cell type:markdown id: tags:

-# Bastankhah Gaussian WakeModel
+# IEA37 Simple Bastankhah Gaussian WakeModel

 %% Cell type:code id: tags:

 ``` python
 import numpy as np
 import matplotlib.pyplot as plt

 from py_wake.examples.data.hornsrev1 import V80, Hornsrev1Site, wt9_x, wt9_y
 from py_wake.wake_models import IEA37SimpleBastankhahGaussian


 windTurbines = V80()
 wake_model = IEA37SimpleBastankhahGaussian(windTurbines)
 ```

 %% Cell type:markdown id: tags:

 The `IEA37SimpleBastankhahGaussian` is a subclass of the general `WakeModel` class, see documentation [here](https://topfarm.pages.windenergy.dtu.dk/PyWake/wake_models/WakeModel.html)

 Implemented according to [IEA37](https://github.com/byuflowlab/iea37-wflo-casestudies/blob/master/iea37-wakemodel.pdf) and is equivalent to `BastankhahGaussian` for $beta=1/\sqrt{8} \sim ct=0.9637188$

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

 - i: turbines ordered by id
 - k: wind speeds
 - l: wind directions

 This means that `WS_ilk[0,1,2]` holds the wind speed at the first turbine for the second wind direction and third wind speed

 %% Cell type:markdown id: tags:

 `WakeModel` contains a method, [`calc_wake`](https://topfarm.pages.windenergy.dtu.dk/PyWake/wake_models/WakeModel.html#py_wake.wake_model.WakeModel.calc_wake), to calculate the effective wind speed, turbulence intensity (not implemented yet), power and thrust coefficient.

 Let us try to calculate the effective wind speed for two V80 turbines separated by 200m in 10m/s and wind direction parallel to a line between the two turbines

 %% Cell type:code id: tags:

 ``` python
 WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk = wake_model.calc_wake(
    WS_ilk=np.array([[[10]],[[10]]]), # 10m/s at both turbines
    TI_ilk=np.array([[[.1]],[[.1]]]), # turbulence intensity (not used)
    dw_iil=np.array([[0,200],[-200,0]])[:,:,np.newaxis], # 200m down stream separation
    hcw_iil=np.array([[0,0],[0,0]])[:,:,np.newaxis], # wt aligned with wind -> zero cross wind distance
    dh_iil=np.array([[0,0],[0,0]])[:,:,np.newaxis], # no hub height difference
    dw_order_indices_dl=np.array([[0,1]]), # down stream order of turbines
    types_i=[0,0]) # both are turbine type 0

 for i in [0,1]:
    print ('Turbine', i)
    print ('Effective wind speed %fm/s'%WS_eff_ilk[i,0,0])
    print ('Power production %.2fW'%power_ilk[i,0,0])
    print ('Trust coefficient %f'%ct_ilk[i,0,0])
    print()
 ```

 %% Output

    Turbine 0
    Effective wind speed 10.000000m/s
    Power production 1341000.00W
    Trust coefficient 0.793000
    
    Turbine 1
    Effective wind speed 6.895002m/s
    Power production 441310.28W
    Trust coefficient 0.804895
    

 %% Cell type:markdown id: tags:

 To calculate this, `calc_wake`, uses two wake-model specific methods, `calc_deficit` and `calc_effective_WS`.

 `calc_deficit` calculates the deficit:

 %% Cell type:code id: tags:

 ``` python
 deficit = wake_model.calc_deficit(
    WS_lk=np.array([[10]]), # wind speed at current turbine
    D_src_l=np.array([80]), # diameter of current turbine
    dw_jl=np.array([[200]]), # down wind distance
    cw_jl=np.array([[0]]), # cross wind distance (both horizontal and vertical)
    ct_lk=np.array([[0.793000]])) # thrust coefficient

 print (deficit[0,0,0])
 ```

 %% Output

    3.1049984549612297

 %% Cell type:code id: tags:

 ``` python
 x = np.arange(-300,301,1.)
 y = np.arange(-100,1001,1)
 X,Y = np.meshgrid(x,y)
 down_stream = Y>0
 deficit_map=np.zeros_like(X)

 deficit_map[down_stream] = wake_model.calc_deficit(
    WS_lk=np.array([[10]]), # wind speed at current turbine
    D_src_l=np.array([80]), # diameter of current turbine
    dw_jl=Y[down_stream][:,np.newaxis], # down wind distance
    cw_jl=np.abs(X[down_stream])[:,np.newaxis], # cross wind distance (both horizontal and vertical)
    ct_lk=np.array([[0.793000]]))[:,0,0] # thrust coefficient

 c = plt.contourf(Y,X,deficit_map,100)
 plt.colorbar(c)
 plt.plot([0,0],[-40,40],'r',lw=3)
 ```

 %% Output

    [<matplotlib.lines.Line2D at 0x1e5b118d780>]



 %% Cell type:markdown id: tags:

 while `calc_effective_WS` calculates the effective wind speed by subtracting the deficits from all upstream turbines from the local wind speed. For `NOJ` it subtracts the square root of the sum of squared deficits

 %% Cell type:code id: tags:

 ``` python
 wake_model.calc_effective_WS(
    WS_lk=np.array([[10]]),
    deficit_jlk=deficit)
 ```

 %% Output

    array([[6.89500155]])

 %% Cell type:markdown id: tags:

 Finally, `WakeModel`, contains the method `wake_map` to find the effective wind speed at arbitrary positions

 %% Cell type:code id: tags:

 ``` python
 #calculate the wake 200m down stream of a V80

 wake_model.wake_map(
    WS_ilk=np.array([[[10]]]), # wind speed at turbine
    WS_eff_ilk=np.array([[[8.92340252]]]),
    dw_ijl=np.array([[[200]]]), # one point 500m down stream
    hcw_ijl=np.array([[[0]]]), # 0m cross wind distance
    dh_ijl=np.array([[[0]]]), # at hub height
    ct_ilk=np.array([[[0.793000]]]),
    types_i=[0],
    WS_jlk=np.array([[[10]]]), # local wind speed at point
 )
 ```

 %% Output

    array([[[6.89500155]]])

 %% Cell type:markdown id: tags:

 For standard purposes, however, we do not call these methods manually. Instead we use the functions in `AEPCalculator`

 %% Cell type:markdown id: tags:

 **Calcculate AEP**

 %% Cell type:code id: tags:

 ``` python
 from py_wake.aep_calculator import AEPCalculator

 site = Hornsrev1Site()
 aep_calc = AEPCalculator(site, windTurbines, wake_model)

 aep_gwh_ilk = aep_calc.calculate_AEP(x_i=[0,0], y_i=[0,-200])
 aep_gwh_noloss_ilk = aep_calc.calculate_AEP_no_wake_loss(x_i=[0,0], y_i=[0,-200])


 # AEP pr turbine
 print ('AEP pr turbine:', aep_gwh_ilk.sum((1,2)))

 # AEP pr wind direction
 plt.plot(aep_gwh_noloss_ilk[:,:,7].sum(0), label='Without wake loss')
 plt.plot(aep_gwh_ilk[:,:,7].sum(0), label='With wake loss')
 plt.title('Wind speed: 10m/s')
 plt.xlabel('Wind direction [deg]')
 plt.ylabel('AEP [GWh]')
 plt.legend()

 # AEP pr wind speed
 plt.figure()
 plt.plot(site.default_ws, aep_gwh_noloss_ilk[:,0].sum(0), label='Without wake loss')
 plt.plot(site.default_ws, aep_gwh_ilk[:,0].sum(0), label='With wake loss')
 plt.title('Wind direction: 0deg')
 plt.xlabel('Wind speed [deg]')
 plt.ylabel('AEP w[GWh]')
 plt.legend()
 ```

 %% Output

    AEP pr turbine: [9.01126077 9.18331398]

    <matplotlib.legend.Legend at 0x1e5b1b14f28>





 %% Cell type:markdown id: tags:

 **Wake map**

 %% Cell type:code id: tags:

 ``` python
 import matplotlib.pyplot as plt
 x,y = wt9_x, wt9_y
 aep_calc.plot_wake_map(wt_x=x, wt_y=y, wd=[0], ws=[10])
 windTurbines.plot(x, y)
 plt.show()
 ```

 %% Output



 %% Cell type:markdown id: tags:

 **Effective wind speed, power and thrust coefficient**

 %% Cell type:code id: tags:

 ``` python
 # wind from 0 deg(index=0) and 10m/s (index=7)

 aep_calc.calculate_AEP(x_i=[0,0], y_i=[0,-200])
 print ("Effective wind speed: wt0: %f\twt1: %f"%tuple(aep_calc.WS_eff_ilk[:,0,7]))
 print ("Power: wt0: %dW\twt1: %dW"%tuple(aep_calc.power_ilk[:,0,7]))
 print ("Thrust coefficient: wt0: %dW\twt1: %dW"%tuple(aep_calc.ct_ilk[:,0,7]))
 print ("Probability of ws=10m/s and wd=0deg: %f"%aep_calc.P_lk[0,7])
 ```

 %% Output

    Effective wind speed: wt0: 10.000000	wt1: 6.895002
    Power: wt0: 1341000W	wt1: 441310W
    Thrust coefficient: wt0: 0W	wt1: 0W
    Probability of ws=10m/s and wd=0deg: 0.000103

 %% Cell type:code id: tags:

 ``` python
 ```