diff --git a/examples/41_compare_yaw_loss.py b/examples/41_compare_yaw_loss.py new file mode 100644 index 000000000..c30b6887d --- /dev/null +++ b/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() diff --git a/floris/simulation/rotor_velocity.py b/floris/simulation/rotor_velocity.py index 25f94d55d..b57f11a34 100644 --- a/floris/simulation/rotor_velocity.py +++ b/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, diff --git a/floris/simulation/turbine/__init__.py b/floris/simulation/turbine/__init__.py index 355f5c2df..f38615129 100644 --- a/floris/simulation/turbine/__init__.py +++ b/floris/simulation/turbine/__init__.py @@ -14,6 +14,7 @@ from floris.simulation.turbine.operation_models import ( CosineLossTurbine, + MITLossTurbine, MixedOperationTurbine, SimpleDeratingTurbine, SimpleTurbine, diff --git a/floris/simulation/turbine/operation_models.py b/floris/simulation/turbine/operation_models.py index 82c11ee70..7244bf79a 100644 --- a/floris/simulation/turbine/operation_models.py +++ b/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): """ diff --git a/floris/simulation/turbine/turbine.py b/floris/simulation/turbine/turbine.py index f9435facb..09291d1ba 100644 --- a/floris/simulation/turbine/turbine.py +++ b/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, }, diff --git a/tests/turbine_operation_models_test.py b/tests/turbine_operation_models_test.py index 446695855..dcbf59757 100644 --- a/tests/turbine_operation_models_test.py +++ b/tests/turbine_operation_models_test.py @@ -3,6 +3,7 @@ from floris.simulation.turbine.operation_models import ( CosineLossTurbine, + MITLossTurbine, MixedOperationTurbine, POWER_SETPOINT_DEFAULT, rotor_velocity_air_density_correction, @@ -49,6 +50,10 @@ def test_submodel_attributes(): assert hasattr(MixedOperationTurbine, "thrust_coefficient") assert hasattr(MixedOperationTurbine, "axial_induction") + assert hasattr(MITLossTurbine, "power") + assert hasattr(MITLossTurbine, "thrust_coefficient") + assert hasattr(MITLossTurbine, "axial_induction") + def test_SimpleTurbine(): n_turbines = 1 @@ -498,3 +503,109 @@ def test_MixedOperationTurbine(): tilt_angles=tilt_angles_nom, tilt_interp=None ) + +def test_MITLossTurbine(): + + # NOTE: These tests should be updated to reflect actual expected behavior + # of the MITLossTurbine model. Currently, match the CosineLossTurbine model. + + n_turbines = 1 + wind_speed = 10.0 + turbine_data = SampleInputs().turbine + + yaw_angles_nom = 0 * np.ones((1, n_turbines)) + tilt_angles_nom = turbine_data["power_thrust_table"]["ref_tilt"] * np.ones((1, n_turbines)) + yaw_angles_test = 20 * np.ones((1, n_turbines)) + tilt_angles_test = 0 * np.ones((1, n_turbines)) + + + # Check that power works as expected + test_power = MITLossTurbine.power( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=turbine_data["power_thrust_table"]["ref_air_density"], # Matches ref_air_density + yaw_angles=yaw_angles_nom, + tilt_angles=tilt_angles_nom, + tilt_interp=None + ) + truth_index = turbine_data["power_thrust_table"]["wind_speed"].index(wind_speed) + baseline_power = turbine_data["power_thrust_table"]["power"][truth_index] * 1000 + assert np.allclose(baseline_power, test_power) + + # Check that yaw and tilt angle have an effect + test_power = MITLossTurbine.power( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=turbine_data["power_thrust_table"]["ref_air_density"], # Matches ref_air_density + yaw_angles=yaw_angles_test, + tilt_angles=tilt_angles_test, + tilt_interp=None + ) + assert test_power < baseline_power + + # Check that a lower air density decreases power appropriately + test_power = MITLossTurbine.power( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=1.1, + yaw_angles=yaw_angles_nom, + tilt_angles=tilt_angles_nom, + tilt_interp=None + ) + assert test_power < baseline_power + + + # Check that thrust coefficient works as expected + test_Ct = MITLossTurbine.thrust_coefficient( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=1.1, # Unused + yaw_angles=yaw_angles_nom, + tilt_angles=tilt_angles_nom, + tilt_interp=None + ) + baseline_Ct = turbine_data["power_thrust_table"]["thrust_coefficient"][truth_index] + assert np.allclose(baseline_Ct, test_Ct) + + # Check that yaw and tilt angle have the expected effect + test_Ct = MITLossTurbine.thrust_coefficient( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=1.1, # Unused + yaw_angles=yaw_angles_test, + tilt_angles=tilt_angles_test, + tilt_interp=None + ) + absolute_tilt = tilt_angles_test - turbine_data["power_thrust_table"]["ref_tilt"] + assert test_Ct == baseline_Ct * cosd(yaw_angles_test) * cosd(absolute_tilt) + + + # Check that thrust coefficient works as expected + test_ai = MITLossTurbine.axial_induction( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=1.1, # Unused + yaw_angles=yaw_angles_nom, + tilt_angles=tilt_angles_nom, + tilt_interp=None + ) + baseline_misalignment_loss = ( + cosd(yaw_angles_nom) + * cosd(tilt_angles_nom - turbine_data["power_thrust_table"]["ref_tilt"]) + ) + baseline_ai = ( + 1 - np.sqrt(1 - turbine_data["power_thrust_table"]["thrust_coefficient"][truth_index]) + ) / 2 / baseline_misalignment_loss + assert np.allclose(baseline_ai, test_ai) + + # Check that yaw and tilt angle have the expected effect + test_ai = MITLossTurbine.axial_induction( + power_thrust_table=turbine_data["power_thrust_table"], + velocities=wind_speed * np.ones((1, n_turbines, 3, 3)), # 1 findex, 1 turbine, 3x3 grid + air_density=1.1, # Unused + yaw_angles=yaw_angles_test, + tilt_angles=tilt_angles_test, + tilt_interp=None + ) + absolute_tilt = tilt_angles_test - turbine_data["power_thrust_table"]["ref_tilt"] + assert test_Ct == baseline_Ct * cosd(yaw_angles_test) * cosd(absolute_tilt)