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]:
mbook = foxes.ModelBook()
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()
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 applies under-relaxation of the variable FV.CT, corresponding to 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. This avoids 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(
mbook,
farm,
states=states,
rotor_model="grid36",
wake_models=["RHB", "Bastankhah025_linear_lim_k004"],
wake_frame="rotor_wd",
partial_wakes_model="auto",
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 'grid36'
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'
Partial wakes 'auto': Applying RotorPoints to ['RHB']
Partial wakes 'auto': Applying PartialAxiwake to ['Bastankhah025_linear_lim_k004']
Initializing model 'Madsen'
Initializing model 'RHB'
Initializing model 'Bastankhah025_linear_lim_k004'
Initializing model 'RotorPoints1'
Initializing model 'PartialAxiwake10'
Initializing model 'auto'
--------------------------------------------------
Running Iterative: calc_farm
--------------------------------------------------
n_states : 1
n_turbines: 48
--------------------------------------------------
states : SingleStateStates (SingleStateStates)
rotor : grid36 (GridRotor)(n=6)
controller: basic_ctrl (BasicFarmController)
partialwks: auto (Mapped)
wake frame: rotor_wd (RotorWD)
--------------------------------------------------
wakes:
0) RHB (RankineHalfBody)
1) Bastankhah025_linear_lim_k004 (Bastankhah2014)(k=0.04, sp=ws_linear_lim)
--------------------------------------------------
turbine models:
0) DTU10MW (PCtFile)
--------------------------------------------------
Initializing model 'SetXYHD'
Initializing model 'SetXYHD_t2f'
Initializing model 'calc_yaw_CentreRotor1'
Initializing model 'CalcOrder'
Initializing model 'SetAmbFarmResults'
Initializing model 'URelax_CT'
Initializing model 'FarmWakesCalculation'
Initializing model 'Iterative_calc'
--------------------------------------------------
Model oder
--------------------------------------------------
00) SetXYHD_t2f
01) basic_ctrl
02) calc_yaw_CentreRotor1
03) grid36
04) CalcOrder
05) basic_ctrl
05.0) Post-rotor: DTU10MW
06) SetAmbFarmResults
07) FarmWakesCalculation
--------------------------------------------------
[########################################] | 100% Completed | 101.41 ms
Input data:
<xarray.Dataset>
Dimensions: (state: 1, turbine: 48, tmodels: 1)
Coordinates:
* state (state) int64 0
* tmodels (tmodels) <U7 'DTU10MW'
Dimensions without coordinates: turbine
Data variables:
weight (state, turbine) float64 dask.array<chunksize=(1, 48), meta=np.ndarray>
tmodel_sels (state, turbine, tmodels) bool dask.array<chunksize=(1, 48, 1), meta=np.ndarray>
Output 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, weight
Chunks: {'state': None, 'point': 4000}
Calculating 1 states for 48 turbines
[########################################] | 100% Completed | 101.70 ms
Algorithm Iterative: Iteration 1
[########################################] | 100% Completed | 101.99 ms
[########################################] | 100% Completed | 102.89 ms
DefaultConv: Convergence check
REWS: delta = 1.414e+00, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 8.024e-03, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 2
[########################################] | 100% Completed | 102.43 ms
[########################################] | 100% Completed | 101.81 ms
DefaultConv: Convergence check
REWS: delta = 6.799e-01, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 1.190e-03, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 3
[########################################] | 100% Completed | 101.99 ms
[########################################] | 100% Completed | 102.00 ms
DefaultConv: Convergence check
REWS: delta = 3.019e-01, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 1.659e-03, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 4
[########################################] | 100% Completed | 101.67 ms
[########################################] | 100% Completed | 161.72 ms
DefaultConv: Convergence check
REWS: delta = 1.254e-01, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 3.777e-04, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 5
[########################################] | 100% Completed | 101.46 ms
[########################################] | 100% Completed | 102.54 ms
DefaultConv: Convergence check
REWS: delta = 4.966e-02, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 3.488e-04, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 6
[########################################] | 100% Completed | 101.51 ms
[########################################] | 100% Completed | 101.65 ms
DefaultConv: Convergence check
REWS: delta = 1.876e-02, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 1.125e-04, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 7
[########################################] | 100% Completed | 101.99 ms
[########################################] | 100% Completed | 104.89 ms
DefaultConv: Convergence check
REWS: delta = 6.906e-03, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 7.697e-05, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 8
[########################################] | 100% Completed | 102.29 ms
[########################################] | 100% Completed | 102.36 ms
DefaultConv: Convergence check
REWS: delta = 2.448e-03, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 3.103e-05, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 9
[########################################] | 100% Completed | 101.63 ms
[########################################] | 100% Completed | 105.81 ms
DefaultConv: Convergence check
REWS: delta = 8.622e-04, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 1.790e-05, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 10
[########################################] | 100% Completed | 101.80 ms
[########################################] | 100% Completed | 104.29 ms
DefaultConv: Convergence check
REWS: delta = 2.922e-04, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 8.086e-06, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 11
[########################################] | 100% Completed | 101.48 ms
[########################################] | 100% Completed | 103.15 ms
DefaultConv: Convergence check
REWS: delta = 1.008e-04, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 4.324e-06, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 12
[########################################] | 100% Completed | 102.14 ms
[########################################] | 100% Completed | 101.80 ms
DefaultConv: Convergence check
REWS: delta = 3.288e-05, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 2.059e-06, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 13
[########################################] | 100% Completed | 102.10 ms
[########################################] | 100% Completed | 103.71 ms
DefaultConv: Convergence check
REWS: delta = 1.141e-05, lim = 1.000e-05 -- FAILED
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 1.067e-06, lim = 1.000e-06 -- FAILED
Algorithm Iterative: Iteration 14
[########################################] | 100% Completed | 101.55 ms
[########################################] | 100% Completed | 102.49 ms
DefaultConv: Convergence check
REWS: delta = 3.541e-06, lim = 1.000e-05 -- OK
TI : delta = 0.000e+00, lim = 1.000e-06 -- OK
CT : delta = 5.226e-07, lim = 1.000e-06 -- OK
Algorithm Iterative: Convergence reached.
<xarray.Dataset>
Dimensions: (state: 1, turbine: 48)
Coordinates:
* state (state) int64 0
Dimensions without coordinates: turbine
Data variables: (12/25)
AMB_CT (state, turbine) float64 0.814 0.814 0.814 ... 0.814 0.814 0.814
AMB_P (state, turbine) float64 3.731e+03 3.731e+03 ... 3.731e+03
AMB_REWS (state, turbine) float64 8.0 8.0 8.0 8.0 8.0 ... 8.0 8.0 8.0 8.0
AMB_REWS2 (state, turbine) float64 8.0 8.0 8.0 8.0 8.0 ... 8.0 8.0 8.0 8.0
AMB_REWS3 (state, turbine) float64 8.0 8.0 8.0 8.0 8.0 ... 8.0 8.0 8.0 8.0
AMB_RHO (state, turbine) float64 1.225 1.225 1.225 ... 1.225 1.225 1.225
... ...
X (state, turbine) float64 0.0 0.0 0.0 0.0 ... 4e+03 4e+03 4e+03
Y (state, turbine) float64 0.0 400.0 800.0 ... 2.65e+03 3.05e+03
YAW (state, turbine) float64 270.0 270.0 270.0 ... 270.0 270.0 270.0
order (state, turbine) int64 0 1 2 3 4 5 6 7 ... 40 41 45 46 43 44 47
weight (state, turbine) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
tname (turbine) <U3 'T0' 'T1' 'T2' 'T3' ... 'T44' '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()
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()
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 = 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()
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: 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": {},
"Ground mirror": {bmodel: [0]},
"Ground and 500m mirror": {bmodel: [0, 500]},
}
with foxes.utils.runners.DaskRunner() as runner:
for i, (mirror, mirrors) in enumerate(mirror_combo.items()):
algo = foxes.algorithms.Iterative(
mbook,
farm,
states=states,
rotor_model="grid36",
wake_models=[bmodel],
wake_frame="rotor_wd",
partial_wakes_model="auto",
wake_mirrors=mirrors,
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.73 ms
[########################################] | 100% Completed | 102.68 ms
[########################################] | 100% Completed | 101.82 ms
[########################################] | 100% Completed | 108.98 ms
Partial wakes 'auto': Applying RotorPoints to ['RHB']
[########################################] | 100% Completed | 101.37 ms
[########################################] | 100% Completed | 101.86 ms
[########################################] | 100% Completed | 103.74 ms
Algorithm Iterative: Iteration 0
[########################################] | 100% Completed | 102.48 ms
[########################################] | 100% Completed | 103.56 ms
[########################################] | 100% Completed | 102.30 ms
[########################################] | 100% Completed | 101.75 ms
Partial wakes 'auto': Applying RotorPoints to ['WakeMirror_RHB']
[########################################] | 100% Completed | 102.07 ms
[########################################] | 100% Completed | 102.64 ms
[########################################] | 100% Completed | 102.22 ms
Algorithm Iterative: Iteration 0
[########################################] | 100% Completed | 101.70 ms
[########################################] | 100% Completed | 105.87 ms
[########################################] | 100% Completed | 102.12 ms
[########################################] | 100% Completed | 104.99 ms
Partial wakes 'auto': Applying RotorPoints to ['WakeMirror1_RHB']
[########################################] | 100% Completed | 101.98 ms
[########################################] | 100% Completed | 102.06 ms
[########################################] | 100% Completed | 104.11 ms
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()
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:
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()
mbook = foxes.ModelBook()
bmodels = {"None": [], "SelfSimilar2020": ["SelfSimilar2020"]}
with foxes.utils.runners.DaskRunner() as runner:
for i, (bname, wakes0) in enumerate(bmodels.items()):
algo = foxes.algorithms.Iterative(
mbook,
farm,
states=states,
rotor_model="grid36",
wake_models=wakes0 + ["IECTI2019_quartic", "Bastankhah2016_linear_lim_k004"],
wake_frame="rotor_wd",
partial_wakes_model="auto",
wake_mirrors={bmodel: [0]},
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 | 102.27 ms
[########################################] | 100% Completed | 102.95 ms
[########################################] | 100% Completed | 101.99 ms
[########################################] | 100% Completed | 105.27 ms
Algorithm Iterative: Iteration 0
[########################################] | 100% Completed | 102.26 ms
[########################################] | 100% Completed | 101.44 ms
[########################################] | 100% Completed | 102.14 ms
[########################################] | 100% Completed | 103.29 ms
