Blockage modelling 1¶

Wind farm blockage¶

This example shows how to use the Rankine-Half-Body model for modelling the blockage effect. We need the following imports:

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import foxes
import foxes.variables as FV
import foxes.constants as FC

First, we setup the model book and the wind farm, here we pick 20 x 10 turbines on a regular grid:

In [2]:

farm = foxes.WindFarm()
foxes.input.farm_layout.add_grid(
    farm=farm,
    xy_base=[0.0, 0.0],
    step_vectors=[[800.0, 50.0], [0., 400.0]],
    steps=[6, 8],
    turbine_models=["DTU10MW"],
    verbosity=0
)

ax = foxes.output.FarmLayoutOutput(farm).get_figure(figsize=(8, 4))
plt.show()

../_images/notebooks_blockage_farm_4_0.png

As wind conditions, we choose a single wind direction state from 270 degrees:

In [3]:

# create ambient wind conditions, a single uniform state:
states = foxes.input.states.SingleStateStates(ws=8.0, wd=270.0, ti=0.04, rho=1.225)

For induction models, the Iterative algorithm should be invoked. This is due to the fact that upstream turbines can be influenced by the induction zone of a turbine, which the Downstream algorithm simply would not notice.

By default the Iterative algorithm adds a smooth cutin region to the thrust coefficient curve. This behaviour can be changed by modifying the default choice of the constructor argument mod_cutin={modify_ct=True, modify_P=False}. Alternatively, under-relaxation of the variable FV.CT could by applied, e.g. by calling aglo.set_urelax(pre_wake={FV.CT: 0.5}) before calculations. This means that for each iteration the freshly calculated value is mixed with the result from the previous iteration. Both approaches avoid wakes switching on and off repeatedly between subsequent iterations in some cases, ensuring convergence also in the vicinity of the cut-in and cut-out regions of the thrust curve. The convergence criteria can also be modified by setting the parameter conv_crit, see API.

In [4]:

algo = foxes.algorithms.Iterative(
    farm,
    states,
    rotor_model="centre",
    wake_models=["RHB", "Bastankhah025_linear_lim_k004"],
    wake_frame="rotor_wd",
    partial_wakes="grid36",
    chunks={FC.STATE: None, FC.POINT: 4000},
    verbosity=1,
)

Notice the RHB model among the list of wake models. This is the Rankine-Half-Body induction model, adding an induction zone for each turbine. We can now calculate the results:

In [5]:

with foxes.utils.runners.DaskRunner() as runner:
    farm_results = runner.run(algo.calc_farm)

print(farm_results)


Algorithm Iterative: Iteration 0


Initializing algorithm 'Iterative'
Initializing model 'Iterative'
Initializing model 'SingleStateStates'
Initializing model 'centre'
Initializing model 'basic_ctrl_prer'
Initializing model 'DTU10MW'
Turbine type 'DTU10MW': Reading static data from context 'power_ct_curve'
Path: /home/jonas/gits/wakes/foxes/foxes/data/power_ct_curves/DTU-10MW-D178d3-H119.csv
Initializing model 'basic_ctrl_postr'
Initializing model 'basic_ctrl'
Initializing model 'rotor_wd'
Initializing model 'Madsen'
Initializing model 'RHB'
Initializing model 'WakeK'
Initializing model 'Bastankhah025_linear_lim_k004'
Initializing model 'GridRotor'
Initializing model 'grid36'

--------------------------------------------------
  Running Iterative: calc_farm
--------------------------------------------------
  n_states : 1
  n_turbines: 48
--------------------------------------------------
  states   : SingleStateStates()
  rotor    : CentreRotor()
  controller: BasicFarmController()
  wake frame: RotorWD()
--------------------------------------------------
  wakes:
    0) RHB: RankineHalfBody(induction=Madsen)
    1) Bastankhah025_linear_lim_k004: Bastankhah2014(ws_linear_lim, induction=Madsen, k=0.04)
--------------------------------------------------
  partial wakes:
    0) RHB: grid36, PartialGrid(rotor_model=GridRotor)
    1) Bastankhah025_linear_lim_k004: grid36, PartialGrid(rotor_model=GridRotor)
--------------------------------------------------
  turbine models:
    0) DTU10MW: PCtFile(D=178.3, H=119.0, P_nominal=10000.0, P_unit=kW, rho=1.225, var_ws_ct=REWS2, var_ws_P=REWS3)
--------------------------------------------------

Initializing model 'InitFarmData'
Initializing model 'SetAmbFarmResults'
Initializing model 'FarmWakesCalculation'
Initializing model 'Iterative_calc'

--------------------------------------------------
  Model oder
--------------------------------------------------
00) basic_ctrl
01) InitFarmData
02) centre
03) basic_ctrl
  03.0) Post-rotor: DTU10MW
04) SetAmbFarmResults
05) FarmWakesCalculation
--------------------------------------------------

[########################################] | 100% Completed | 101.67 ms

Input data:

 <xarray.Dataset> Size: 468B
Dimensions:      (state: 1, turbine: 48, tmodels: 1)
Coordinates:
  * state        (state) int64 8B 0
  * tmodels      (tmodels) <U7 28B 'DTU10MW'
Dimensions without coordinates: turbine
Data variables:
    weight       (state, turbine) float64 384B dask.array<chunksize=(1, 48), meta=np.ndarray>
    tmodel_sels  (state, turbine, tmodels) bool 48B dask.array<chunksize=(1, 48, 1), meta=np.ndarray>


Farm variables: AMB_CT, AMB_P, AMB_REWS, AMB_REWS2, AMB_REWS3, AMB_RHO, AMB_TI, AMB_WD, AMB_YAW, CT, D, H, P, REWS, REWS2, REWS3, RHO, TI, WD, X, Y, YAW, order, order_inv, order_ssel, weight

Output variables: AMB_CT, AMB_P, AMB_REWS, AMB_REWS2, AMB_REWS3, AMB_RHO, AMB_TI, AMB_WD, AMB_YAW, CT, D, H, P, REWS, REWS2, REWS3, RHO, TI, WD, X, Y, YAW, order, order_inv, order_ssel, weight

Chunks: {'state': None, 'point': 4000, 'target': 4000}


Calculating 1 states for 48 turbines
[########################################] | 100% Completed | 103.39 ms

Algorithm Iterative: Iteration 1

[########################################] | 100% Completed | 101.50 ms
[########################################] | 100% Completed | 101.61 ms

DefaultConv: Convergence check
  REWS: delta = 7.281e-01, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 1.153e-02, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 2

[########################################] | 100% Completed | 101.29 ms
[########################################] | 100% Completed | 104.23 ms

DefaultConv: Convergence check
  REWS: delta = 1.334e-01, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 4.417e-03, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 3

[########################################] | 100% Completed | 101.78 ms
[########################################] | 100% Completed | 102.30 ms

DefaultConv: Convergence check
  REWS: delta = 1.750e-02, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 3.497e-04, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 4

[########################################] | 100% Completed | 101.09 ms
[########################################] | 100% Completed | 102.45 ms

DefaultConv: Convergence check
  REWS: delta = 1.682e-03, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 2.523e-05, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 5

[########################################] | 100% Completed | 170.14 ms
[########################################] | 100% Completed | 103.89 ms

DefaultConv: Convergence check
  REWS: delta = 1.512e-04, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 2.853e-06, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 6

[########################################] | 100% Completed | 101.71 ms
[########################################] | 100% Completed | 104.75 ms

DefaultConv: Convergence check
  REWS: delta = 1.804e-05, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 2.983e-07, lim = 1.000e-07  --  FAILED

Algorithm Iterative: Iteration 7

[########################################] | 100% Completed | 102.28 ms
[########################################] | 100% Completed | 101.83 ms

DefaultConv: Convergence check
  REWS: delta = 1.859e-06, lim = 1.000e-06  --  FAILED
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 3.149e-08, lim = 1.000e-07  --  OK

Algorithm Iterative: Iteration 8

[########################################] | 100% Completed | 101.62 ms
[########################################] | 100% Completed | 102.30 ms

DefaultConv: Convergence check
  REWS: delta = 1.958e-07, lim = 1.000e-06  --  OK
  TI  : delta = 0.000e+00, lim = 1.000e-07  --  OK
  CT  : delta = 3.232e-09, lim = 1.000e-07  --  OK

Algorithm Iterative: Convergence reached.

[########################################] | 100% Completed | 101.46 ms
[########################################] | 100% Completed | 101.48 ms
<xarray.Dataset> Size: 8kB
Dimensions:     (state: 1, turbine: 48)
Coordinates:
  * state       (state) int64 8B 0
Dimensions without coordinates: turbine
Data variables: (12/20)
    AMB_P       (state, turbine) float64 384B 3.731e+03 3.731e+03 ... 3.731e+03
    AMB_REWS    (state, turbine) float64 384B 8.0 8.0 8.0 8.0 ... 8.0 8.0 8.0
    AMB_RHO     (state, turbine) float64 384B 1.225 1.225 1.225 ... 1.225 1.225
    AMB_TI      (state, turbine) float64 384B 0.04 0.04 0.04 ... 0.04 0.04 0.04
    AMB_WD      (state, turbine) float64 384B 270.0 270.0 270.0 ... 270.0 270.0
    CT          (state, turbine) float64 384B 0.8165 0.8171 ... 0.9116 0.9042
    ...          ...
    YAW         (state, turbine) float64 384B 270.0 270.0 270.0 ... 270.0 270.0
    order       (state, turbine) int64 384B 0 1 2 3 4 5 6 ... 42 43 44 45 46 47
    order_inv   (state, turbine) int64 384B 0 1 2 3 4 5 6 ... 42 43 44 45 46 47
    order_ssel  (state, turbine) int64 384B 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
    weight      (state, turbine) float64 384B 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0
    tname       (turbine) <U3 576B 'T0' 'T1' 'T2' 'T3' ... 'T45' 'T46' 'T47'

The following visualizes the flow field around the most south-westerly turbine in the wind farm:

In [6]:

o = foxes.output.FlowPlots2D(algo, farm_results)
g = o.gen_states_fig_xy(
    FV.WS,
    figsize=(6, 5),
    resolution=2,
    xmin=-300, xmax=500,
    ymin=-200, ymax=200,
    levels=20,
    quiver_n=15,
    quiver_pars={"scale": 0.25, "alpha": 0.3}
)
fig = next(g)
plt.show()

../_images/notebooks_blockage_farm_12_0.png

The induction zone in front of the rotor is clearly visible. For the complete wind farm, the individual induction zones merge:

In [7]:

o = foxes.output.FlowPlots2D(algo, farm_results)
g = o.gen_states_fig_xy(
    FV.WS,
    figsize=(5, 5),
    resolution=10,
    xspace=2000,
    yspace=1000,
    levels=30,
)
fig = next(g)
plt.show()

../_images/notebooks_blockage_farm_14_0.png

The individual contributions to the blockage zone are visible when plotting the wind speed along a straight line in y direction, at fixed x in front of the first row of turbines:

In [8]:

H = algo.mbook.turbine_types["DTU10MW"].H
n_points = 10000
points = np.zeros((1, n_points, 3))
points[:, :, 1] = np.linspace(-1000, 4000, n_points)[None, :]
points[:, :, 2] = H

xlist = [-500, -400, -300, -200]
fig, ax = plt.subplots(figsize=(10, 4))
for x in xlist:
    points[:, :, 0] = x
    point_results = algo.calc_points(farm_results, points)
    ax.plot(point_results[FV.WS][0, :], points[0, :, 1], label=f"x = {x} m")
ax.set_ylabel("y [m]")
ax.set_xlabel("Wind speed [m/s]")
ax.legend()
plt.show()

../_images/notebooks_blockage_farm_16_0.png

Now we investigate the effect of mirroring turbines to account for the speed up effect due to the constrictions from the ground and boundary layer. For this we create a wind farm which consists of only one turbine:

In [9]:

farm = foxes.WindFarm()
farm.add_turbine(foxes.Turbine(xy=[0., 0.], turbine_models=["DTU10MW"]))

Turbine 0, T0: xy=(0.00, 0.00), DTU10MW

Here we calculate farm results for 3 scenarios using only the Rankine Half Body model: 1) No ground effect 2) Using the ground as mirror 3) Using the height 500m as a mirror for the boundary layer

A flow plot for each scenario is generated, along with a windspeed profile comparison.

In [10]:

# set up figure
fig, axs = plt.subplots(1, 3, figsize=(9, 3), sharex=True, sharey=True)
ax = axs.ravel()
vmin = states.ws * 0.5
vmax = None
npoints = 10000
profile_pts = np.ndarray(shape=(npoints, 2, 3))

bmodel = "RHB"
mirror_combo = {
    "No mirror": None,
    "Ground mirror": "ground_mirror",
    "Ground and 500m mirror": "blh_mirror_h500",
}

with foxes.utils.runners.DaskRunner() as runner:
    for i, (mirror, gmodel) in enumerate(mirror_combo.items()):

        algo = foxes.algorithms.Iterative(
            farm,
            states,
            rotor_model="centre",
            wake_models=[bmodel],
            wake_frame="rotor_wd",
            partial_wakes="grid36",
            ground_models=gmodel,
            chunks={FC.STATE: None, FC.POINT: 4000},
            verbosity=0,
        )

        farm_results = runner.run(algo.calc_farm)

        o = foxes.output.FlowPlots2D(algo, farm_results, runner=runner)
        g = o.gen_states_fig_xz(
            var=FV.WS, resolution=2,
            fig=fig, ax=ax[i],
            xmin=-300, xmax=500,
            zmin=0, zmax=500,
            levels=20,
            vmin=vmin, vmax=vmax,
            title=mirror,
        )
        fig = next(g)

        # get point data for profile plot and store in df
        points = np.zeros((1, npoints, 3))
        points[:, :, 0] = -200
        points[:, :, 1] = 0
        points[:, :, 2] = np.linspace(0, 500, npoints)[None, :]
        point_results = runner.run(algo.calc_points, args=(farm_results, points))
        profile_pts[:, 0, i] = point_results[FV.WS][0, :] # x points, WS
        profile_pts[:, 1, i] = points[0, :, 2] #  y points, z values

fig.tight_layout()
plt.show()


Algorithm Iterative: Iteration 0

[########################################] | 100% Completed | 101.32 ms
[########################################] | 100% Completed | 106.25 ms
[########################################] | 100% Completed | 102.10 ms
[########################################] | 100% Completed | 103.93 ms
[########################################] | 100% Completed | 101.48 ms
[########################################] | 100% Completed | 104.09 ms
[########################################] | 100% Completed | 101.41 ms
[########################################] | 100% Completed | 102.15 ms
[########################################] | 100% Completed | 103.62 ms
[########################################] | 100% Completed | 103.03 ms

Algorithm Iterative: Iteration 0

[########################################] | 100% Completed | 101.95 ms
[                                        ] | 0% Completed | 304.93 us

/tmp/ipykernel_2197/3121316729.py:51: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  profile_pts[:, 0, i] = point_results[FV.WS][0, :] # x points, WS

[########################################] | 100% Completed | 102.77 ms
[########################################] | 100% Completed | 101.53 ms
[########################################] | 100% Completed | 101.86 ms
[########################################] | 100% Completed | 101.52 ms
[########################################] | 100% Completed | 105.89 ms
[########################################] | 100% Completed | 101.63 ms
[########################################] | 100% Completed | 102.65 ms
[########################################] | 100% Completed | 103.36 ms
[########################################] | 100% Completed | 102.51 ms

Algorithm Iterative: Iteration 0

[########################################] | 100% Completed | 101.38 ms
[                                        ] | 0% Completed | 158.00 us

/tmp/ipykernel_2197/3121316729.py:51: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  profile_pts[:, 0, i] = point_results[FV.WS][0, :] # x points, WS

[########################################] | 100% Completed | 102.43 ms
[########################################] | 100% Completed | 101.54 ms
[########################################] | 100% Completed | 107.01 ms
[########################################] | 100% Completed | 101.16 ms
[########################################] | 100% Completed | 103.61 ms
[########################################] | 100% Completed | 101.66 ms
[########################################] | 100% Completed | 101.13 ms
[########################################] | 100% Completed | 101.95 ms
[########################################] | 100% Completed | 101.99 ms

/tmp/ipykernel_2197/3121316729.py:51: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  profile_pts[:, 0, i] = point_results[FV.WS][0, :] # x points, WS

../_images/notebooks_blockage_farm_20_6.png

Now we plot the wind speed profile in front of the row of turbines at x=-200 and y=0. Here we observe the additional blockage effect caused by the mirrored turbines:

In [11]:

fig, ax = plt.subplots()
lines = ["-","--","-.",":"]
for i, mirror, in enumerate(mirror_combo):
    ax.plot(profile_pts[:,0,i],profile_pts[:,1,i], label=mirror, linestyle=lines[i])
ax.set_title(f'Profile comparison')
ax.set_ylabel("z")
ax.set_xlabel("Wind speed [m/s]")
plt.legend()
plt.show()

../_images/notebooks_blockage_farm_22_0.png

While the ground mirror seems to have a noticable effect in this example, the reflection from the imagined top of the boundary layer at 500 m seems to be small at the rotor relevant heights.

Let’s finish by looking at the flow in a more realistic modelling setting. As an example, we apply the SelfSimilar2020 induction model with ground effects to a stable ambient wind speed profile, and combine this with the Bastankhah2016 wake deficit model without ground mirroring:

In [12]:

states = foxes.input.states.SingleStateStates(
    ws=8.0, wd=270.0, ti=0.04, rho=1.225, z0=1e-4, MOL=200, H=120,
    profiles={FV.WS: "ABLLogStableWsProfile"})

fig, axs = plt.subplots(2, 1, figsize=(8, 5), sharex=True, sharey=True)
ax = axs.ravel()

bmodels = {"None": [], "SelfSimilar2020": ["SelfSimilar2020"]}
with foxes.utils.runners.DaskRunner() as runner:
    for i, (bname, wakes0) in enumerate(bmodels.items()):

        algo = foxes.algorithms.Iterative(
            farm,
            states,
            rotor_model="grid36",
            wake_models=wakes0 + ["IECTI2019_quartic", "Bastankhah2016_linear_lim_k004"],
            wake_frame="rotor_wd",
            partial_wakes="rotor_points",
            ground_models={"SelfSimilar2020": "ground_mirror"},
            chunks={FC.STATE: None, FC.POINT: 4000},
            verbosity=0,
        )

        farm_results = runner.run(algo.calc_farm)

        o = foxes.output.FlowPlots2D(algo, farm_results, runner=runner)
        g = o.gen_states_fig_xz(
            var=FV.WS, resolution=2,
            xmin=-1000, xmax=3000,
            zmin=0, zmax=1000,
            levels=20,
            fig=fig,
            ax=ax[i],
            title=f"Induction: {bname}"
        )
        fig = next(g)

fig.tight_layout()
plt.show()


Algorithm Iterative: Iteration 0

[########################################] | 100% Completed | 101.30 ms
[########################################] | 100% Completed | 101.86 ms
[########################################] | 100% Completed | 101.36 ms
[########################################] | 100% Completed | 104.59 ms
[########################################] | 100% Completed | 102.76 ms
[########################################] | 100% Completed | 102.04 ms

Algorithm Iterative: Iteration 0

[########################################] | 100% Completed | 102.19 ms
[########################################] | 100% Completed | 102.00 ms
[########################################] | 100% Completed | 102.21 ms
[########################################] | 100% Completed | 101.54 ms
[########################################] | 100% Completed | 101.46 ms
[########################################] | 100% Completed | 101.78 ms

../_images/notebooks_blockage_farm_25_1.png