Eddy viscosity wake model

docs/references.bib

@@ -273,6 +273,19 @@ @Article{bay_2022
 DOI = {10.5194/wes-2022-17}
 }
 
+@article{gunn2019evmodel,
+	title = {Improvements to the Eddy Viscosity Wind Turbine Wake Model},
+	volume = {1222},
+	issn = {1742-6596},
+	url = {https://dx.doi.org/10.1088/1742-6596/1222/1/012003},
+	doi = {10.1088/1742-6596/1222/1/012003},
+	pages = {012003},
+	number = {1},
+	journaltitle = {Journal of Physics: Conference Series},
+	author = {Gunn, K.},
+	note = {Publisher: {IOP} Publishing},
+}
+
 @article{Pedersen_2022_turbopark2,
 url = {https://dx.doi.org/10.1088/1742-6596/2265/2/022063},
 year = {2022},
@@ -299,3 +312,16 @@ @article{SinnerFleming2024grs
 title = {Robust wind farm layout optimization},
 journal = {Journal of Physics: Conference Series},
 }
+
+@proceedings{Kuo2014soed,
+    author = {Kuo, Jim Y. J. and Romero, David A. and Amon, Cristina H.},
+    title = "{A Novel Wake Interaction Model for Wind Farm Layout Optimization}",
+    volume = {Volume 6B: Energy},
+    series = {ASME International Mechanical Engineering Congress and Exposition},
+    pages = {V06BT07A074},
+    year = {2014},
+    month = {11},
+    doi = {10.1115/IMECE2014-39073},
+    url = {https://doi.org/10.1115/IMECE2014-39073},
+    eprint = {https://asmedigitalcollection.asme.org/IMECE/proceedings-pdf/IMECE2014/46521/V06BT07A074/4267418/v06bt07a074-imece2014-39073.pdf},
+}

docs/wake_models.ipynb

@@ -202,6 +202,26 @@
     "model_plot(\"../examples/inputs/turboparkgauss_cubature.yaml\", include_wake_deflection=False)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Eddy Viscosity\n",
+    "\n",
+    "The eddy viscosity model is a wake model originally presented by {cite:t}`ainslie1988calculating` and more recently advanced by {cite:t}`gunn2019evmodel`. The definition of the wake model requires solving an ordinary differential equation (ODE) to determine the wake centerline deficit at various downstream distances. Following {cite:t}`gunn2019evmodel` and others, the FLORIS implementation makes the self-similarity assumption to simplify the wake centerline ODE. The wake profile is Gaussian.\n",
+    "\n",
+    "Further, the wake width for individual turbine wakes is \"expanded\" upon encountering another turbine, and the eddy viscosity model is designed to be run with the \"sum of energy deficit\" (SOED) wake combination model."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_plot(\"../examples/inputs/ev.yaml\", include_wake_deflection=False)"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -360,6 +380,15 @@
     "combination_plot(\"sosfs\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Sum of Energy Deficit (SOED)\n",
+    "\n",
+    "This model combines the wakes by conserving momentum, as described in {cite:t}`Kuo2014soed`."
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",

examples/examples_eddy_viscosity/001_demonstrate_eddy_viscosity_model.py

@@ -0,0 +1,124 @@
+"""Example: Eddy viscosity model
+This example demonstrates the wake profile of the eddy viscosity wake model,
+presented originally by Ainslie (1988) and updated by Gunn (2019).
+Links:
+- Ainslie (1988): https://doi.org/10.1016/0167-6105(88)90037-2
+- Gunn (2019): https://dx.doi.org/10.1088/1742-6596/1222/1/012003
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import floris.flow_visualization as flowviz
+import floris.layout_visualization as layoutviz
+from floris import FlorisModel
+
+
+# Instantiate FLORIS model
+fmodel = FlorisModel("../inputs/ev.yaml")
+
+## Plot the flow velocity profiles
+# Set up a two-turbine farm
+D = 126
+fmodel.set(layout_x=[0, 500, 1000], layout_y=[0, 0, 0])
+
+fig, ax = plt.subplots(1, 1)
+fig.set_size_inches(10, 4)
+ax.scatter(fmodel.layout_x, fmodel.layout_y, color="red", label="Turbine location")
+
+# Set a single wind condition of westerly wind
+wd_array = np.array([270])
+ws_array = 8.0 * np.ones_like(wd_array)
+ti_array = 0.06 * np.ones_like(wd_array)
+fmodel.set(wind_directions=wd_array, wind_speeds=ws_array, turbulence_intensities=ti_array)
+
+# Create points to sample the flow in the turbines' wakes, both at the centerline and
+# across the wake's width
+x_locs = np.linspace(0, 2000, 400)
+y_locs = np.linspace(-100, 100, 41)
+points_x, points_y = np.meshgrid(x_locs, y_locs)
+points_x = points_x.flatten()
+points_y = points_y.flatten()
+points_z = 90 * np.ones_like(points_x)
+
+# Collect the points
+u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z)
+u_at_points = u_at_points.reshape((len(y_locs), len(x_locs), 1))
+
+# Plot the flow velocities
+for y_idx, y in enumerate(y_locs):
+    a = 1-np.abs(y/100)
+    if y_idx == (len(y_locs)-1)/2:
+        label = r"$r = 0D$ (centerline)"
+    elif y_idx == 3*(len(y_locs)-1)/4:
+        label = r"$r \approx 0.4D$"
+    elif y_idx == len(y_locs)-1:
+        label = r"$r \approx 0.8D$"
+    else:
+        label=None
+    ax.plot(x_locs, u_at_points[y_idx, :, 0].flatten(), color="black", alpha=a, label=label)
+ax.grid()
+ax.legend()
+ax.set_xlabel("x location [m]")
+ax.set_ylabel("Wind Speed [m/s]")
+
+## Visualize the flow in aligned and slightly misaligned conditions for a 9-turbine farm
+fig, ax = plt.subplots(2,2)
+fig.set_size_inches(10, 10)
+D = 126.0
+x_locs = np.array([0, 5*D, 10*D])
+y_locs = x_locs
+points_x, points_y = np.meshgrid(x_locs, y_locs)
+fmodel.set(
+    layout_x = points_x.flatten(),
+    layout_y = points_y.flatten()
+)
+
+# Aligned
+layoutviz.plot_turbine_rotors(fmodel, ax=ax[0,0])
+layoutviz.plot_turbine_labels(fmodel, ax=ax[0,0])
+
+fmodel.set(
+    wind_speeds=[8.0, 8.0],
+    wind_directions=[270.0, 270.0+15.0],
+    turbulence_intensities=[0.06, 0.06]
+)
+fmodel.run()
+cut_plane = fmodel.calculate_horizontal_plane(height=90, findex_for_viz=0)
+
+flowviz.visualize_cut_plane(cut_plane, ax=ax[0,0])
+ax[0,0].set_title("Aligned flow")
+
+# Plot and print the power output of the turbines in each row
+np.set_printoptions(formatter={"float": "{0:0.3f}".format})
+print("Aligned case:")
+for i in range(3):
+    idxs = [3*i, 3*i+1, 3*i+2]
+    ax[1,0].scatter([0,1,2], fmodel.get_turbine_powers()[0, idxs]/1e6, label="Column {0}".format(i))
+    print(idxs, " -- ", fmodel.get_turbine_powers()[0, idxs]/1e6)
+ax[1,0].grid()
+ax[1,0].legend()
+ax[1,0].set_xlabel("Turbine in column")
+ax[1,0].set_ylabel("Power [MW]")
+ax[1,0].set_ylim([0.5, 1.8])
+
+# Misaligned
+layoutviz.plot_turbine_rotors(fmodel, yaw_angles=(-15.0)*np.ones(9), ax=ax[0,1])
+layoutviz.plot_turbine_labels(fmodel, ax[0,1])
+
+cut_plane = fmodel.calculate_horizontal_plane(height=90, findex_for_viz=1)
+flowviz.visualize_cut_plane(cut_plane, ax=ax[0,1])
+ax[0,1].set_title("Misaligned flow")
+
+print("\nMisaligned case:")
+for i in range(3):
+    idxs = [3*i, 3*i+1, 3*i+2]
+    ax[1,1].scatter([0,1,2], fmodel.get_turbine_powers()[1, idxs]/1e6, label="Column {0}".format(i))
+    print(idxs, " -- ", fmodel.get_turbine_powers()[1, idxs]/1e6)
+ax[1,1].grid()
+ax[1,1].legend()
+ax[1,1].set_xlabel("Turbine in column")
+ax[1,1].set_ylabel("Power [MW]")
+ax[1,1].set_ylim([0.5, 1.8])
+
+plt.show()

examples/examples_eddy_viscosity/002_reproduce_published_results.py

@@ -0,0 +1,153 @@
+"""Example: Reproduce published eddy viscosity results
+This example attempts to reproduce the results of Ainslie (1988) and Gunn (2019)
+using the FLORIS implementation of the eddy viscosity model.
+
+Links:
+- Ainslie (1988): https://doi.org/10.1016/0167-6105(88)90037-2
+- Gunn (2019): https://dx.doi.org/10.1088/1742-6596/1222/1/012003
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel
+from floris.turbine_library import build_cosine_loss_turbine_dict
+
+
+# Build a constant CT turbine model for use in comparisons (not realistic)
+u_0 = 8.0 # wind speed [m/s]
+
+# Load the EV model
+fmodel = FlorisModel("../inputs/ev.yaml")
+
+## First, attempt to reproduce Figs. 3 and 4 from Ainslie (1988).
+
+# It is not clear exactly how the parametrization of Gunn, which is the one
+# that is implemented in the EddyViscosityVelocity model, can be matched to
+# the parameters given by Ainslie. Rather than try to do this, we simply
+# generate plots similar to Figs. 3 and 4. using the default parameters following
+# Gunn (2019).
+
+# Generate figure to plot on
+fig, ax = plt.subplots(1, 2, sharex=True, sharey=True)
+fig.set_size_inches(8, 4)
+
+HH = 50
+D = 50
+
+for i, wd_std_ev in enumerate([0.0, np.rad2deg(0.1)]):
+    fmodel.set_param(["wake", "wake_velocity_parameters", "eddy_viscosity", "wd_std_ev"], wd_std_ev)
+    for C_T, ls in zip([0.5, 0.7, 0.9], ["-", "--", ":"]):
+        const_CT_turb = build_cosine_loss_turbine_dict(
+            turbine_data_dict={
+                "wind_speed":[0.0, 30.0],
+                "power":[0.0, 1.0], # Not realistic but won't be used here
+                "thrust_coefficient":[C_T, C_T]
+            },
+            turbine_name="ConstantCT",
+            rotor_diameter=D,
+            hub_height=HH,
+            ref_tilt=0.0,
+        )
+
+        # Load the EV model and set a constant CT turbine
+        fmodel.set(
+            layout_x=[0],
+            layout_y=[0],
+            turbine_type=[const_CT_turb],
+            wind_speeds=[u_0],
+            wind_directions=[270],
+            turbulence_intensities=[0.14], # As specified by Ainslie
+            wind_shear=0.0,
+            reference_wind_height=HH
+        )
+
+        points_x = np.linspace(2*D, 10*D, 1000)
+        points_y = np.zeros_like(points_x)
+        points_z = HH * np.ones_like(points_x)
+
+        u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z)
+
+        # Plot results (not different axis scales in Ainslie)
+        ax[i].plot(
+            points_x/D, 1-u_at_points[0, :]/u_0, color="k", linestyle=ls, label=rf"$C_T$ = {C_T}"
+        )
+        ax[i].set_title(r"WD std. dev. $\sigma_{\theta} = $"+"{:.1f} deg".format(wd_std_ev))
+
+ax[0].set_xlabel("Downstream distance [D]")
+ax[1].set_xlabel("Downstream distance [D]")
+ax[0].set_ylabel("Centerline velocity deficit [-]")
+ax[0].set_ylim([0, 1])
+ax[0].legend()
+ax[0].grid()
+ax[1].legend()
+ax[1].grid()
+
+
+## Second, reproduce Figure 1b (right) from Gunn (2019)
+
+# Reset to the default model parameters
+fmodel = FlorisModel("../inputs/ev.yaml")
+
+# Adjustments
+C_T2 = 0.34
+D2 = 1.52
+
+y_offset = -20
+HH = 0.8
+
+# Match depends on ambient turbulence intensity. 7.5% appears close. Also, switch off meandering.
+TI = 0.075
+fmodel.set_param(["wake", "wake_velocity_parameters", "eddy_viscosity", "wd_std_ev"], 0.0)
+
+turb_1 = build_cosine_loss_turbine_dict(
+    turbine_data_dict={
+        "wind_speed":[0.0, 30.0],
+        "power":[0.0, 1.0], # Not realistic but won't be used here
+        "thrust_coefficient":[0.8, 0.8]
+    },
+    turbine_name="ConstantCT1",
+    rotor_diameter=1,
+    hub_height=HH,
+    ref_tilt=0.0,
+)
+turb_2 = build_cosine_loss_turbine_dict(
+    turbine_data_dict={
+        "wind_speed":[0.0, 30.0],
+        "power":[0.0, 1.0], # Not realistic but won't be used here
+        "thrust_coefficient":[C_T2, C_T2]
+    },
+    turbine_name="ConstantCT2",
+    rotor_diameter=D2,
+    hub_height=HH,
+    ref_tilt=0.0,
+)
+
+fmodel.set(
+    layout_x=[0, 1.5],
+    layout_y=[0, y_offset],
+    turbine_type=[turb_1, turb_2],
+    wind_speeds=[u_0],
+    wind_directions=[270.0],
+    turbulence_intensities=[TI],
+    wind_shear=0.0,
+    reference_wind_height=HH
+)
+
+n_pts = 1000
+points_x = np.concatenate((np.linspace(0, 10, n_pts), np.linspace(1.5, 11.5, n_pts)))
+points_y = np.concatenate((np.zeros(n_pts), y_offset*np.ones(n_pts)))
+points_z = HH * np.ones_like(points_x)
+u_at_points = fmodel.sample_flow_at_points(points_x, points_y, points_z)
+
+fig, ax = plt.subplots(1,1)
+ax.plot([0.0, 0.0], [0.3, 1.1], color="black")
+ax.plot([1.5, 1.5], [0.3, 1.1], color="black")
+ax.plot(points_x[:n_pts], u_at_points[0, :n_pts]/u_0, color="C0")
+ax.plot(points_x[n_pts:], u_at_points[0, n_pts:]/u_0, color="C2", linestyle="--")
+ax.grid()
+ax.set_xlabel(r"$X$")
+ax.set_ylabel(r"$\tilde{U}_c$")
+ax.set_ylim([0.3, 1.1])
+
+plt.show()