From da61128340ac0e520b4583eb759bf771512332f6 Mon Sep 17 00:00:00 2001 From: misi9170 <39596329+misi9170@users.noreply.github.com> Date: Wed, 13 Dec 2023 15:06:05 -0500 Subject: [PATCH] Add utility to build turbine yaml from absolute power, thrust curves (#729) * Structure of file built out. * Main functionality built out. * Test built. * passes tests. * isort, ruff. * Issues with 3.8, 3.9 tests. * isort issue when removing import. * newline? * Reverting; __future__ import to handle testing on 3.8, 3.9 * Now defaults to providing dictionary back; option to print to yaml. * Test updated; ruff, isort. * Improving documentation and exposition of function arguments.' * Update comments, syntax, and typos * Ensure path independence in unit tests --------- Co-authored-by: Rafael M Mudafort --- docs/examples.md | 6 + examples/32_specify_turbine_power_curve.py | 78 +++++++++ floris/turbine_library/__init__.py | 1 + floris/turbine_library/turbine_previewer.py | 2 + floris/turbine_library/turbine_utilities.py | 166 ++++++++++++++++++++ tests/turbine_unit_test.py | 90 +++++++++++ 6 files changed, 343 insertions(+) create mode 100644 examples/32_specify_turbine_power_curve.py create mode 100644 floris/turbine_library/turbine_utilities.py diff --git a/docs/examples.md b/docs/examples.md index e87627cd4..cdd0ae7b0 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -172,6 +172,12 @@ This example builds on example 30. Specifically, fictional data for varying Cp/C wave period, Ts, and wave height, Hs, is used to show the difference in power performance for different wave heights. +### 32_specify_turbine_power_curve.py + +This example demonstrates how to generate a turbine dictionary or yaml input file based on +a specified power and thrust curve. The power and thrust curves may be specified as power +and thrust coefficients or as absolute values. + ## Optimization These examples demonstrate use of the optimization routines diff --git a/examples/32_specify_turbine_power_curve.py b/examples/32_specify_turbine_power_curve.py new file mode 100644 index 000000000..03fbf9978 --- /dev/null +++ b/examples/32_specify_turbine_power_curve.py @@ -0,0 +1,78 @@ +# Copyright 2021 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 floris.tools.visualization as wakeviz +from floris.tools import FlorisInterface +from floris.turbine_library.turbine_utilities import build_turbine_dict + + +""" +This example demonstrates how to specify a turbine model based on a power +and thrust curve for the wind turbine, as well as possible physical parameters +(which default to the parameters of the NREL 5MW reference turbine). + +Note that it is also possible to have a .yaml created, if the file_path +argument to build_turbine_dict is set. +""" + +# Generate an example turbine power and thrust curve for use in the FLORIS model +turbine_data_dict = { + "wind_speed":[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20], + "power_absolute":[0, 30, 200, 500, 1000, 2000, 4000, 4000, 4000, 4000, 4000], + "thrust_coefficient":[0, 0.9, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.25, 0.2] +} + +turbine_dict = build_turbine_dict( + turbine_data_dict, + "example_turbine", + file_path=None, + generator_efficiency=1, + hub_height=90, + pP=1.88, + pT=1.88, + rotor_diameter=126, + TSR=8, + air_density=1.225, + ref_tilt_cp_ct=5 +) + +fi = FlorisInterface("inputs/gch.yaml") +wind_speeds = np.linspace(1, 15, 100) +# Replace the turbine(s) in the FLORIS model with the created one +fi.reinitialize( + layout_x=[0], + layout_y=[0], + wind_speeds=wind_speeds, + turbine_type=[turbine_dict] +) +fi.calculate_wake() + +powers = fi.get_farm_power() + +fig, ax = plt.subplots(1,1) + +ax.scatter(wind_speeds, powers[0,:]/1000, color="C0", s=5, label="Test points") +ax.scatter(turbine_data_dict["wind_speed"], turbine_data_dict["power_absolute"], + color="red", s=20, label="Specified points") + +ax.grid() +ax.set_xlabel("Wind speed [m/s]") +ax.set_ylabel("Power [kW]") +ax.legend() + +plt.show() diff --git a/floris/turbine_library/__init__.py b/floris/turbine_library/__init__.py index 828c50eb2..933615b0c 100644 --- a/floris/turbine_library/__init__.py +++ b/floris/turbine_library/__init__.py @@ -1 +1,2 @@ from floris.turbine_library.turbine_previewer import TurbineInterface, TurbineLibrary +from floris.turbine_library.turbine_utilities import build_turbine_dict diff --git a/floris/turbine_library/turbine_previewer.py b/floris/turbine_library/turbine_previewer.py index c7d1edf72..a582355b6 100644 --- a/floris/turbine_library/turbine_previewer.py +++ b/floris/turbine_library/turbine_previewer.py @@ -12,6 +12,8 @@ # See https://floris.readthedocs.io for documentation +from __future__ import annotations + from pathlib import Path import attrs diff --git a/floris/turbine_library/turbine_utilities.py b/floris/turbine_library/turbine_utilities.py new file mode 100644 index 000000000..c862c21bd --- /dev/null +++ b/floris/turbine_library/turbine_utilities.py @@ -0,0 +1,166 @@ +import os.path + +import numpy as np +import yaml + + +def build_turbine_dict( + turbine_data_dict, + turbine_name, + file_path=None, + generator_efficiency=1.0, + hub_height=90.0, + pP=1.88, + pT=1.88, + rotor_diameter=126.0, + TSR=8.0, + air_density=1.225, + ref_tilt_cp_ct=5.0 +): + """ + Tool for formatting a full turbine dict from data formatted as a + dictionary. + + Default value for turbine physical parameters are from the NREL 5MW reference + wind turbine. + + Returns a turbine dictionary object as expected by FLORIS. Optionally, + prints the dictionary to a yaml to be included in a FLORIS wake model yaml. + + turbine_data is a dictionary that contains keys specifying the + turbine power and thrust as a function of wind speed. The following keys + are possible: + - wind_speed [m/s] + - power_absolute [kW] + - power_coefficient [-] + - thrust_absolute [kN] + - thrust_coefficient [-] + Of these, wind_speed is required. One of power_absolute and power_coefficient + must be specified; and one of thrust_absolute and thrust_coefficient must be + specified. If both _absolute and _coefficient versions are specified, the + _coefficient entry will be used and the _absolute entry ignored. + + Args: + turbine_data_dict (dict): Dictionary containing performance of the wind + turbine as a function of wind speed. Described in more detail above. + turbine_name (string): Name of the turbine, which will be used for the + turbine_type field as well as the filename. + file_path (): Path for placement of the produced yaml. Defaults to None, + in which case no yaml is written. + generator_efficiency (float): Generator efficiency [-]. Defaults to 1.0. + hub_height (float): Hub height [m]. Defaults to 90.0. + pP (float): Cosine exponent for power loss to yaw [-]. Defaults to 1.88. + pT (float): Cosine exponent for thrust loss to yaw [-]. Defaults to 1.88. + rotor_diameter (float). Rotor diameter [m]. Defaults to 126.0. + TSR (float). Turbine optimal tip-speed ratio [-]. Defaults to 8.0. + air_density (float). Air density used to specify power and thrust + curves [kg/m^3]. Defaults to 1.225. + ref_tilt_cp_ct (float). Rotor tilt (due to shaft tilt and/or platform + tilt) used when defining the power and thrust curves [deg]. Defaults + to 5.0. + + Returns: + turbine_dict (dict): Formatted turbine dictionary as expected by FLORIS. + """ + + # Check that necessary columns are specified + if "wind_speed" not in turbine_data_dict: + raise KeyError("wind_speed column must be specified.") + u = np.array(turbine_data_dict["wind_speed"]) + A = np.pi * rotor_diameter**2/4 + + # Construct the Cp curve + if "power_coefficient" in turbine_data_dict: + if "power_absolute" in turbine_data_dict: + print( + "Found both power_absolute and power_coefficient." + "Ignoring power_absolute." + ) + Cp = np.array(turbine_data_dict["power_coefficient"]) + + elif "power_absolute" in turbine_data_dict: + P = np.array(turbine_data_dict["power_absolute"]) + if _find_nearest_value_for_wind_speed(P, u, 10) > 20000 or \ + _find_nearest_value_for_wind_speed(P, u, 10) < 1000: + print( + "Unusual power value detected. Please check that power_absolute", + "is specified in kW." + ) + + validity_mask = (P != 0) | (u != 0) + Cp = np.zeros_like(P, dtype=float) + + Cp[validity_mask] = (P[validity_mask]*1000) / \ + (0.5*air_density*A*u[validity_mask]**3) + + else: + raise KeyError( + "Either power_absolute or power_coefficient must be specified." + ) + + # Construct Ct curve + if "thrust_coefficient" in turbine_data_dict: + if "thrust_absolute" in turbine_data_dict: + print( + "Found both thrust_absolute and thrust_coefficient." + "Ignoring thrust_absolute." + ) + Ct = np.array(turbine_data_dict["thrust_coefficient"]) + + elif "thrust_absolute" in turbine_data_dict: + T = np.array(turbine_data_dict["thrust_absolute"]) + if _find_nearest_value_for_wind_speed(T, u, 10) > 3000 or \ + _find_nearest_value_for_wind_speed(T, u, 10) < 100: + print( + "Unusual thrust value detected. Please check that thrust_absolute", + "is specified in kN." + ) + + validity_mask = (T != 0) | (u != 0) + Ct = np.zeros_like(T) + + Ct[validity_mask] = (T[validity_mask]*1000)/\ + (0.5*air_density*A*u[validity_mask]**2) + + else: + raise KeyError( + "Either thrust_absolute or thrust_coefficient must be specified." + ) + + # Build the turbine dict + power_thrust_dict = { + "wind_speed": u.tolist(), + "power": Cp.tolist(), + "thrust": Ct.tolist() + } + + turbine_dict = { + "turbine_type": turbine_name, + "generator_efficiency": generator_efficiency, + "hub_height": hub_height, + "pP": pP, + "pT": pT, + "rotor_diameter": rotor_diameter, + "TSR": TSR, + "ref_density_cp_ct": air_density, + "ref_tilt_cp_ct": ref_tilt_cp_ct, + "power_thrust_table": power_thrust_dict + } + + # Create yaml file + if file_path is not None: + full_name = os.path.join(file_path, turbine_name+".yaml") + yaml.dump( + turbine_dict, + open(full_name, "w"), + sort_keys=False + ) + + print(full_name, "created.") + + return turbine_dict + +def _find_nearest_value_for_wind_speed(test_vals, ws_vals, ws): + errs = np.absolute(ws_vals-ws) + idx = errs.argmin() + return test_vals[idx] diff --git a/tests/turbine_unit_test.py b/tests/turbine_unit_test.py index 9704483b0..65d433b34 100644 --- a/tests/turbine_unit_test.py +++ b/tests/turbine_unit_test.py @@ -13,9 +13,13 @@ # See https://floris.readthedocs.io for documentation +import os +from pathlib import Path + import attr import numpy as np import pytest +import yaml from scipy.interpolate import interp1d from floris.simulation import ( @@ -31,6 +35,7 @@ compute_tilt_angles_for_floating_turbines, PowerThrustTable, ) +from floris.turbine_library import build_turbine_dict from tests.conftest import SampleInputs, WIND_SPEEDS @@ -563,3 +568,88 @@ def test_asdict(sample_inputs_fixture: SampleInputs): dict2 = new_turb.as_dict() assert dict1 == dict2 + +def test_build_turbine_dict(): + + orig_file_path = Path(__file__).resolve().parent / "data" / "nrel_5MW_custom.yaml" + test_turb_name = "test_turbine_export" + test_file_path = "." + + in_dict = yaml.safe_load( open(orig_file_path, "r") ) + + # Mocked up turbine data + turbine_data_dict = { + "wind_speed":in_dict["power_thrust_table"]["wind_speed"], + "power_coefficient":in_dict["power_thrust_table"]["power"], + "thrust_coefficient":in_dict["power_thrust_table"]["thrust"] + } + + build_turbine_dict( + turbine_data_dict, + test_turb_name, + file_path=test_file_path, + generator_efficiency=in_dict["generator_efficiency"], + hub_height=in_dict["hub_height"], + pP=in_dict["pP"], + pT=in_dict["pT"], + rotor_diameter=in_dict["rotor_diameter"], + TSR=in_dict["TSR"], + air_density=in_dict["ref_density_cp_ct"], + ref_tilt_cp_ct=in_dict["ref_tilt_cp_ct"] + ) + + test_dict = yaml.safe_load( + open(os.path.join(test_file_path, test_turb_name+".yaml"), "r") + ) + + # Correct intended difference for test; assert equal + test_dict["turbine_type"] = in_dict["turbine_type"] + assert list(in_dict.keys()) == list(test_dict.keys()) + assert in_dict == test_dict + + # Now, in absolute values + Cp = np.array(in_dict["power_thrust_table"]["power"]) + Ct = np.array(in_dict["power_thrust_table"]["thrust"]) + ws = np.array(in_dict["power_thrust_table"]["wind_speed"]) + + P = 0.5 * in_dict["ref_density_cp_ct"] * (np.pi * in_dict["rotor_diameter"]**2/4) \ + * Cp * ws**3 + T = 0.5 * in_dict["ref_density_cp_ct"] * (np.pi * in_dict["rotor_diameter"]**2/4) \ + * Ct * ws**2 + + turbine_data_dict = { + "wind_speed":in_dict["power_thrust_table"]["wind_speed"], + "power_absolute": P/1000, + "thrust_absolute": T/1000 + } + + build_turbine_dict( + turbine_data_dict, + test_turb_name, + file_path=test_file_path, + generator_efficiency=in_dict["generator_efficiency"], + hub_height=in_dict["hub_height"], + pP=in_dict["pP"], + pT=in_dict["pT"], + rotor_diameter=in_dict["rotor_diameter"], + TSR=in_dict["TSR"], + air_density=in_dict["ref_density_cp_ct"], + ref_tilt_cp_ct=in_dict["ref_tilt_cp_ct"] + ) + + test_dict = yaml.safe_load( + open(os.path.join(test_file_path, test_turb_name+".yaml"), "r") + ) + + test_dict["turbine_type"] = in_dict["turbine_type"] + assert list(in_dict.keys()) == list(test_dict.keys()) + for k in in_dict.keys(): + if type(in_dict[k]) is dict: + for k2 in in_dict[k].keys(): + assert np.allclose(in_dict[k][k2], test_dict[k][k2]) + elif type(in_dict[k]) is str: + assert in_dict[k] == test_dict[k] + else: + assert np.allclose(in_dict[k], test_dict[k]) + + os.remove( os.path.join(test_file_path, test_turb_name+".yaml") )