Torque 2024 - 4 turbine cases#

These cases place 4 turbines in a line with a spacing of 5 rotor diameters, and the turbine diameters are 240 meters. The layout is:

  • x: [0.0, 1200.0, 2400.0, 3600.0]

  • y: [0.0, 0.0, 0.0, 0.0]

The wind speed is 9.8 m/s and the wind direction is 270 degrees.

The following analytical wake models are included in FLORIS, FOXES, and PyWake.

Velocity Models

Wake Model

FLORIS

FOXES

PyWake

Jensen 1983

Larson 2009

Bastankhah / Porte Agel 2014

Bastankhah / Porte Agel 2016

Niayifar / Porté-Agel 2016

IEA Task37 Bastankhah 2018

Carbajo Fuertes / Markfort / Porté-Agel 2018

Blondel / Cathelain 2020

Zong / Porté-Agel 2020

Cumulative Curl 2022

TurbOPark (Nygaard 2022)

•*

Empirical Gauss 2023

Deflection Models

Wake Model

FLORIS

FOXES

PyWake

Jímenez 2010

Bastankhah / Porté-Agel 2016

Larsen et al 2020

Empirical Gauss 2023

# stdlib
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np

# wcomp
from wcomp import WCompFloris, WCompPyWake, WCompFoxes
from wcomp.plotting import plot_plane

# plot settings
PROFILE_LINEWIDTH = 1.0
ERROR_LINEWIDTH = 1.5

SMALL_SIZE = 12
MEDIUM_SIZE = 14
BIGGER_SIZE = 16
# plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('axes', labelsize=MEDIUM_SIZE)
plt.rc('axes', titlesize=BIGGER_SIZE)
plt.rc('legend', fontsize=MEDIUM_SIZE)
plt.rc('figure', titlesize=BIGGER_SIZE)

# Constants for all cases. Copy to a particular code block to change, as needed.
CASE_DIR = Path('cases_torque2024/four_turbine')
this_case = CASE_DIR / Path('jensen/wind_energy_system.yaml')
floris_case = WCompFloris(this_case)
ROTOR_D = floris_case.rotor_diameter
XMIN = -1 * ROTOR_D
XMAX = (20+15) * ROTOR_D
YMIN = -2 * ROTOR_D
YMAX =  2 * ROTOR_D

SAVE_FIGS = False
SAVE_PATH = Path('figures')
if SAVE_FIGS:
    SAVE_PATH.mkdir(exist_ok=True)

Case Layout#

y_turbine = np.array([-ROTOR_D/2, ROTOR_D/2])
x_streamwise = np.array([XMIN, XMAX])
y_streamwise = np.array([0.0, 0.0])
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot([0, 0], y_turbine, color="black", linestyle='-', linewidth=3)
ax.plot([5*ROTOR_D, 5*ROTOR_D], y_turbine, color="black", linestyle='-', linewidth=3)
ax.plot([10*ROTOR_D, 10*ROTOR_D], y_turbine, color="black", linestyle='-', linewidth=3)
ax.plot([15*ROTOR_D, 15*ROTOR_D], y_turbine, color="black", linestyle='-', linewidth=3, label="Turbine locations")
ax.plot(x_streamwise, y_streamwise, '-.', color='black', linewidth=2, label="Streamwise velocity profile")
ax.set_title("Four-turbine velocity profiles")
ax.set_xlabel("X (m)")
ax.set_ylabel("Y (m)")
ax.set_ylim([-1000, 1000])
ax.axis('equal')
ax.grid()
ax.legend()

if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('4t config.png'))
_images/e172f87ca3b3f60388bdf0173b95b02ded578a3c5803cbd1b736f1bb084ddd1a.png

Velocity models#

Jensen#

This is the Jensen model as described in the paper: A note on wind generator interaction (1983) https://backend.orbit.dtu.dk/ws/portalfiles/portal/55857682/ris_m_2411.pdf

There, the wake width is given as 0.1 on pages 6 and 8. PyWake has 0.1 as default FLORIS has 0.05 as default, but it is changed to 0.1 here. FOXES has a set of predefined Jensen wake models available, but you are expected to create your own, as needed. That is done here with k=0.1.

Implemented in: FLORIS, FOXES, PyWake

this_case = CASE_DIR / Path('jensen/wind_energy_system.yaml')
floris_1turbine = WCompFloris(this_case)
foxes_1turbine = WCompFoxes(this_case)
pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
ax.plot([0, 0], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([5*ROTOR_D, 5*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([10*ROTOR_D, 10*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([15*ROTOR_D, 15*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH, label="Turbine locations")
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
x3, y3 = lines[2].get_data()
e1 = np.abs(y1 - y2)
e2 = np.abs(y2 - y3)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FLORIS - FOXES|")
ax.plot(x1, e2, color="black", linestyle=':', linewidth=ERROR_LINEWIDTH, label="|FOXES - PyWake|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('jensen 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))
print(np.linalg.norm(e2, ord=2))


# fig, ax = plt.subplots(3, 1)
# fig.suptitle("Cross section velocity profile")
# X_D = [1, 5, 10]
# for i, D_X in enumerate(X_D):
#     plt.axes(ax[i])
#     floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     ax[i].set_title(f"{D_X} D")
#     ax[i].set_ylabel("U (m/s)")
#     ax[i].set_ybound(lower=0.0, upper=12.0)
#     ax[i].grid()
#     if i < len(X_D) - 1:
#         ax[i].xaxis.set_ticklabels([])
#     else:
#         ax[i].set_xlabel("Y (m)")
#         ax[i].legend()
Turbine 0, T0: windio_turbine
Turbine 1, T1: windio_turbine
Turbine 2, T2: windio_turbine
Turbine 3, T3: windio_turbine
Error norm:
2.6329332064277464
3.5667386439657554
_images/58c13ac627b0c26679f2e910a85920badc2bb251169e5fc4d9bfc2dea40d6b1c.png

Bastankhah / Porte Agel 2014#

A new analytical model for wind-turbine wakes https://www.sciencedirect.com/science/article/abs/pii/S0960148114000317

k*, the wake growth rate d-sigma/dx, is the only parameter

Implemented in: FOXES, PyWake

this_case = CASE_DIR / Path('bastankhah2014/wind_energy_system.yaml')
# floris_1turbine = WCompFloris(this_case)  # Not implemented
foxes_1turbine = WCompFoxes(this_case)
pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
# floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
ax.plot([0, 0], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([5*ROTOR_D, 5*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([10*ROTOR_D, 10*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([15*ROTOR_D, 15*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH, label="Turbine locations")
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
e1 = np.abs(y1 - y2)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FOXES - PyWake|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('b2014 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))


# fig, ax = plt.subplots(3, 1)
# fig.suptitle("Cross section velocity profile")
# X_D = [1, 5, 10]
# for i, D_X in enumerate(X_D):
#     plt.axes(ax[i])
#     # floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     ax[i].set_title(f"{D_X} D")
#     ax[i].set_ylabel("U (m/s)")
#     ax[i].set_ybound(lower=0.0, upper=12.0)
#     ax[i].grid()
#     if i < len(X_D) - 1:
#         ax[i].xaxis.set_ticklabels([])
#     else:
#         ax[i].set_xlabel("Y (m)")
#         ax[i].legend()
Turbine 0, T0: windio_turbine
Turbine 1, T1: windio_turbine
Turbine 2, T2: windio_turbine
Turbine 3, T3: windio_turbine
Error norm:
4.925589047967199
_images/68d136ed8e18aee6fb062546407bfb5934d34f04d2dd2f3b521bd5cdebb2a43f.png

Bastankhah / Porte Agel 2016#

Experimental and theoretical study of wind turbine wakes in yawed conditions https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/experimental-and-theoretical-study-of-wind-turbine-wakes-in-yawed-conditions/400859134E38F340C8811FD7AAD0CE05

ky, kz are wake growth rates in the y and z directions Eq 6.16 has parameters alpha and beta and below says that they are 0.58 and 0.077. Equations 7.1 - 7.4 summarize the model but use parameters alpha_star and beta_star, which are 2.32 and 0.154. These are related to alpha and beta by:

  • alpha_star = 4 * alpha = 4 * 0.58 = 2.32

  • beta_star = 2 * beta = 2 * 0.077 = 0.154

In the near wake, the distribution away from the wake centerline is given by 6.13.

Implemented in: FOXES, PyWake

NOTE: k is set to 0.022 as mentioned in the last paragraph of section 4.

NOTE: FLORIS has ky and kz parameters, but FOXES uses the same, k, for both. This simplification is mentioned in the paper in section 7 Model Predictions. For simplicity, FLORIS uses k for both ky and kz.

NOTE: FLORIS sets ky and kz in the deflection model based on Eq 15 from Niayifar / Porte Agel 2016 https://www.mdpi.com/1996-1073/9/9/741 ky = kz = ka * I + kb ka = 0.3837 kb = 0.003678 Where as Bastankhah / Porte Agel 2016 sets ky and kz directly as in the velocity model.

this_case = CASE_DIR / Path('bastankhah2016/wind_energy_system.yaml')
floris_1turbine = WCompFloris(this_case)
foxes_1turbine = WCompFoxes(this_case)
# pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
# pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
ax.plot([0, 0], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([5*ROTOR_D, 5*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([10*ROTOR_D, 10*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([15*ROTOR_D, 15*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH, label="Turbine locations")
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
e1 = np.abs(y1 - y2)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FLORIS - FOXES|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('b2016 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))


# fig, ax = plt.subplots(3, 1)
# fig.suptitle("Cross section velocity profile")
# X_D = [1, 5, 10]
# for i, D_X in enumerate(X_D):
#     plt.axes(ax[i])
#     floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     # pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     ax[i].set_title(f"{D_X} D")
#     ax[i].set_ylabel("U (m/s)")
#     ax[i].set_ybound(lower=0.0, upper=12.0)
#     ax[i].grid()
#     if i < len(X_D) - 1:
#         ax[i].xaxis.set_ticklabels([])
#     else:
#         ax[i].set_xlabel("Y (m)")
#         ax[i].legend()
Turbine 0, T0: windio_turbine
Turbine 1, T1: windio_turbine
Turbine 2, T2: windio_turbine
Turbine 3, T3: windio_turbine
Error norm:
13.013058028397586
_images/2af1c4b00d042fd84a02dc1f27db4d86fcdbc62357946a583a80fbc0197204fb.png

TurbOPark#

Modelling cluster wakes and wind farm blockage https://iopscience.iop.org/article/10.1088/1742-6596/1618/6/062072

The only parameter is A with value = 0.04

Implemented in: FLORIS, FOXES, PyWake, but FLORIS doesn't support plotting so it's excluded.

this_case = CASE_DIR / Path('turbopark/wind_energy_system.yaml')
# floris_1turbine = WCompFloris(this_case)
foxes_1turbine = WCompFoxes(this_case)
pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
# floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
ax.plot([0, 0], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([5*ROTOR_D, 5*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([10*ROTOR_D, 10*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([15*ROTOR_D, 15*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH, label="Turbine locations")
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
e1 = np.abs(y1 - y2)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FOXES - PyWake|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('turbopark 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))


fig, ax = plt.subplots(3, 1)
fig.suptitle("Cross section velocity profile")
X_D = [1, 5, 10]
for i, D_X in enumerate(X_D):
    plt.axes(ax[i])
    # floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
    foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
    pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
    ax[i].set_title(f"{D_X} D")
    ax[i].set_ylabel("U (m/s)")
    ax[i].set_ybound(lower=0.0, upper=12.0)
    ax[i].grid()
    if i < len(X_D) - 1:
        ax[i].xaxis.set_ticklabels([])
    else:
        ax[i].set_xlabel("Y (m)")
        ax[i].legend()
Turbine 0, T0: windio_turbine
Turbine 1, T1: windio_turbine
Turbine 2, T2: windio_turbine
Turbine 3, T3: windio_turbine
Error norm:
3.3160717557378208
_images/a1467bb33ec43f246b1ac9860ed3c8b359a280af42e015b9a8eeb3f67e38860a.png _images/48a4e99ca0bc953aa50a97c75589e6a8569c1c36661ed07253ae79644ea060e7.png

Wake deflection models#

Here we yaw the turbines 10 degrees to activate the wake deflection models.

Jimenez#

Application of a LES technique to characterize the wake deflection of a wind turbine in yaw https://onlinelibrary.wiley.com/doi/10.1002/we.380

Using the Jensen velocity deficit model, we compare the Jimenez deflection models.

beta is the only parameter, and the paper suggest values between 0.9 and 0.125. We use 0.1.

NOTE: The FLORIS implementation does not correlate to the paper, and there are additional parameters. ad, bd are default to 0 so that was left. I've set kd to beta.

Implemented in: FLORIS and PyWake

this_case = CASE_DIR / Path('jimenez/wind_energy_system.yaml')
floris_1turbine = WCompFloris(this_case)
# foxes_1turbine = WCompFoxes(this_case)
pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
# foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
ax.plot([0, 0], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([5*ROTOR_D, 5*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([10*ROTOR_D, 10*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH)
ax.plot([15*ROTOR_D, 15*ROTOR_D], [0, 10], color="black", linestyle='-', linewidth=PROFILE_LINEWIDTH, label="Turbine locations")
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
e1 = np.abs(y1 - y2)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FLORIS - PyWake|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('jimenez 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))


# fig, ax = plt.subplots(3, 1)
# fig.suptitle("Cross section velocity profile")
# X_D = [1, 5, 10]
# for i, D_X in enumerate(X_D):
#     plt.axes(ax[i])
#     floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     # foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     ax[i].set_title(f"{D_X} D")
#     ax[i].set_ylabel("U (m/s)")
#     ax[i].set_ybound(lower=0.0, upper=12.0)
#     ax[i].grid()
#     if i < len(X_D) - 1:
#         ax[i].xaxis.set_ticklabels([])
#     else:
#         ax[i].set_xlabel("Y (m)")
#         ax[i].legend()


fig, ax = plt.subplots(3, 1)
fig.suptitle("Horizontal streamwise velocity contour")

plt.axes(ax[0])
ax[0].xaxis.set_ticklabels([])
floris_plane = floris_1turbine.horizontal_contour(wind_direction=270)

plt.axes(ax[1])
ax[1].xaxis.set_ticklabels([])
# foxes_plane = foxes_1turbine.horizontal_contour(wind_direction=270)

plt.axes(ax[2])
ax[1].xaxis.set_ticklabels([])
pywake_plane = pywake_1turbine.horizontal_contour(wind_direction=270)


fig, ax = plt.subplots(2, 1)
fig.suptitle("Diff: FLORIS as basis")

plt.axes(ax[0])
diff = floris_plane - pywake_plane
abs_diff = np.abs(diff.values)
min_speed = -1 * np.max(abs_diff)
max_speed = np.max(abs_diff)
plot_plane(
    diff,
    min_speed=min_speed,
    max_speed=max_speed,
    cmap='PuOr',
    clevels=100,
    color_bar=True,
    title="FLORIS - PyWake"
)

# plt.axes(ax[1])
# diff = floris_plane - foxes_plane
# abs_diff = np.abs(diff.values)
# min_speed = -1 * np.max(abs_diff)
# max_speed = np.max(abs_diff)
# plot_plane(
#     diff,
#     ax=ax[1],
#     min_speed=min_speed,
#     max_speed=max_speed,
#     cmap='PuOr',
#     clevels=100,
#     color_bar=True,
#     title="FLORIS - FOXES"
# )
Error norm:
4.274196853600336
_images/208dbcb3eb93081bbf10fa409967fd194c705602965ab53c890b4c5dea65f779.png _images/b5e964f394151c625f5b9095bf6e24130807c55b1f5868ec4db50d0a8a57132f.png _images/83404f9ccbade7035cb9a333644eec3dd1cebfd21fa7fb4000cdd2d6a4d7c222.png

Bastankhah / Porte Agel 2016#

Experimental and theoretical study of wind turbine wakes in yawed conditions https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/experimental-and-theoretical-study-of-wind-turbine-wakes-in-yawed-conditions/400859134E38F340C8811FD7AAD0CE05

ky, kz are wake growth rates in the y and z directions Eq 6.16 has parameters alpha and beta and below says that they are 0.58 and 0.077. Equations 7.1 - 7.4 summarize the model but use parameters alpha_star and beta_star, which are 2.32 and 0.154. These are related to alpha and beta by:

  • alpha_star = 4 * alpha = 4 * 0.58 = 2.32

  • beta_star = 2 * beta = 2 * 0.077 = 0.154

In the near wake, the distribution away from the wake centerline is given by 6.13.

Implemented in: FLORIS, FOXES PyWake implements the velocity model, but not the deflection model.

NOTE: k is set to 0.022 as mentioned in the last paragraph of section 4.

NOTE: FLORIS has ky and kz parameters, but FOXES uses the same, k, for both. This simplification is mentioned in the paper in section 7 Model Predictions. For simplicity, FLORIS uses k for both ky and kz.

NOTE: FLORIS sets ky and kz in the deflection model based on Eq 15 from Niayifar / Porte Agel 2016 https://www.mdpi.com/1996-1073/9/9/741 ky = kz = ka * I + kb ka = 0.3837 kb = 0.003678 Where as Bastankhah / Porte Agel 2016 sets ky and kz directly as in the velocity model.

this_case = CASE_DIR / Path('bastankhah2016_deflection/wind_energy_system.yaml')
floris_1turbine = WCompFloris(this_case)
foxes_1turbine = WCompFoxes(this_case)
# pywake_1turbine = WCompPyWake(this_case)

fig, ax = plt.subplots(figsize=(12.8, 4.8))
floris_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
foxes_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
# pywake_1turbine.streamwise_profile_plot(wind_direction=270, y_coordinate=0.0, xmin=XMIN, xmax=XMAX)
lines = ax.lines
x1, y1 = lines[0].get_data()
x2, y2 = lines[1].get_data()
e1 = np.abs(y1 - y2)
emax = np.max([np.max(e1), np.max(e2)])
ax.plot(x1, e1, color="black", linestyle='-.', linewidth=ERROR_LINEWIDTH, label="|FLORIS - FOXES|")
ax.set_title("Four-turbine streamwise velocity profile")
ax.set_xlabel("X (m)")
ax.set_ylabel('U (m/s)')
ax.set_ybound(lower=0.0)
ax.legend()
ax.grid()
if SAVE_FIGS: fig.savefig(SAVE_PATH / Path('b2016 deflection 4t.png'))

print("Error norm:")
print(np.linalg.norm(e1, ord=2))


# fig, ax = plt.subplots(3, 1)
# fig.suptitle("Cross section velocity profile")
# X_D = [1, 5, 10]
# for i, D_X in enumerate(X_D):
#     plt.axes(ax[i])
#     floris_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     foxes_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     # pywake_1turbine.xsection_profile_plot(wind_direction=270, x_coordinate=D_X * ROTOR_D, ymin=YMIN, ymax=YMAX)
#     ax[i].set_title(f"{D_X} D")
#     ax[i].set_ylabel("U (m/s)")
#     ax[i].set_ybound(lower=0.0, upper=12.0)
#     ax[i].grid()
#     if i < len(X_D) - 1:
#         ax[i].xaxis.set_ticklabels([])
#     else:
#         ax[i].set_xlabel("Y (m)")
#         ax[i].legend()


fig, ax = plt.subplots(3, 1)
fig.suptitle("Horizontal streamwise velocity contour")

plt.axes(ax[0])
# ax[0].xaxis.set_ticklabels([])
# pywake_plane = pywake_1turbine.horizontal_contour(wind_direction=270)

plt.axes(ax[1])
ax[1].xaxis.set_ticklabels([])
floris_plane = floris_1turbine.horizontal_contour(wind_direction=270)

plt.axes(ax[2])
ax[1].xaxis.set_ticklabels([])
foxes_plane = foxes_1turbine.horizontal_contour(wind_direction=270)


fig, ax = plt.subplots(2, 1)
fig.suptitle("Diff: FLORIS as basis")

plt.axes(ax[0])
# diff = floris_plane - pywake_plane
# abs_diff = np.abs(diff.values)
# min_speed = -1 * np.max(abs_diff)
# max_speed = np.max(abs_diff)
# plot_plane(
#     diff,
#     ax=ax[0],
#     min_speed=min_speed,
#     max_speed=max_speed,
#     cmap='PuOr',
#     clevels=100,
#     color_bar=True,
#     title="FLORIS - PyWake"
# )

plt.axes(ax[1])
diff = floris_plane - foxes_plane
abs_diff = np.abs(diff.values)
min_speed = -1 * np.max(abs_diff)
max_speed = np.max(abs_diff)
plot_plane(
    diff,
    min_speed=min_speed,
    max_speed=max_speed,
    cmap='PuOr',
    clevels=100,
    color_bar=True,
    title="FLORIS - FOXES"
)
Turbine 0, T0: windio_turbine
Turbine 1, T1: windio_turbine
Turbine 2, T2: windio_turbine
Turbine 3, T3: windio_turbine
Error norm:
12.375514689019521
_images/efa552e941a63ae5b825f07d66e06a3ff3091c591fa99956118f4a794faf9973.png _images/f62d23a39d2b71220f010a48e705e31b192e2139803be90662f90cbe75f1a546.png _images/93c5067e1b6a4c5e2c1451bd9b6fd1f6f88e948e7a59e4852521b5c503f4946d.png