Develop MIT Yaw Loss model in V4

examples/41_compare_yaw_loss.py

@@ -0,0 +1,87 @@
+# Copyright 2024 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+import matplotlib.pyplot as plt
+import numpy as np
+import yaml
+
+from floris.tools import FlorisInterface
+
+
+"""
+Example to test out derating of turbines and mixed derating and yawing. Will be refined before
+release. TODO: Demonstrate shutting off turbines also, once developed.
+"""
+
+# Parameters
+N  = 101 # How many steps to cover yaw range in
+yaw_max = 30 # Maximum yaw 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)
+
+
+
+# Now loop over the operational models to compare
+op_models = ["cosine-loss", "mit-loss"]
+results = {}
+
+for op_model in op_models:
+
+    print(f"Evaluating model: {op_model}")
+
+    # Grab model of FLORIS
+    fi = FlorisInterface("inputs/gch.yaml")
+
+    # Initialize to a simple 1 turbine case with n_findex = N
+    fi.reinitialize(layout_x=[0],
+                    layout_y=[0],
+                    wind_directions=270 * np.ones(N),
+                    wind_speeds=8 * np.ones(N),
+                    )
+
+    with open(str(
+        fi.floris.as_dict()["farm"]["turbine_library_path"] /
+        (fi.floris.as_dict()["farm"]["turbine_type"][0] + ".yaml")
+    )) as t:
+        turbine_type = yaml.safe_load(t)
+    turbine_type["power_thrust_model"] = op_model
+
+    # Change the turbine type
+    fi.reinitialize(turbine_type=[turbine_type])
+
+    # Calculate the power
+    fi.calculate_wake(yaw_angles=yaw_angles)
+    turbine_power = fi.get_turbine_powers().squeeze()
+
+    # Save the results
+    results[op_model] = turbine_power
+
+# Plot the results
+fig, ax = plt.subplots()
+
+colors = ["C0", "k", "r"]
+linestyles = ["solid", "dashed", "dotted"]
+for key, c, ls in zip(results, colors, linestyles):
+    central_power = results[key][yaw_angles.squeeze() == 0]
+    ax.plot(yaw_angles.squeeze(), results[key]/central_power, label=key, color=c, linestyle=ls)
+
+ax.grid(True)
+ax.legend()
+ax.set_xlabel("Yaw angle [deg]")
+ax.set_ylabel("Normalized turbine power [deg]")
+
+plt.show()

floris/simulation/rotor_velocity.py

@@ -42,6 +42,18 @@ def rotor_velocity_yaw_correction(
 
     return rotor_effective_velocities
 
+def mit_rotor_velocity_yaw_correction(
+    pP: float,
+    yaw_angles: NDArrayFloat,
+    rotor_effective_velocities: NDArrayFloat,
+) -> NDArrayFloat:
+    # Compute the rotor effective velocity adjusting for yaw settings
+    pW = pP / 3.0  # Convert from pP to w
+    # TODO: cosine loss hard coded
+    rotor_effective_velocities = rotor_effective_velocities * cosd(yaw_angles) ** pW
+
+    return rotor_effective_velocities
+
 def rotor_velocity_tilt_correction(
     tilt_angles: NDArrayFloat,
     ref_tilt: NDArrayFloat,

floris/simulation/turbine/__init__.py

@@ -14,6 +14,7 @@
 
 from floris.simulation.turbine.operation_models import (
     CosineLossTurbine,
+    MITLossTurbine,
     MixedOperationTurbine,
     SimpleDeratingTurbine,
     SimpleTurbine,

floris/simulation/turbine/operation_models.py

@@ -30,6 +30,7 @@
 from floris.simulation.rotor_velocity import (
     average_velocity,
     compute_tilt_angles_for_floating_turbines,
+    mit_rotor_velocity_yaw_correction,
     rotor_velocity_tilt_correction,
     rotor_velocity_yaw_correction,
 )
@@ -318,6 +319,151 @@ def axial_induction(
         misalignment_loss = cosd(yaw_angles) * cosd(tilt_angles - power_thrust_table["ref_tilt"])
         return 0.5 / misalignment_loss * (1 - np.sqrt(1 - thrust_coefficient * misalignment_loss))
 
+
+
+@define
+class MITLossTurbine(BaseOperationModel):
+    """
+    Static class defining an actuator disk turbine model that may be misaligned with the flow.
+    Nonzero tilt and yaw angles are handled via cosine relationships, with the power lost to yawing
+    defined by the pP exponent. This turbine submodel is the default, and matches the turbine
+    model in FLORIS v3.
+
+    As with all turbine submodules, implements only static power() and thrust_coefficient() methods,
+    which are called by power() and thrust_coefficient() on turbine.py, respectively. This class is
+    not intended to be instantiated; it simply defines a library of static methods.
+
+    TODO: Should the turbine submodels each implement axial_induction()?
+    """
+
+    def power(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        air_density: float,
+        yaw_angles: NDArrayFloat,
+        tilt_angles: NDArrayFloat,
+        tilt_interp: NDArrayObject,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        correct_cp_ct_for_tilt: bool = False,
+        **_ # <- Allows other models to accept other keyword arguments
+    ):
+        # Construct power interpolant
+        power_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["power"],
+            fill_value=0.0,
+            bounds_error=False,
+        )
+
+        # Compute the power-effective wind speed across the rotor
+        rotor_average_velocities = average_velocity(
+            velocities=velocities,
+            method=average_method,
+            cubature_weights=cubature_weights,
+        )
+
+        rotor_effective_velocities = rotor_velocity_air_density_correction(
+            velocities=rotor_average_velocities,
+            air_density=air_density,
+            ref_air_density=power_thrust_table["ref_air_density"]
+        )
+
+        rotor_effective_velocities = mit_rotor_velocity_yaw_correction(
+            pP=power_thrust_table["pP"],
+            yaw_angles=yaw_angles,
+            rotor_effective_velocities=rotor_effective_velocities,
+        )
+
+        rotor_effective_velocities = rotor_velocity_tilt_correction(
+            tilt_angles=tilt_angles,
+            ref_tilt=power_thrust_table["ref_tilt"],
+            pT=power_thrust_table["pT"],
+            tilt_interp=tilt_interp,
+            correct_cp_ct_for_tilt=correct_cp_ct_for_tilt,
+            rotor_effective_velocities=rotor_effective_velocities,
+        )
+
+        # Compute power
+        power = power_interpolator(rotor_effective_velocities) * 1e3 # Convert to W
+
+        return power
+
+    def thrust_coefficient(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        yaw_angles: NDArrayFloat,
+        tilt_angles: NDArrayFloat,
+        tilt_interp: NDArrayObject,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        correct_cp_ct_for_tilt: bool = False,
+        **_ # <- Allows other models to accept other keyword arguments
+    ):
+        # Construct thrust coefficient interpolant
+        thrust_coefficient_interpolator = interp1d(
+            power_thrust_table["wind_speed"],
+            power_thrust_table["thrust_coefficient"],
+            fill_value=0.0001,
+            bounds_error=False,
+        )
+
+        # Compute the effective wind speed across the rotor
+        rotor_average_velocities = average_velocity(
+            velocities=velocities,
+            method=average_method,
+            cubature_weights=cubature_weights,
+        )
+
+        # TODO: Do we need an air density correction here?
+        thrust_coefficient = thrust_coefficient_interpolator(rotor_average_velocities)
+        thrust_coefficient = np.clip(thrust_coefficient, 0.0001, 0.9999)
+
+        # Apply tilt and yaw corrections
+        # Compute the tilt, if using floating turbines
+        old_tilt_angles = copy.deepcopy(tilt_angles)
+        tilt_angles = compute_tilt_angles_for_floating_turbines(
+            tilt_angles=tilt_angles,
+            tilt_interp=tilt_interp,
+            rotor_effective_velocities=rotor_average_velocities,
+        )
+        # Only update tilt angle if requested (if the tilt isn't accounted for in the Ct curve)
+        tilt_angles = np.where(correct_cp_ct_for_tilt, tilt_angles, old_tilt_angles)
+
+        thrust_coefficient = (
+            thrust_coefficient
+            * cosd(yaw_angles)
+            * cosd(tilt_angles - power_thrust_table["ref_tilt"])
+        )
+
+        return thrust_coefficient
+
+    def axial_induction(
+        power_thrust_table: dict,
+        velocities: NDArrayFloat,
+        yaw_angles: NDArrayFloat,
+        tilt_angles: NDArrayFloat,
+        tilt_interp: NDArrayObject,
+        average_method: str = "cubic-mean",
+        cubature_weights: NDArrayFloat | None = None,
+        correct_cp_ct_for_tilt: bool = False,
+        **_ # <- Allows other models to accept other keyword arguments
+    ):
+
+        thrust_coefficient = CosineLossTurbine.thrust_coefficient(
+            power_thrust_table=power_thrust_table,
+            velocities=velocities,
+            yaw_angles=yaw_angles,
+            tilt_angles=tilt_angles,
+            tilt_interp=tilt_interp,
+            average_method=average_method,
+            cubature_weights=cubature_weights,
+            correct_cp_ct_for_tilt=correct_cp_ct_for_tilt
+        )
+
+        misalignment_loss = cosd(yaw_angles) * cosd(tilt_angles - power_thrust_table["ref_tilt"])
+        return 0.5 / misalignment_loss * (1 - np.sqrt(1 - thrust_coefficient * misalignment_loss))
+
 @define
 class SimpleDeratingTurbine(BaseOperationModel):
     """

floris/simulation/turbine/turbine.py

@@ -27,6 +27,7 @@
 from floris.simulation import BaseClass
 from floris.simulation.turbine import (
     CosineLossTurbine,
+    MITLossTurbine,
     MixedOperationTurbine,
     SimpleDeratingTurbine,
     SimpleTurbine,
@@ -47,6 +48,7 @@
     "power_thrust_model": {
         "simple": SimpleTurbine,
         "cosine-loss": CosineLossTurbine,
+        "mit-loss": MITLossTurbine,
         "simple-derating": SimpleDeratingTurbine,
         "mixed": MixedOperationTurbine,
     },