Add MIT yaw correction model (3rd pass)

docs/operation_models_user.ipynb

@@ -460,18 +460,117 @@
     "ax[1].set_ylabel(\"Power [kW]\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "92912bf7",
+   "metadata": {},
+   "source": [
+    "### Unified Momentum Model\n",
+    "\n",
+    "User-level name: `\"unified-momentum\"`\n",
+    "\n",
+    "Underlying class: `UnifiedMomentumModelTurbine`\n",
+    "\n",
+    "Required data on `power_thrust_table`:\n",
+    "- `ref_air_density` (scalar)\n",
+    "- `ref_tilt` (scalar)\n",
+    "- `wind_speed` (list)\n",
+    "- `power` (list)\n",
+    "- `thrust_coefficient` (list)\n",
+    "\n",
+    "An extension of the classical one-dimensional momentum theory to model the induction of an\n",
+    "actuator disk is presented in {cite:t}`HeckJohlasHowland2023_yawed_adm` to directly account\n",
+    "for power and thrust loss due to yaw misalignment rather than using an empirical correction\n",
+    "as in the cosine loss model. Analytical expressions for the induction, thrust, initial wake\n",
+    "velocities and power are developed as a function of the yaw angle and thrust coefficient.\n",
+    "\n",
+    "Note that the low thrust limit of the Unified Momentum Model is presently implemented in FLORIS, which returns the equations derived and validated in Heck et al. (2023).\n",
+    "This low thrust limit will be accurate for thrust coefficients approximately less than 0.9.\n",
+    "\n",
+    "This section recreates key validation figures discussed in the paper through FLORIS."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "92912bf7",
+   "id": "8bbac518",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_points = 20\n",
+    "\n",
+    "fmodel.set_operation_model(\"unified-momentum\")\n",
+    "fmodel.set(layout_x=[0.0], layout_y=[0.0])\n",
+    "fmodel.reset_operation()\n",
+    "fmodel.set(\n",
+    "    wind_data=TimeSeries(\n",
+    "        wind_speeds=np.array(n_points * [11.0]),\n",
+    "        wind_directions=np.array(n_points * [270.0]),\n",
+    "        turbulence_intensities=0.06\n",
+    "    )\n",
+    ")\n",
+    "yaw_angles = np.linspace(0, 50, n_points)\n",
+    "cos_reference = np.cos(np.radians(yaw_angles))\n",
+    "cos3_reference = np.cos(np.radians(yaw_angles))**3\n",
+    "\n",
+    "fmodel.set(yaw_angles=np.reshape(yaw_angles, (-1,1)))\n",
+    "fmodel.run()\n",
+    "\n",
+    "powers = fmodel.get_turbine_powers()\n",
+    "power_ratio_umm = powers[:,0] / powers[0,0]\n",
+    "\n",
+    "fig, ax = plt.subplots(1,1)\n",
+    "ax.plot(yaw_angles, power_ratio_umm, label=\"Unified momentum model\", color=\"black\")\n",
+    "ax.plot(yaw_angles, cos_reference, label=\"$\\cos(\\gamma)$\", linestyle=\":\", color=\"purple\")\n",
+    "ax.plot(yaw_angles, cos3_reference, label=\"$\\cos^3(\\gamma)$\", linestyle=\":\", color=\"orange\")\n",
+    "ax.grid()\n",
+    "ax.legend()\n",
+    "ax.set_title(\"Figure 2 (a): Power ratio vs yaw angle\")\n",
+    "ax.set_xlabel(\"Yaw angle, $\\gamma$ (degrees)\")\n",
+    "ax.set_ylabel(\"Power ratio, $P(\\gamma)/P(0)$\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "30e54ab2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from floris.core.turbine.unified_momentum_model import Heck, LimitedHeck, MomentumSolution\n",
+    "\n",
+    "n_points = 20\n",
+    "yaw_angles = np.linspace(0, 50, n_points)\n",
+    "\n",
+    "ct_prime = 1.33\n",
+    "\n",
+    "heck = Heck()\n",
+    "ai_umm = np.array([heck(ct_prime, np.radians(yaw)).an for yaw in yaw_angles])\n",
+    "heck_no_spanwise = LimitedHeck()\n",
+    "ai_no_spanwise = np.array([heck_no_spanwise(ct_prime, np.radians(yaw)).an for yaw in yaw_angles])\n",
+    "\n",
+    "fig, ax = plt.subplots(1,1)\n",
+    "ax.plot(yaw_angles, ai_umm / ai_umm[0], label=\"Yaw-dependent UMM\", color=\"black\")\n",
+    "ax.plot(yaw_angles, ai_no_spanwise/ ai_no_spanwise[0], label=\"Low outlet spanwise velocity limit\", linestyle=\"--\", color=\"blue\")\n",
+    "ax.grid()\n",
+    "ax.legend()\n",
+    "ax.set_title(\"Figure 3: Normalized rotor-normal, rotor-averaged induction for the yawed UMM\")\n",
+    "ax.set_xlabel(\"Yaw angle, $\\gamma$ (degrees)\")\n",
+    "ax.set_ylabel(\"Axial induction ratio, $a_n(\\gamma) / a_n(0)$\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "00172dfa",
    "metadata": {},
    "outputs": [],
    "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "floris",
    "language": "python",
    "name": "python3"
   },
@@ -485,7 +584,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.13.2"
   }
  },
  "nbformat": 4,

docs/references.bib

@@ -299,3 +299,15 @@ @article{SinnerFleming2024grs
 title = {Robust wind farm layout optimization},
 journal = {Journal of Physics: Conference Series},
 }
+
+@article{HeckJohlasHowland2023_yawed_adm,
+doi = {10.1017/jfm.2023.129},
+url = {https://doi.org/10.1017/jfm.2023.129},
+year = {2023},
+month = {mar},
+publisher={Cambridge University Press},
+volume = {959},
+author = {K.S. Heck, H.M. Johlas and M.F. Howland},
+title = {Modelling the induction, thrust and power of a yaw-misaligned actuator disk},
+journal = {Journal of Fluid Mechanics},
+}

examples/examples_turbine/004_compare_yaw_loss.py

@@ -0,0 +1,87 @@
+"""
+Example: Change operation model and compare power loss in yaw.
+
+This example illustrates how to define different operational models and compares
+the power loss resulting from yaw misalignment across these various models.
+"""
+
+import itertools
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel, TimeSeries
+
+
+# Parameters
+N = 101  # How many steps to cover yaw range in
+yaw_max = 30  # Maximum yaw angle to test
+
+# Set up the yaw angle sweep
+yaw_angles = np.zeros((N, 1))
+yaw_angles[:, 0] = np.linspace(-yaw_max, yaw_max, N)
+# print(yaw_angles.shape)
+
+
+def evaluate_yawed_power(wsp: float, op_model: str) -> float:
+    print(f"Evaluating model: {op_model}   wind speed: {wsp} m/s")
+
+    # Grab model of FLORIS
+    fmodel = FlorisModel("../inputs/gch.yaml")
+
+    # Run N cases by setting up a TimeSeries (which is just several independent simulations)
+    wind_directions = np.ones(N) * 270.0
+    fmodel.set(
+        wind_data=TimeSeries(
+            wind_speeds=wsp,
+            wind_directions=wind_directions,
+            turbulence_intensities=0.06,
+        )
+    )
+
+    yaw_angles = np.array(
+        [(yaw, 0.0, 0.0) for yaw in np.linspace(-yaw_max, yaw_max, N)]
+    )
+    fmodel.set_operation_model(op_model)
+    fmodel.set(yaw_angles=yaw_angles)
+    fmodel.run()
+
+    # Save the power output results in kW
+    return fmodel.get_turbine_powers()[:, 0] / 1000
+
+
+# Loop over the operational models and wind speeds to compare
+op_models = ["simple", "cosine-loss", "unified-momentum"]
+wind_speeds = [11.0, 11.5, 15.0]
+results = {}
+for op_model, wsp in itertools.product(op_models, wind_speeds):
+
+    # Save the power output results in kW
+    results[(op_model, wsp)] = evaluate_yawed_power(wsp, op_model)
+# Plot the results
+fig, axes = plt.subplots(1, len(wind_speeds), sharey=True)
+
+colors = ["C0", "k", "r"]
+linestyles = ["solid", "dashed", "dotted"]
+for wsp, ax in zip(wind_speeds, axes):
+    ax.set_title(f"wsp: {wsp} m/s")
+    ax.set_xlabel("Yaw angle [deg]")
+    ax.grid(True)
+    for op_model, c, ls in zip(op_models, colors, linestyles):
+
+        upstream_yaw_angle = yaw_angles[:, 0]
+        central_power = results[(op_model, wsp)][upstream_yaw_angle == 0]
+        ax.plot(
+            upstream_yaw_angle,
+            results[(op_model, wsp)] / central_power,
+            label=op_model,
+            color=c,
+            linestyle=ls,
+        )
+
+ax.grid(True)
+ax.legend()
+axes[0].set_xlabel("Yaw angle [deg]")
+axes[0].set_ylabel("Normalized turbine power [-]")
+
+plt.show()

floris/core/turbine/__init__.py

@@ -1,4 +1,5 @@
 
+from floris.core.turbine.unified_momentum_model import UnifiedMomentumModelTurbine
 from floris.core.turbine.operation_models import (
     AWCTurbine,
     CosineLossTurbine,

floris/core/turbine/turbine.py

@@ -15,6 +15,7 @@
 from floris.core.turbine import (
     AWCTurbine,
     CosineLossTurbine,
+    UnifiedMomentumModelTurbine,
     MixedOperationTurbine,
     PeakShavingTurbine,
     SimpleDeratingTurbine,
@@ -41,6 +42,7 @@
         "mixed": MixedOperationTurbine,
         "awc": AWCTurbine,
         "peak-shaving": PeakShavingTurbine,
+        "unified-momentum": UnifiedMomentumModelTurbine,
     },
 }