From 32112f83de569f75d379b6dfcfe83a86f2cc789f Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 19 Jun 2023 18:19:28 +0200 Subject: [PATCH 01/46] First commit --- floris/turbine_library/t1_vawt_200kW.yaml | 201 ++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 floris/turbine_library/t1_vawt_200kW.yaml diff --git a/floris/turbine_library/t1_vawt_200kW.yaml b/floris/turbine_library/t1_vawt_200kW.yaml new file mode 100644 index 000000000..9cf737908 --- /dev/null +++ b/floris/turbine_library/t1_vawt_200kW.yaml @@ -0,0 +1,201 @@ + +### +# An ID for this type of turbine definition. +# This is not currently used, but it will be enabled in the future. This should typically +# match the root name of the file. +turbine_type: 't1_vawt_200kW' + +### +# Setting for generator losses to power. +generator_efficiency: 1.0 + +### +# Hub height. +hub_height: 40.0 + +### +# Cosine exponent for power loss due to yaw misalignment. +#pP: 1.88 + +### +# Cosine exponent for power loss due to tilt. +#pT: 1.88 + +### +# A vertical-axis wind turbine is characterized by the +# rotor diameter and blade length (in the vertical direction) +rotor_diameter: 26.0 +blade_length: 24.0 + +### +# Tip speed ratio defined as linear blade tip speed normalized by the incoming wind speed. +TSR: 3.8 + +### +# The air density at which the Cp and Ct curves are defined. +#ref_density_cp_ct: 1.225 + +### +# The tilt angle at which the Cp and Ct curves are defined. This is used to capture +# the effects of a floating platform on a turbine's power and wake. +#ref_tilt_cp_ct: 5.0 + +### +# Cp and Ct as a function of wind speed for the turbine's full range of operating conditions. +#power_thrust_table: +# power: +# - 0.0 +# - 0.000000 +# - 0.000000 +# - 0.178085 +# - 0.289075 +# - 0.349022 +# - 0.384728 +# - 0.406059 +# - 0.420228 +# - 0.428823 +# - 0.433873 +# - 0.436223 +# - 0.436845 +# - 0.436575 +# - 0.436511 +# - 0.436561 +# - 0.436517 +# - 0.435903 +# - 0.434673 +# - 0.433230 +# - 0.430466 +# - 0.378869 +# - 0.335199 +# - 0.297991 +# - 0.266092 +# - 0.238588 +# - 0.214748 +# - 0.193981 +# - 0.175808 +# - 0.159835 +# - 0.145741 +# - 0.133256 +# - 0.122157 +# - 0.112257 +# - 0.103399 +# - 0.095449 +# - 0.088294 +# - 0.081836 +# - 0.075993 +# - 0.070692 +# - 0.065875 +# - 0.061484 +# - 0.057476 +# - 0.053809 +# - 0.050447 +# - 0.047358 +# - 0.044518 +# - 0.041900 +# - 0.039483 +# - 0.0 +# - 0.0 +# thrust: +# - 0.0 +# - 0.0 +# - 0.0 +# - 0.99 +# - 0.99 +# - 0.97373036 +# - 0.92826162 +# - 0.89210543 +# - 0.86100905 +# - 0.835423 +# - 0.81237673 +# - 0.79225789 +# - 0.77584769 +# - 0.7629228 +# - 0.76156073 +# - 0.76261984 +# - 0.76169723 +# - 0.75232027 +# - 0.74026851 +# - 0.72987175 +# - 0.70701647 +# - 0.54054532 +# - 0.45509459 +# - 0.39343381 +# - 0.34250785 +# - 0.30487242 +# - 0.27164979 +# - 0.24361964 +# - 0.21973831 +# - 0.19918151 +# - 0.18131868 +# - 0.16537679 +# - 0.15103727 +# - 0.13998636 +# - 0.1289037 +# - 0.11970413 +# - 0.11087113 +# - 0.10339901 +# - 0.09617888 +# - 0.09009926 +# - 0.08395078 +# - 0.0791188 +# - 0.07448356 +# - 0.07050731 +# - 0.06684119 +# - 0.06345518 +# - 0.06032267 +# - 0.05741999 +# - 0.05472609 +# - 0.0 +# - 0.0 +# wind_speed: +# - 0.0 +# - 2.0 +# - 2.5 +# - 3.0 +# - 3.5 +# - 4.0 +# - 4.5 +# - 5.0 +# - 5.5 +# - 6.0 +# - 6.5 +# - 7.0 +# - 7.5 +# - 8.0 +# - 8.5 +# - 9.0 +# - 9.5 +# - 10.0 +# - 10.5 +# - 11.0 +# - 11.5 +# - 12.0 +# - 12.5 +# - 13.0 +# - 13.5 +# - 14.0 +# - 14.5 +# - 15.0 +# - 15.5 +# - 16.0 +# - 16.5 +# - 17.0 +# - 17.5 +# - 18.0 +# - 18.5 +# - 19.0 +# - 19.5 +# - 20.0 +# - 20.5 +# - 21.0 +# - 21.5 +# - 22.0 +# - 22.5 +# - 23.0 +# - 23.5 +# - 24.0 +# - 24.5 +# - 25.0 +# - 25.01 +# - 25.02 +# - 50.0 From a10cfb0005b4c53e7529507ae71ac570a3cad549 Mon Sep 17 00:00:00 2001 From: Gustav_V Date: Wed, 21 Jun 2023 16:53:07 +0200 Subject: [PATCH 02/46] First test to input blade_length for a VAWT --- floris/simulation/farm.py | 16 ++ floris/simulation/floris.py | 1 + floris/simulation/turbine.py | 2 + floris/turbine_library/t1_vawt_200kW.yaml | 201 ---------------------- 4 files changed, 19 insertions(+), 201 deletions(-) delete mode 100644 floris/turbine_library/t1_vawt_200kW.yaml diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index 9aba71fa6..fe2eff6fe 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -72,6 +72,7 @@ class Farm(BaseClass): turbine_fCts: tuple = field(init=False, default=[]) turbine_type_map_sorted: NDArrayObject = field(init=False, default=[]) rotor_diameters_sorted: NDArrayFloat = field(init=False, default=[]) + blade_lengths_sorted: NDArrayFloat = field(init=False, default=[]) TSRs_sorted: NDArrayFloat = field(init=False, default=[]) pPs: NDArrayFloat = field(init=False, default=[]) pPs_sorted: NDArrayFloat = field(init=False, default=[]) @@ -209,6 +210,11 @@ def construct_rotor_diameters(self): turb['rotor_diameter'] for turb in self.turbine_definitions ]) + def construct_blade_lengths(self): + self.blade_lengths = np.array([ + turb['blade_length'] for turb in self.turbine_definitions + ]) + def construct_turbine_TSRs(self): self.TSRs = np.array([turb['TSR'] for turb in self.turbine_definitions]) @@ -271,6 +277,11 @@ def expand_farm_properties( sorted_coord_indices, axis=2 ) + self.blade_lengths_sorted = np.take_along_axis( + self.blade_lengths * template_shape, + sorted_coord_indices, + axis=2 + ) self.TSRs_sorted = np.take_along_axis( self.TSRs * template_shape, sorted_coord_indices, @@ -361,6 +372,11 @@ def finalize(self, unsorted_indices): unsorted_indices[:,:,:,0,0], axis=2 ) + self.blade_lengths = np.take_along_axis( + self.blade_lengths_sorted, + unsorted_indices[:,:,:,0,0], + axis=2 + ) self.TSRs = np.take_along_axis( self.TSRs_sorted, unsorted_indices[:,:,:,0,0], diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 8ca8854e7..8c44c1a81 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -76,6 +76,7 @@ def __attrs_post_init__(self) -> None: self.farm.construct_turbine_power_interps() self.farm.construct_hub_heights() self.farm.construct_rotor_diameters() + self.farm.construct_blade_lengths() self.farm.construct_turbine_TSRs() self.farm.construct_turbine_pPs() self.farm.construct_turbine_pTs() diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index cdf457878..3c16039de 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -587,6 +587,7 @@ class Turbine(BaseClass): Parameters: rotor_diameter (:py:obj: float): The rotor diameter (m). + blade_length (:py:obj: float): The blade length (m). hub_height (:py:obj: float): The hub height (m). pP (:py:obj: float): The cosine exponent relating the yaw misalignment angle to power. @@ -617,6 +618,7 @@ class Turbine(BaseClass): turbine_type: str = field() rotor_diameter: float = field() + blade_length: float = field() hub_height: float = field() pP: float = field() pT: float = field() diff --git a/floris/turbine_library/t1_vawt_200kW.yaml b/floris/turbine_library/t1_vawt_200kW.yaml deleted file mode 100644 index 9cf737908..000000000 --- a/floris/turbine_library/t1_vawt_200kW.yaml +++ /dev/null @@ -1,201 +0,0 @@ - -### -# An ID for this type of turbine definition. -# This is not currently used, but it will be enabled in the future. This should typically -# match the root name of the file. -turbine_type: 't1_vawt_200kW' - -### -# Setting for generator losses to power. -generator_efficiency: 1.0 - -### -# Hub height. -hub_height: 40.0 - -### -# Cosine exponent for power loss due to yaw misalignment. -#pP: 1.88 - -### -# Cosine exponent for power loss due to tilt. -#pT: 1.88 - -### -# A vertical-axis wind turbine is characterized by the -# rotor diameter and blade length (in the vertical direction) -rotor_diameter: 26.0 -blade_length: 24.0 - -### -# Tip speed ratio defined as linear blade tip speed normalized by the incoming wind speed. -TSR: 3.8 - -### -# The air density at which the Cp and Ct curves are defined. -#ref_density_cp_ct: 1.225 - -### -# The tilt angle at which the Cp and Ct curves are defined. This is used to capture -# the effects of a floating platform on a turbine's power and wake. -#ref_tilt_cp_ct: 5.0 - -### -# Cp and Ct as a function of wind speed for the turbine's full range of operating conditions. -#power_thrust_table: -# power: -# - 0.0 -# - 0.000000 -# - 0.000000 -# - 0.178085 -# - 0.289075 -# - 0.349022 -# - 0.384728 -# - 0.406059 -# - 0.420228 -# - 0.428823 -# - 0.433873 -# - 0.436223 -# - 0.436845 -# - 0.436575 -# - 0.436511 -# - 0.436561 -# - 0.436517 -# - 0.435903 -# - 0.434673 -# - 0.433230 -# - 0.430466 -# - 0.378869 -# - 0.335199 -# - 0.297991 -# - 0.266092 -# - 0.238588 -# - 0.214748 -# - 0.193981 -# - 0.175808 -# - 0.159835 -# - 0.145741 -# - 0.133256 -# - 0.122157 -# - 0.112257 -# - 0.103399 -# - 0.095449 -# - 0.088294 -# - 0.081836 -# - 0.075993 -# - 0.070692 -# - 0.065875 -# - 0.061484 -# - 0.057476 -# - 0.053809 -# - 0.050447 -# - 0.047358 -# - 0.044518 -# - 0.041900 -# - 0.039483 -# - 0.0 -# - 0.0 -# thrust: -# - 0.0 -# - 0.0 -# - 0.0 -# - 0.99 -# - 0.99 -# - 0.97373036 -# - 0.92826162 -# - 0.89210543 -# - 0.86100905 -# - 0.835423 -# - 0.81237673 -# - 0.79225789 -# - 0.77584769 -# - 0.7629228 -# - 0.76156073 -# - 0.76261984 -# - 0.76169723 -# - 0.75232027 -# - 0.74026851 -# - 0.72987175 -# - 0.70701647 -# - 0.54054532 -# - 0.45509459 -# - 0.39343381 -# - 0.34250785 -# - 0.30487242 -# - 0.27164979 -# - 0.24361964 -# - 0.21973831 -# - 0.19918151 -# - 0.18131868 -# - 0.16537679 -# - 0.15103727 -# - 0.13998636 -# - 0.1289037 -# - 0.11970413 -# - 0.11087113 -# - 0.10339901 -# - 0.09617888 -# - 0.09009926 -# - 0.08395078 -# - 0.0791188 -# - 0.07448356 -# - 0.07050731 -# - 0.06684119 -# - 0.06345518 -# - 0.06032267 -# - 0.05741999 -# - 0.05472609 -# - 0.0 -# - 0.0 -# wind_speed: -# - 0.0 -# - 2.0 -# - 2.5 -# - 3.0 -# - 3.5 -# - 4.0 -# - 4.5 -# - 5.0 -# - 5.5 -# - 6.0 -# - 6.5 -# - 7.0 -# - 7.5 -# - 8.0 -# - 8.5 -# - 9.0 -# - 9.5 -# - 10.0 -# - 10.5 -# - 11.0 -# - 11.5 -# - 12.0 -# - 12.5 -# - 13.0 -# - 13.5 -# - 14.0 -# - 14.5 -# - 15.0 -# - 15.5 -# - 16.0 -# - 16.5 -# - 17.0 -# - 17.5 -# - 18.0 -# - 18.5 -# - 19.0 -# - 19.5 -# - 20.0 -# - 20.5 -# - 21.0 -# - 21.5 -# - 22.0 -# - 22.5 -# - 23.0 -# - 23.5 -# - 24.0 -# - 24.5 -# - 25.0 -# - 25.01 -# - 25.02 -# - 50.0 From 2a28abc4cf3282d7082a7ba3e4c13e82e064d41a Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 22 Jun 2023 13:07:48 +0200 Subject: [PATCH 03/46] First attempt of inputs to super_gaussian_vawt velocity model --- floris/simulation/solver.py | 1 + floris/simulation/wake.py | 4 +- floris/simulation/wake_velocity/__init__.py | 1 + .../wake_velocity/super_gaussian_vawt.py | 151 ++++++++++++++++++ 4 files changed, 156 insertions(+), 1 deletion(-) create mode 100644 floris/simulation/wake_velocity/super_gaussian_vawt.py diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 5f1767699..bc7cf9bef 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -132,6 +132,7 @@ def sequential_solver( yaw_angle_i = farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = farm.hub_heights_sorted[: ,:, i:i+1, None, None] rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] + #blade_length_i = farm.blade_lengths_sorted[: ,:, i:i+1, None, None] TSR_i = farm.TSRs_sorted[: ,:, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) diff --git a/floris/simulation/wake.py b/floris/simulation/wake.py index a141c94be..8ba32f12a 100644 --- a/floris/simulation/wake.py +++ b/floris/simulation/wake.py @@ -38,6 +38,7 @@ GaussVelocityDeficit, JensenVelocityDeficit, NoneVelocityDeficit, + SuperGaussianVAWTVelocityDeficit, TurbOParkVelocityDeficit, ) @@ -65,7 +66,8 @@ "gauss": GaussVelocityDeficit, "jensen": JensenVelocityDeficit, "turbopark": TurbOParkVelocityDeficit, - "empirical_gauss": EmpiricalGaussVelocityDeficit + "empirical_gauss": EmpiricalGaussVelocityDeficit, + "super_gaussian_vawt": SuperGaussianVAWTVelocityDeficit }, } diff --git a/floris/simulation/wake_velocity/__init__.py b/floris/simulation/wake_velocity/__init__.py index f551f5be8..596b38be3 100644 --- a/floris/simulation/wake_velocity/__init__.py +++ b/floris/simulation/wake_velocity/__init__.py @@ -18,4 +18,5 @@ from floris.simulation.wake_velocity.gauss import GaussVelocityDeficit from floris.simulation.wake_velocity.jensen import JensenVelocityDeficit from floris.simulation.wake_velocity.none import NoneVelocityDeficit +from floris.simulation.wake_velocity.super_gaussian_vawt import SuperGaussianVAWTVelocityDeficit from floris.simulation.wake_velocity.turbopark import TurbOParkVelocityDeficit diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py new file mode 100644 index 000000000..fe3442097 --- /dev/null +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -0,0 +1,151 @@ +# 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. + +from typing import Any, Dict + +import numexpr as ne +import numpy as np +from attrs import define, field + +from floris.simulation import ( + BaseModel, + Farm, + FlowField, + Grid, + Turbine, +) + + +@define +class SuperGaussianVAWTVelocityDeficit(BaseModel): + """ + The Empirical Gauss velocity model has a Gaussian profile + (see :cite:`bastankhah2016experimental` and + :cite:`King2019Controls`) throughout and expands in a (smoothed) + piecewise linear fashion. + + parameter_dictionary (dict): Model-specific parameters. + Default values are used when a parameter is not included + in `parameter_dictionary`. Possible key-value pairs include: + + - **wake_expansion_rates** (*list*): List of expansion + rates for the Gaussian wake width. Must be of length 1 + or greater. + - **breakpoints_D** (*list*): List of downstream + locations, specified in terms of rotor diameters, where + the expansion rates go into effect. Must be one element + shorter than wake_expansion_rates. May be empty. + - **sigma_0_D** (*float*): Initial width of the Gaussian + wake at the turbine location, specified as a multiplier + of the rotor diameter. + - **smoothing_length_D** (*float*): Distance over which + the corners in the piece-wise linear wake expansion rate + are smoothed (specified as a multiplier of the rotor + diameter). + - **mixing_gain_deflection** (*float*): Gain to set the + increase in wake expansion due to wake-induced mixing. + + References: + .. bibliography:: /references.bib + :style: unsrt + :filter: docname in docnames + """ + wake_expansion_coeff_y: float = field(default=0.50) + wake_expansion_coeff_z: float = field(default=0.50) + ay: float = field(default=0.95) + az: float = field(default=4.5) + by: float = field(default=0.35) + bz: float = field(default=0.70) + cy: float = field(default=2.4) + cz: float = field(default=2.4) + + def prepare_function( + self, + grid: Grid, + flow_field: FlowField, + ) -> Dict[str, Any]: + + kwargs = { + "x": grid.x_sorted, + "y": grid.y_sorted, + "z": grid.z_sorted, + "wind_veer": flow_field.wind_veer + } + return kwargs + + def function( + self, + x_i: np.ndarray, + y_i: np.ndarray, + z_i: np.ndarray, + axial_induction_i: np.ndarray, + deflection_field_i: np.ndarray, + yaw_angle_i: np.ndarray, + tilt_angle_i: np.ndarray, + turbulence_intensity_i: np.ndarray, + ct_i: np.ndarray, + hub_height_i: float, + rotor_diameter_i: np.ndarray, + #blade_length_i: np.ndarray, + # enforces the use of the below as keyword arguments and adherence to the + # unpacking of the results from prepare_function() + *, + x: np.ndarray, + y: np.ndarray, + z: np.ndarray, + wind_veer: float + ) -> None: + """ + Calculates the velocity deficits in the wake. + + Args: + x_i (np.array): Streamwise direction grid coordinates of + the ith turbine (m). + y_i (np.array): Cross stream direction grid coordinates of + the ith turbine (m). + z_i (np.array): Vertical direction grid coordinates of + the ith turbine (m) [not used]. + axial_induction_i (np.array): Axial induction factor of the + ith turbine (-) [not used]. + yaw_angle_i (np.array): Yaw angle of the ith turbine (deg). + tilt_angle_i (np.array): Tilt angle of the ith turbine + (deg). + ct_i (np.array): Thrust coefficient for the ith turbine (-). + hub_height_i (float): Hub height for the ith turbine (m). + rotor_diameter_i (np.array): Rotor diameter for the ith + turbine (m). + + x (np.array): Streamwise direction grid coordinates of the + flow field domain (m). + y (np.array): Cross stream direction grid coordinates of the + flow field domain (m). + z (np.array): Vertical direction grid coordinates of the + flow field domain (m). + wind_veer (np.array): Wind veer (deg). + + Returns: + np.array: Velocity deficits (-). + """ + + # Initial wake widths + + # No specific near, far wakes in this model + downstream_mask = np.array(x > x_i + 0.1) + + # Wake expansion in the lateral (y) and the vertical (z) + # TODO: could compute shared components in sigma_z, sigma_y + # with one function call. + wake_deficit = 0.5 + + velocity_deficit = wake_deficit * downstream_mask + + return velocity_deficit From d3c526b6c3c440755c3724bca3b5e90c42dc5e67 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 22 Jun 2023 14:57:14 +0200 Subject: [PATCH 04/46] Made vawt_blade_length an optional turbine input --- floris/simulation/farm.py | 20 ++++++++++--------- floris/simulation/floris.py | 2 +- floris/simulation/solver.py | 6 +++++- floris/simulation/turbine.py | 3 +-- .../wake_velocity/super_gaussian_vawt.py | 5 +---- 5 files changed, 19 insertions(+), 17 deletions(-) diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index fe2eff6fe..a7b5187d5 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -72,7 +72,7 @@ class Farm(BaseClass): turbine_fCts: tuple = field(init=False, default=[]) turbine_type_map_sorted: NDArrayObject = field(init=False, default=[]) rotor_diameters_sorted: NDArrayFloat = field(init=False, default=[]) - blade_lengths_sorted: NDArrayFloat = field(init=False, default=[]) + vawt_blade_lengths_sorted: NDArrayFloat = field(init=False, default=[]) TSRs_sorted: NDArrayFloat = field(init=False, default=[]) pPs: NDArrayFloat = field(init=False, default=[]) pPs_sorted: NDArrayFloat = field(init=False, default=[]) @@ -210,10 +210,12 @@ def construct_rotor_diameters(self): turb['rotor_diameter'] for turb in self.turbine_definitions ]) - def construct_blade_lengths(self): - self.blade_lengths = np.array([ - turb['blade_length'] for turb in self.turbine_definitions - ]) + def construct_vawt_blade_lengths(self): + for turb in self.turbine_definitions: + try: + self.vawt_blade_lengths = turb['vawt_blade_length'] + except KeyError: + self.vawt_blade_lengths = 0.0 def construct_turbine_TSRs(self): self.TSRs = np.array([turb['TSR'] for turb in self.turbine_definitions]) @@ -277,8 +279,8 @@ def expand_farm_properties( sorted_coord_indices, axis=2 ) - self.blade_lengths_sorted = np.take_along_axis( - self.blade_lengths * template_shape, + self.vawt_blade_lengths_sorted = np.take_along_axis( + self.vawt_blade_lengths * template_shape, sorted_coord_indices, axis=2 ) @@ -372,8 +374,8 @@ def finalize(self, unsorted_indices): unsorted_indices[:,:,:,0,0], axis=2 ) - self.blade_lengths = np.take_along_axis( - self.blade_lengths_sorted, + self.vawt_blade_lengths = np.take_along_axis( + self.vawt_blade_lengths_sorted, unsorted_indices[:,:,:,0,0], axis=2 ) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 8c44c1a81..fc9a8326b 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -76,7 +76,7 @@ def __attrs_post_init__(self) -> None: self.farm.construct_turbine_power_interps() self.farm.construct_hub_heights() self.farm.construct_rotor_diameters() - self.farm.construct_blade_lengths() + self.farm.construct_vawt_blade_lengths() self.farm.construct_turbine_TSRs() self.farm.construct_turbine_pPs() self.farm.construct_turbine_pTs() diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index bc7cf9bef..cc184276f 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -132,7 +132,7 @@ def sequential_solver( yaw_angle_i = farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = farm.hub_heights_sorted[: ,:, i:i+1, None, None] rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] - #blade_length_i = farm.blade_lengths_sorted[: ,:, i:i+1, None, None] + vawt_blade_length_i = farm.vawt_blade_lengths_sorted[: ,:, i:i+1, None, None] TSR_i = farm.TSRs_sorted[: ,:, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) @@ -205,6 +205,7 @@ def sequential_solver( ct_i, hub_height_i, rotor_diameter_i, + vawt_blade_length_i, **deficit_model_args ) @@ -271,6 +272,7 @@ def full_flow_sequential_solver( turbine_grid_farm.construct_turbine_power_interps() turbine_grid_farm.construct_hub_heights() turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_vawt_blade_lengths() turbine_grid_farm.construct_turbine_TSRs() turbine_grid_farm.construct_turbine_pPs() turbine_grid_farm.construct_turbine_pTs() @@ -361,6 +363,7 @@ def full_flow_sequential_solver( yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + vawt_blade_length_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) @@ -421,6 +424,7 @@ def full_flow_sequential_solver( ct_i, hub_height_i, rotor_diameter_i, + vawt_blade_length_i, **deficit_model_args ) diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index 3c16039de..9aec5b31b 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -587,7 +587,6 @@ class Turbine(BaseClass): Parameters: rotor_diameter (:py:obj: float): The rotor diameter (m). - blade_length (:py:obj: float): The blade length (m). hub_height (:py:obj: float): The hub height (m). pP (:py:obj: float): The cosine exponent relating the yaw misalignment angle to power. @@ -618,7 +617,6 @@ class Turbine(BaseClass): turbine_type: str = field() rotor_diameter: float = field() - blade_length: float = field() hub_height: float = field() pP: float = field() pT: float = field() @@ -629,6 +627,7 @@ class Turbine(BaseClass): power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) floating_tilt_table = field(default=None) floating_correct_cp_ct_for_tilt = field(default=None) + vawt_blade_length: float = field(default=None) # rloc: float = float_attrib() # TODO: goes here or on the Grid? # use_points_on_perimeter: bool = bool_attrib() diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index fe3442097..1e274b72e 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -90,12 +90,11 @@ def function( axial_induction_i: np.ndarray, deflection_field_i: np.ndarray, yaw_angle_i: np.ndarray, - tilt_angle_i: np.ndarray, turbulence_intensity_i: np.ndarray, ct_i: np.ndarray, hub_height_i: float, rotor_diameter_i: np.ndarray, - #blade_length_i: np.ndarray, + vawt_blade_length_i: np.ndarray, # enforces the use of the below as keyword arguments and adherence to the # unpacking of the results from prepare_function() *, @@ -117,8 +116,6 @@ def function( axial_induction_i (np.array): Axial induction factor of the ith turbine (-) [not used]. yaw_angle_i (np.array): Yaw angle of the ith turbine (deg). - tilt_angle_i (np.array): Tilt angle of the ith turbine - (deg). ct_i (np.array): Thrust coefficient for the ith turbine (-). hub_height_i (float): Hub height for the ith turbine (m). rotor_diameter_i (np.array): Rotor diameter for the ith From 2ed197ed788887b6171f26c26de6b7c3f983b046 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 22 Jun 2023 15:15:04 +0200 Subject: [PATCH 05/46] Bug fix --- floris/simulation/solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index cc184276f..330c3a17a 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -363,7 +363,7 @@ def full_flow_sequential_solver( yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] - vawt_blade_length_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + vawt_blade_length_i = turbine_grid_farm.vawt_blade_lengths_sorted[:, :, i:i+1, None, None] TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) From fb3a75dfcc63ef072a1b4fff4b204df8a5e58efb Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 22 Jun 2023 21:18:50 +0200 Subject: [PATCH 06/46] Untested draft of super_gaussian_vawt wake_velocity model --- .../wake_velocity/super_gaussian_vawt.py | 62 ++++++++++++++++--- pyproject.toml | 1 + 2 files changed, 53 insertions(+), 10 deletions(-) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 1e274b72e..7b14f35de 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -15,6 +15,7 @@ import numexpr as ne import numpy as np from attrs import define, field +from scipy.special import gamma from floris.simulation import ( BaseModel, @@ -133,16 +134,57 @@ def function( np.array: Velocity deficits (-). """ - # Initial wake widths + # Model parameters need to be unpacked for use in Numexpr + ay, by, cy = self.ay, self.by, self.cy + az, bz, cz = self.az, self.bz, self.cz - # No specific near, far wakes in this model + # No specific near or far wakes in this model downstream_mask = np.array(x > x_i + 0.1) - # Wake expansion in the lateral (y) and the vertical (z) - # TODO: could compute shared components in sigma_z, sigma_y - # with one function call. - wake_deficit = 0.5 - - velocity_deficit = wake_deficit * downstream_mask - - return velocity_deficit + # Nondimensional coordinates. Is z_i always equal to hub_height_i? + x_tilde = ne.evaluate("downstream_mask * (x - x_i) / rotor_diameter_i") + x_tilde_H = ne.evaluate("downstream_mask * (x - x_i) / vawt_blade_length_i") + y_tilde = ne.evaluate("downstream_mask * (y - y_i) / rotor_diameter_i") + z_tilde = ne.evaluate("downstream_mask * (z - z_i) / vawt_blade_length_i") + + # Wake expansion rates + ky_star = self.wake_expansion_coeff_y * turbulence_intensity_i + kz_star = self.wake_expansion_coeff_z * turbulence_intensity_i + + # Relative initial wake expansion + fac = np.sqrt(1 - ct_i) + beta = 0.5 * (1 + fac) / fac + # Initial wake width + # Obs: Change ny nz to have x_tilde + ny_0 = ay + cy + nz_0 = az + cz + eta_0 = 1/ny_0 + 1/nz_0 + epsilon = beta * ny_0 * nz_0 / (2 ** (2 * eta_0 + 2) * gamma(1 / ny_0) * gamma(1 / nz_0)) + epsilon **= 1 / (2 * eta_0) + epsilon_y, epsilon_z = epsilon, epsilon + + # Characteristic wake widths grow linearly in the streamwise direction + sigma_y_tilde = ne.evaluate("ky_star * x_tilde + epsilon_y") + sigma_z_tilde = ne.evaluate("kz_star * x_tilde_H + epsilon_z") + + # Gaussian exponents which determine the wake shape + ny = ne.evaluate("ay * exp(-by * x_tilde) + cy") + nz = ne.evaluate("az * exp(-bz * x_tilde_H) + cz") + + # Velocity deficit given by the maximum velocity C at a given streamwise distance, + # and two Gaussian shape functions fy and fz + eta = 1 / ny + 1 / nz + ny_inv = ne.evaluate("1 / ny") + nz_inv = ne.evaluate("1 / nz") + gamma_y, gamma_z = gamma(ny_inv), gamma(nz_inv) + fac = ne.evaluate( + "8 * sigma_y_tilde ** (2 / ny) * sigma_z_tilde ** (2 / nz) * gamma_y * gamma_z" + ) + C = ne.evaluate( + "2 ** (eta - 1) - sqrt(2 ** (2 * eta - 2) - ct_i * ny * nz / fac)" + ) + + fy = ne.evaluate("exp(-y_tilde ** ny / (2 * sigma_y_tilde ** 2))") + fz = ne.evaluate("exp(-z_tilde ** nz / (2 * sigma_z_tilde ** 2))") + + return ne.evaluate("C * fy * fz * downstream_mask") diff --git a/pyproject.toml b/pyproject.toml index 39683f439..220dbb00d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -120,6 +120,7 @@ dummy-variable-rgx = "^(_+|(_+[a-zA-Z0-9_]*[a-zA-Z0-9]+?))$" "floris/simulation/wake_velocity/jensen.py" = ["F841"] "floris/simulation/wake_velocity/gauss.py" = ["F841"] "floris/simulation/wake_velocity/empirical_gauss.py" = ["F841"] +"floris/simulation/wake_velocity/super_gaussian_vawt.py" = ["F841"] # I001 unsorted-imports: ignore because the import order is meaningful to navigate # import dependencies From da7fb6dace2b4787dec20590d51c2e883b23dc48 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 25 Jun 2023 19:28:58 +0200 Subject: [PATCH 07/46] Can't have negative coordinates to the power of a float --- floris/simulation/wake_velocity/super_gaussian_vawt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 7b14f35de..0cf5d77dc 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -184,7 +184,7 @@ def function( "2 ** (eta - 1) - sqrt(2 ** (2 * eta - 2) - ct_i * ny * nz / fac)" ) - fy = ne.evaluate("exp(-y_tilde ** ny / (2 * sigma_y_tilde ** 2))") - fz = ne.evaluate("exp(-z_tilde ** nz / (2 * sigma_z_tilde ** 2))") + fy = ne.evaluate("exp(-abs(y_tilde) ** ny / (2 * sigma_y_tilde ** 2))") + fz = ne.evaluate("exp(-abs(z_tilde) ** nz / (2 * sigma_z_tilde ** 2))") return ne.evaluate("C * fy * fz * downstream_mask") From af9e1b08d205e8b18841936f42f4ee0f14970c90 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 25 Jun 2023 20:12:41 +0200 Subject: [PATCH 08/46] Comments --- floris/simulation/wake_velocity/super_gaussian_vawt.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 0cf5d77dc..e93d1c03d 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -154,8 +154,8 @@ def function( # Relative initial wake expansion fac = np.sqrt(1 - ct_i) beta = 0.5 * (1 + fac) / fac + # Initial wake width - # Obs: Change ny nz to have x_tilde ny_0 = ay + cy nz_0 = az + cz eta_0 = 1/ny_0 + 1/nz_0 @@ -171,8 +171,9 @@ def function( ny = ne.evaluate("ay * exp(-by * x_tilde) + cy") nz = ne.evaluate("az * exp(-bz * x_tilde_H) + cz") - # Velocity deficit given by the maximum velocity C at a given streamwise distance, - # and two Gaussian shape functions fy and fz + # At any streamwise location x behind the turbine, the velocity deficit in the + # y-z plane is given by multiplying the maximum deficit C = C(x) with two Gaussian + # shape functions fy = fy(y) and fz = fz(z) eta = 1 / ny + 1 / nz ny_inv = ne.evaluate("1 / ny") nz_inv = ne.evaluate("1 / nz") From eeb0b988d83fe569314e4cdf6704736b7697016f Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 26 Jun 2023 13:04:14 +0200 Subject: [PATCH 09/46] Created a basic vawt-solver with no yaw-related features and no turbulence model --- floris/simulation/__init__.py | 2 + floris/simulation/floris.py | 17 +- floris/simulation/solver.py | 219 +++++++++++++++++- floris/simulation/turbine.py | 2 +- .../wake_velocity/super_gaussian_vawt.py | 3 - 5 files changed, 233 insertions(+), 10 deletions(-) diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 12d30aab8..3d2574e49 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -52,8 +52,10 @@ from .solver import ( cc_solver, empirical_gauss_solver, + vawt_solver, full_flow_cc_solver, full_flow_empirical_gauss_solver, + full_flow_vawt_solver, full_flow_sequential_solver, full_flow_turbopark_solver, sequential_solver, diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index fc9a8326b..8247ec33c 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -32,6 +32,7 @@ full_flow_empirical_gauss_solver, full_flow_sequential_solver, full_flow_turbopark_solver, + full_flow_vawt_solver, Grid, PointsGrid, sequential_solver, @@ -39,6 +40,7 @@ TurbineCubatureGrid, TurbineGrid, turbopark_solver, + vawt_solver, WakeModelManager, ) from floris.utilities import load_yaml @@ -215,7 +217,7 @@ def steady_state_atmospheric_condition(self): # <> # start = time.time() - if vel_model in ["gauss", "cc", "turbopark", "jensen"] and \ + if vel_model in ["gauss", "cc", "turbopark", "jensen", "super_gaussian_vawt"] and \ self.farm.correct_cp_ct_for_tilt.any(): self.logger.warn( "The current model does not account for vertical wake deflection due to " + @@ -244,6 +246,13 @@ def steady_state_atmospheric_condition(self): self.grid, self.wake ) + elif vel_model=="super_gaussian_vawt": + vawt_solver( + self.farm, + self.flow_field, + self.grid, + self.wake + ) else: sequential_solver( self.farm, @@ -274,6 +283,8 @@ def solve_for_viz(self): full_flow_turbopark_solver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, self.grid, self.wake) + elif vel_model=="super_gaussian_vawt": + full_flow_vawt_solver(self.farm, self.flow_field, self.grid, self.wake) else: full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) @@ -306,10 +317,12 @@ def solve_for_points(self, x, y, z): if vel_model == "cc" or vel_model == "turbopark": raise NotImplementedError( "solve_for_points is currently only available with the "+\ - "gauss, jensen, and empirical_guass models." + "gauss, jensen, empirical_guass, and super_gaussian_vawt models." ) elif vel_model == "empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) + elif vel_model == "super_gaussian_vawt": + full_flow_vawt_solver(self.farm, self.flow_field, field_grid, self.wake) else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 330c3a17a..b2a53c0ec 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -132,7 +132,6 @@ def sequential_solver( yaw_angle_i = farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = farm.hub_heights_sorted[: ,:, i:i+1, None, None] rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] - vawt_blade_length_i = farm.vawt_blade_lengths_sorted[: ,:, i:i+1, None, None] TSR_i = farm.TSRs_sorted[: ,:, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) @@ -205,7 +204,6 @@ def sequential_solver( ct_i, hub_height_i, rotor_diameter_i, - vawt_blade_length_i, **deficit_model_args ) @@ -363,7 +361,6 @@ def full_flow_sequential_solver( yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None] hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] - vawt_blade_length_i = turbine_grid_farm.vawt_blade_lengths_sorted[:, :, i:i+1, None, None] TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) @@ -424,7 +421,6 @@ def full_flow_sequential_solver( ct_i, hub_height_i, rotor_diameter_i, - vawt_blade_length_i, **deficit_model_args ) @@ -1499,3 +1495,218 @@ def full_flow_empirical_gauss_solver( flow_field.u_sorted = flow_field.u_initial_sorted - wake_field flow_field.v_sorted += v_wake flow_field.w_sorted += w_wake + + + +def vawt_solver( + farm: Farm, + flow_field: FlowField, + grid: TurbineGrid, + model_manager: WakeModelManager +) -> None: + # Algorithm + # For each turbine, calculate its effect on every downstream turbine. + # For the current turbine, we are calculating the deficit that it adds to downstream turbines. + # Integrate this into the main data structure. + # Move on to the next turbine. + # + # Note that no turbulence model has been implmented in this solver. + + # <> + deficit_model_args = model_manager.velocity_model.prepare_function(grid, flow_field) + + # This is u_wake + wake_field = np.zeros_like(flow_field.u_initial_sorted) + + turbine_turbulence_intensity = ( + flow_field.turbulence_intensity + * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1)) + ) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for i in range(grid.n_turbines): + + # Get the current turbine quantities + x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4)) + x_i = x_i[:, :, :, None, None] + y_i = np.mean(grid.y_sorted[:, :, i:i+1], axis=(3, 4)) + y_i = y_i[:, :, :, None, None] + z_i = np.mean(grid.z_sorted[:, :, i:i+1], axis=(3, 4)) + z_i = z_i[:, :, :, None, None] + + ct_i = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + tilt_angle=farm.tilt_angles_sorted, + ref_tilt_cp_ct=farm.ref_tilt_cp_cts_sorted, + fCt=farm.turbine_fCts, + tilt_interp=farm.turbine_fTilts, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights + ) + # Since we are filtering for the i'th turbine in the Ct function, + # get the first index here (0:1) + ct_i = ct_i[:, :, 0:1, None, None] + + turbulence_intensity_i = turbine_turbulence_intensity[:, :, i:i+1] + hub_height_i = farm.hub_heights_sorted[: ,:, i:i+1, None, None] + rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] + vawt_blade_length_i = farm.vawt_blade_lengths_sorted[: ,:, i:i+1, None, None] + + if model_manager.enable_secondary_steering: + raise NotImplementedError( + "Secondary steering not available for this model.") + + if model_manager.enable_transverse_velocities: + raise NotImplementedError( + "Transverse velocities not available for this model.") + + if model_manager.enable_yaw_added_recovery: + raise NotImplementedError( + "Yaw added recovery not available for this model.") + + # Model calculations + # NOTE: exponential + velocity_deficit = model_manager.velocity_model.function( + x_i, + y_i, + z_i, + turbulence_intensity_i, + ct_i, + hub_height_i, + rotor_diameter_i, + vawt_blade_length_i, + **deficit_model_args + ) + + wake_field = model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + flow_field.u_sorted = flow_field.u_initial_sorted - wake_field + + flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity + flow_field.turbulence_intensity_field_sorted_avg = np.mean( + turbine_turbulence_intensity, + axis=(3,4) + )[:, :, :, None, None] + + +def full_flow_vawt_solver( + farm: Farm, + flow_field: FlowField, + flow_field_grid: FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid, + model_manager: WakeModelManager +) -> None: + + # Get the flow quantities and turbine performance + turbine_grid_farm = copy.deepcopy(farm) + turbine_grid_flow_field = copy.deepcopy(flow_field) + + turbine_grid_farm.construct_turbine_map() + turbine_grid_farm.construct_turbine_fCts() + turbine_grid_farm.construct_turbine_power_interps() + turbine_grid_farm.construct_hub_heights() + turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_vawt_blade_lengths() + turbine_grid_farm.construct_turbine_TSRs() + turbine_grid_farm.construct_turbine_pPs() + turbine_grid_farm.construct_turbine_pTs() + turbine_grid_farm.construct_turbine_ref_density_cp_cts() + turbine_grid_farm.construct_turbine_ref_tilt_cp_cts() + turbine_grid_farm.construct_turbine_fTilts() + turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() + turbine_grid_farm.construct_coordinates() + turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) + + turbine_grid = TurbineGrid( + turbine_coordinates=turbine_grid_farm.coordinates, + reference_turbine_diameter=turbine_grid_farm.rotor_diameters, + wind_directions=turbine_grid_flow_field.wind_directions, + wind_speeds=turbine_grid_flow_field.wind_speeds, + grid_resolution=3, + time_series=turbine_grid_flow_field.time_series, + ) + turbine_grid_farm.expand_farm_properties( + turbine_grid_flow_field.n_wind_directions, + turbine_grid_flow_field.n_wind_speeds, + turbine_grid.sorted_coord_indices + ) + turbine_grid_flow_field.initialize_velocity_field(turbine_grid) + turbine_grid_farm.initialize(turbine_grid.sorted_indices) + vawt_solver(turbine_grid_farm, turbine_grid_flow_field, turbine_grid, model_manager) + + ### Referring to the quantities from above, calculate the wake in the full grid + + # Use full flow_field here to use the full grid in the wake models + deficit_model_args = model_manager.velocity_model.prepare_function( + flow_field_grid, + flow_field + ) + + wake_field = np.zeros_like(flow_field.u_initial_sorted) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for i in range(flow_field_grid.n_turbines): + + # Get the current turbine quantities + x_i = np.mean(turbine_grid.x_sorted[:, :, i:i+1], axis=(3, 4)) + x_i = x_i[:, :, :, None, None] + y_i = np.mean(turbine_grid.y_sorted[:, :, i:i+1], axis=(3, 4)) + y_i = y_i[:, :, :, None, None] + z_i = np.mean(turbine_grid.z_sorted[:, :, i:i+1], axis=(3, 4)) + z_i = z_i[:, :, :, None, None] + + ct_i = Ct( + velocities=turbine_grid_flow_field.u_sorted, + yaw_angle=turbine_grid_farm.yaw_angles_sorted, + tilt_angle=turbine_grid_farm.tilt_angles_sorted, + ref_tilt_cp_ct=turbine_grid_farm.ref_tilt_cp_cts_sorted, + fCt=turbine_grid_farm.turbine_fCts, + tilt_interp=turbine_grid_farm.turbine_fTilts, + correct_cp_ct_for_tilt=turbine_grid_farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=turbine_grid_farm.turbine_type_map_sorted, + ix_filter=[i], + ) + # Since we are filtering for the i'th turbine in the Ct function, + # get the first index here (0:1) + ct_i = ct_i[:, :, 0:1, None, None] + + turbulence_intensity_i = \ + turbine_grid_flow_field.turbulence_intensity_field_sorted_avg[:, :, i:i+1] + hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] + rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + vawt_blade_length_i = turbine_grid_farm.vawt_blade_lengths_sorted[:, :, i:i+1, None, None] + + if model_manager.enable_secondary_steering: + raise NotImplementedError( + "Secondary steering not available for this model.") + + if model_manager.enable_transverse_velocities: + raise NotImplementedError( + "Transverse velocities not available for this model.") + + # Model calculations + # NOTE: exponential + velocity_deficit = model_manager.velocity_model.function( + x_i, + y_i, + z_i, + turbulence_intensity_i, + ct_i, + hub_height_i, + rotor_diameter_i, + vawt_blade_length_i, + **deficit_model_args + ) + + wake_field = model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + flow_field.u_sorted = flow_field.u_initial_sorted - wake_field diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index 9aec5b31b..604cb6ef8 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -627,7 +627,7 @@ class Turbine(BaseClass): power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) floating_tilt_table = field(default=None) floating_correct_cp_ct_for_tilt = field(default=None) - vawt_blade_length: float = field(default=None) + vawt_blade_length = field(default=None) # rloc: float = float_attrib() # TODO: goes here or on the Grid? # use_points_on_perimeter: bool = bool_attrib() diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index e93d1c03d..f879e20aa 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -88,9 +88,6 @@ def function( x_i: np.ndarray, y_i: np.ndarray, z_i: np.ndarray, - axial_induction_i: np.ndarray, - deflection_field_i: np.ndarray, - yaw_angle_i: np.ndarray, turbulence_intensity_i: np.ndarray, ct_i: np.ndarray, hub_height_i: float, From c7c0fe37cd1e7bdeb21110ca3cc490b7583610c8 Mon Sep 17 00:00:00 2001 From: vallbog Date: Wed, 28 Jun 2023 12:04:46 +0200 Subject: [PATCH 10/46] Moved from other branch --- floris/simulation/__init__.py | 1 + floris/simulation/floris.py | 72 ++++++++++++++++++++++++++++++++ floris/simulation/grid.py | 44 ++++++++++++++++++- floris/tools/floris_interface.py | 44 +++++++++++++++++++ 4 files changed, 160 insertions(+), 1 deletion(-) diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 12d30aab8..c87f39d6e 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -44,6 +44,7 @@ FlowFieldPlanarGrid, Grid, PointsGrid, + VelocityProfileGrid, TurbineGrid, TurbineCubatureGrid ) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 8ca8854e7..00bf7ef25 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -39,6 +39,7 @@ TurbineCubatureGrid, TurbineGrid, turbopark_solver, + VelocityProfileGrid, WakeModelManager, ) from floris.utilities import load_yaml @@ -314,6 +315,77 @@ def solve_for_points(self, x, y, z): return self.flow_field.u_sorted[:,:,:,0,0] # Remove turbine grid dimensions + def solve_for_velocity_deficit_profiles( + self, + direction, + downstream_dists, + profile_range, + resolution, + ref_turbine_diameter, + x_inertial_start, + y_inertial_start + ): + +# field_grid = VelocityProfileGrid( +# direction, +# downstream_dists, +# profile_range, +# resolution, +# ref_turbine_diameter, +# x_inertial_start, +# y_inertial_start, +# self.grid.x_center_of_rotation, +# self.grid.y_center_of_rotation, +# turbine_coordinates=self.farm.coordinates, +# reference_turbine_diameter=self.farm.rotor_diameters, +# grid_resolution=1, +# wind_directions=self.flow_field.wind_directions, +# wind_speeds=self.flow_field.wind_speeds, +# time_series=self.flow_field.time_series +# ) + +# Apparently need to use keywords if multi-init + VelocityProfileGrid( + direction=direction, + downstream_dists=downstream_dists, + profile_range=profile_range, + resolution=resolution, + ref_turbine_diameter=ref_turbine_diameter, + x_inertial_start=x_inertial_start, + y_inertial_start=y_inertial_start, + x_center_of_rotation=self.grid.x_center_of_rotation, + y_center_of_rotation=self.grid.y_center_of_rotation, + turbine_coordinates=self.farm.coordinates, + reference_turbine_diameter=self.farm.rotor_diameters, + grid_resolution=1, + wind_directions=self.flow_field.wind_directions, + wind_speeds=self.flow_field.wind_speeds, + time_series=self.flow_field.time_series + ) + + + +# field_grid = VelocityProfileGrid( +# turbine_coordinates=self.farm.coordinates, +# reference_turbine_diameter=self.farm.rotor_diameters, +# wind_directions=self.flow_field.wind_directions, +# wind_speeds=self.flow_field.wind_speeds, +# grid_resolution=1, +# time_series=self.flow_field.time_series, +# x_center_of_rotation=self.grid.x_center_of_rotation, +# y_center_of_rotation=self.grid.y_center_of_rotation +# ) + +# Worked with minimal example of VelocityProfileGrid +# field_grid = VelocityProfileGrid( +# turbine_coordinates=self.farm.coordinates, +# reference_turbine_diameter=self.farm.rotor_diameters, +# grid_resolution=1, +# wind_directions=self.flow_field.wind_directions, +# wind_speeds=self.flow_field.wind_speeds, +# time_series=self.flow_field.time_series, +# ) + def finalize(self): # Once the wake calculation is finished, unsort the values to match # the user-supplied order of things. diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 8e2325252..f9fef2408 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -86,6 +86,7 @@ class Grid(ABC): cubature_weights: NDArrayFloat = field(init=False, default=None) def __attrs_post_init__(self) -> None: + print(self.grid_resolution, type(self.grid_resolution)) self.turbine_coordinates_array = np.array([c.elements for c in self.turbine_coordinates]) @turbine_coordinates.validator @@ -118,7 +119,7 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter # TODO move this to the grid types and off of the base class """Check that grid resolution is given as int or Vec3 with int components.""" if isinstance(value, int) and \ - isinstance(self, (TurbineGrid, TurbineCubatureGrid, PointsGrid)): + isinstance(self, (TurbineGrid, TurbineCubatureGrid, PointsGrid, VelocityProfileGrid)): return elif isinstance(value, Iterable) and isinstance(self, FlowFieldPlanarGrid): assert type(value[0]) is int @@ -744,3 +745,44 @@ def set_grid(self) -> None: self.x_sorted = x[:,:,:,None,None] self.y_sorted = y[:,:,:,None,None] self.z_sorted = z[:,:,:,None,None] + +@define +class VelocityProfileGrid(Grid): + """ + docstr + """ + direction: str + downstream_dists: NDArrayFloat + profile_range: NDArrayFloat + resolution: float + ref_turbine_diameter: float + x_inertial_start: float + y_inertial_start: float + x_center_of_rotation: float | None = field(default=None) + y_center_of_rotation: float | None = field(default=None) + + def __attrs_post_init__(self) -> None: + super().__attrs_post_init__() + self.set_grid() + + def set_grid(self) -> None: + """ + Set points for calculation based on a series of user-supplied coordinates. + """ + print(self.downstream_dists) + +#@define +#class VelocityProfileGrid(Grid): +## x_center_of_rotation: float | None = field(default=None) +## y_center_of_rotation: float | None = field(default=None) +# +# def __attrs_post_init__(self) -> None: +# super().__attrs_post_init__() +# self.set_grid() +# +# def set_grid(self) -> None: +# """ +# Set points for calculation based on a series of user-supplied coordinates. +# """ +# print("Hi") +# # print(self.x_center_of_rotation) diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 4e4cc352e..29c6189af 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -950,6 +950,50 @@ def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloa return self.floris.solve_for_points(x, y, z) + def sample_velocity_deficit_profiles( + self, + direction='y', + downstream_dists=None, + profile_range=None, + resolution=100, + ref_turbine_diameter=None, + x_inertial_start=None, + y_inertial_start=None + ): + """ + Extract velocity deficit profiles + + Args: + x (1DArrayFloat | list): x-locations of points where flow is desired. + y (1DArrayFloat | list): y-locations of points where flow is desired. + z (1DArrayFloat | list): z-locations of points where flow is desired. + + Returns: + 3DArrayFloat containing wind speed with dimensions + (# of wind directions, # of wind speeds, # of sample points) + """ + + if downstream_dists is None: + downstream_dists = ref_turbine_diameter * np.array([3, 5, 7, 9]) + + if profile_range is None: + profile_range = ref_turbine_diameter * np.array([-2, 2]) + + if direction not in ('y', 'z'): + raise ValueError("'direction' must be either y or z") + + velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( + direction, + downstream_dists, + profile_range, + resolution, + ref_turbine_diameter, + x_inertial_start, + y_inertial_start + ) + + return velocity_deficit_profiles + @property def layout_x(self): """ From 6688ee29eddc831bfaf1b0aee6d76ce7fa88dca6 Mon Sep 17 00:00:00 2001 From: vallbog Date: Wed, 28 Jun 2023 19:24:37 +0200 Subject: [PATCH 11/46] Preparing inputs in FlorisInterface --- floris/simulation/grid.py | 6 +-- floris/tools/floris_interface.py | 90 ++++++++++++++++++++++++++------ 2 files changed, 76 insertions(+), 20 deletions(-) diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index f9fef2408..9328c0d40 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -752,9 +752,9 @@ class VelocityProfileGrid(Grid): docstr """ direction: str - downstream_dists: NDArrayFloat - profile_range: NDArrayFloat - resolution: float + downstream_dists: NDArrayFloat = field(converter=floris_array_converter) + profile_range: NDArrayFloat = field(converter=floris_array_converter) + resolution: int ref_turbine_diameter: float x_inertial_start: float y_inertial_start: float diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 29c6189af..f21535dae 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -952,47 +952,103 @@ def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloa def sample_velocity_deficit_profiles( self, - direction='y', - downstream_dists=None, - profile_range=None, + direction: NDArrayFloat = 'y', + downstream_dists: NDArrayFloat | list = None, + profile_range: NDArrayFloat | list = None, resolution=100, - ref_turbine_diameter=None, - x_inertial_start=None, - y_inertial_start=None + wind_direction: float = None, + homogeneous_wind_speed: float = None, + ref_rotor_diameter: float = None, + x_inertial_start: float = None, + y_inertial_start: float = None, + reference_height: float = None ): """ Extract velocity deficit profiles Args: - x (1DArrayFloat | list): x-locations of points where flow is desired. - y (1DArrayFloat | list): y-locations of points where flow is desired. - z (1DArrayFloat | list): z-locations of points where flow is desired. + downstream_dists: Streamwise location for each velocity profile. Default starting point + is at (x, y) = (0, 0). + + reference_height: If 'direction' is y, then 'reference_height' defines the height of the + xy-plane in which the velocity profiles are sampled. If 'direction' is z, then the + velocity is sampled along the vertical direction with the 'profile_range' being + relative to the 'reference_height'. Returns: 3DArrayFloat containing wind speed with dimensions (# of wind directions, # of wind speeds, # of sample points) """ + if direction not in ('y', 'z'): + raise ValueError("'direction' must be either y or z") + + if ref_rotor_diameter is None: + rotor_diameters = self.floris.farm.rotor_diameters + if np.all(rotor_diameters == rotor_diameters[0]): + ref_rotor_diameter = rotor_diameters[0] + else: + raise ValueError( + "'ref_rotor_diameter' needs to be defined manually since there are multiple" + "turbine types." + ) + if downstream_dists is None: - downstream_dists = ref_turbine_diameter * np.array([3, 5, 7, 9]) + downstream_dists = ref_rotor_diameter * np.array([3, 5, 7, 9]) if profile_range is None: - profile_range = ref_turbine_diameter * np.array([-2, 2]) + profile_range = ref_rotor_diameter * np.array([-2, 2]) - if direction not in ('y', 'z'): - raise ValueError("'direction' must be either y or z") + wind_directions_copy = np.array(self.floris.flow_field.wind_directions, copy=True) + wind_speeds_copy = np.array(self.floris.flow_field.wind_speeds, copy=True) + + if wind_direction is None: + if len(wind_directions_copy) == 1: + wind_direction = wind_directions_copy[0] + else: + raise ValueError( + "Multiple wind directions detected. Specify a single 'wind_direction' for" + "which to sample the velocity profiles." + ) - velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( + if homogeneous_wind_speed is None: + if len(wind_speeds_copy) == 1: + homogeneous_wind_speed = wind_speeds_copy[0] + # Maybe add info msg that wind_speed is homogeneous + else: + raise ValueError( + "Multiple wind speeds detected. Specify a single 'homogeneous_wind_speed' for" + "which to sample the velocity profiles." + ) + + if reference_height is None: + reference_height = self.floris.flow_field.reference_wind_height + + print( direction, downstream_dists, profile_range, resolution, - ref_turbine_diameter, + wind_direction, + homogeneous_wind_speed, + ref_rotor_diameter, x_inertial_start, - y_inertial_start + y_inertial_start, + reference_height ) - return velocity_deficit_profiles + +# velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( +# direction, +# downstream_dists, +# profile_range, +# resolution, +# ref_turbine_diameter, +# x_inertial_start, +# y_inertial_start +# ) +# +# return velocity_deficit_profiles @property def layout_x(self): From 55281b636fede946fc55709642fa74dd0751e126 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 29 Jun 2023 15:31:31 +0200 Subject: [PATCH 12/46] Created VelocityProfileGrid --- examples/29_plot_velocity_deficit_profiles.py | 40 ++++++++++++ floris/simulation/floris.py | 48 ++------------ floris/simulation/grid.py | 63 +++++++++++++------ floris/tools/floris_interface.py | 50 +++++++++++---- 4 files changed, 126 insertions(+), 75 deletions(-) create mode 100644 examples/29_plot_velocity_deficit_profiles.py diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py new file mode 100644 index 000000000..f89deabe4 --- /dev/null +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -0,0 +1,40 @@ +# 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 + + +""" +docstr +""" + +fi = FlorisInterface("inputs/gch.yaml") +fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) + +fi.sample_velocity_deficit_profiles(direction = 'y', resolution=10, homogeneous_wind_speed=7.0) + +#horizontal_plane = fi.calculate_horizontal_plane( +# x_resolution=200, +# y_resolution=100, +# height=90.0, +#) +# +#wakeviz.visualize_cut_plane(horizontal_plane, title="Horizontal plane") +# +#wakeviz.show_plots() diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 00bf7ef25..381b50775 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -321,38 +321,21 @@ def solve_for_velocity_deficit_profiles( downstream_dists, profile_range, resolution, - ref_turbine_diameter, + ref_rotor_diameter, x_inertial_start, - y_inertial_start + y_inertial_start, + reference_height ): -# field_grid = VelocityProfileGrid( -# direction, -# downstream_dists, -# profile_range, -# resolution, -# ref_turbine_diameter, -# x_inertial_start, -# y_inertial_start, -# self.grid.x_center_of_rotation, -# self.grid.y_center_of_rotation, -# turbine_coordinates=self.farm.coordinates, -# reference_turbine_diameter=self.farm.rotor_diameters, -# grid_resolution=1, -# wind_directions=self.flow_field.wind_directions, -# wind_speeds=self.flow_field.wind_speeds, -# time_series=self.flow_field.time_series -# ) - -# Apparently need to use keywords if multi-init VelocityProfileGrid( direction=direction, downstream_dists=downstream_dists, profile_range=profile_range, resolution=resolution, - ref_turbine_diameter=ref_turbine_diameter, + ref_rotor_diameter=ref_rotor_diameter, x_inertial_start=x_inertial_start, y_inertial_start=y_inertial_start, + reference_height=reference_height, x_center_of_rotation=self.grid.x_center_of_rotation, y_center_of_rotation=self.grid.y_center_of_rotation, turbine_coordinates=self.farm.coordinates, @@ -365,27 +348,6 @@ def solve_for_velocity_deficit_profiles( -# field_grid = VelocityProfileGrid( -# turbine_coordinates=self.farm.coordinates, -# reference_turbine_diameter=self.farm.rotor_diameters, -# wind_directions=self.flow_field.wind_directions, -# wind_speeds=self.flow_field.wind_speeds, -# grid_resolution=1, -# time_series=self.flow_field.time_series, -# x_center_of_rotation=self.grid.x_center_of_rotation, -# y_center_of_rotation=self.grid.y_center_of_rotation -# ) - -# Worked with minimal example of VelocityProfileGrid -# field_grid = VelocityProfileGrid( -# turbine_coordinates=self.farm.coordinates, -# reference_turbine_diameter=self.farm.rotor_diameters, -# grid_resolution=1, -# wind_directions=self.flow_field.wind_directions, -# wind_speeds=self.flow_field.wind_speeds, -# time_series=self.flow_field.time_series, -# ) - def finalize(self): # Once the wake calculation is finished, unsort the values to match # the user-supplied order of things. diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 9328c0d40..3a93d8746 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -86,7 +86,6 @@ class Grid(ABC): cubature_weights: NDArrayFloat = field(init=False, default=None) def __attrs_post_init__(self) -> None: - print(self.grid_resolution, type(self.grid_resolution)) self.turbine_coordinates_array = np.array([c.elements for c in self.turbine_coordinates]) @turbine_coordinates.validator @@ -755,9 +754,10 @@ class VelocityProfileGrid(Grid): downstream_dists: NDArrayFloat = field(converter=floris_array_converter) profile_range: NDArrayFloat = field(converter=floris_array_converter) resolution: int - ref_turbine_diameter: float + ref_rotor_diameter: float x_inertial_start: float y_inertial_start: float + reference_height: float x_center_of_rotation: float | None = field(default=None) y_center_of_rotation: float | None = field(default=None) @@ -769,20 +769,45 @@ def set_grid(self) -> None: """ Set points for calculation based on a series of user-supplied coordinates. """ - print(self.downstream_dists) - -#@define -#class VelocityProfileGrid(Grid): -## x_center_of_rotation: float | None = field(default=None) -## y_center_of_rotation: float | None = field(default=None) -# -# def __attrs_post_init__(self) -> None: -# super().__attrs_post_init__() -# self.set_grid() -# -# def set_grid(self) -> None: -# """ -# Set points for calculation based on a series of user-supplied coordinates. -# """ -# print("Hi") -# # print(self.x_center_of_rotation) + res = self.resolution + nProfiles = len(self.downstream_dists) + # downsteam_dists are defined from the following starting point + coordinates_inertial_start = np.array( + [[self.x_inertial_start, self.y_inertial_start, self.reference_height]] + ) + + x_start, y_start, _, _, _ = rotate_coordinates_rel_west( + self.wind_directions, + coordinates_inertial_start, + x_center_of_rotation=self.x_center_of_rotation, + y_center_of_rotation=self.y_center_of_rotation + ) + x_start, y_start = x_start[0,0,0], y_start[0,0,0] + + downstream_dists_transpose = np.atleast_2d(self.downstream_dists).T + x = (x_start + downstream_dists_transpose) * np.ones((nProfiles, res)) + + if self.direction == 'y': + y_single_profile = np.linspace( + y_start + self.profile_range[0], + y_start + self.profile_range[1], + res + ) + y = y_single_profile * np.ones((nProfiles, res)) + z = self.reference_height * np.ones((nProfiles, res)) + elif self.direction == 'z': + z_min_profile = self.reference_height + self.profile_range[0] + if z_min_profile <= 0.0: + # Arbitrary small value to avoid possible errors at z = 0.0 + z_min_profile = 0.1 + z_single_profile = np.linspace( + z_min_profile, + self.reference_height + self.profile_range[1], + res + ) + z = z_single_profile * np.ones((nProfiles, res)) + y = y_start * np.ones((nProfiles, res)) + + self.x_sorted = x.flatten()[None,None,:,None,None] + self.y_sorted = y.flatten()[None,None,:,None,None] + self.z_sorted = z.flatten()[None,None,:,None,None] diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index f21535dae..fc4814ea0 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -1001,6 +1001,7 @@ def sample_velocity_deficit_profiles( wind_directions_copy = np.array(self.floris.flow_field.wind_directions, copy=True) wind_speeds_copy = np.array(self.floris.flow_field.wind_speeds, copy=True) + wind_shear_copy = self.floris.flow_field.wind_shear if wind_direction is None: if len(wind_directions_copy) == 1: @@ -1021,33 +1022,56 @@ def sample_velocity_deficit_profiles( "which to sample the velocity profiles." ) + if x_inertial_start is None: + x_inertial_start = 0.0 + + if y_inertial_start is None: + y_inertial_start = 0.0 + if reference_height is None: reference_height = self.floris.flow_field.reference_wind_height - print( + self.reinitialize( + wind_directions=[wind_direction], + wind_speeds=[homogeneous_wind_speed], + wind_shear=0.0 + ) + +# print( +# direction, +# downstream_dists, +# profile_range, +# resolution, +# wind_direction, +# homogeneous_wind_speed, +# ref_rotor_diameter, +# x_inertial_start, +# y_inertial_start, +# reference_height, +# self.floris.flow_field.wind_directions, +# self.floris.flow_field.wind_speeds, +# self.floris.flow_field.wind_shear +# ) + + velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( direction, downstream_dists, profile_range, resolution, - wind_direction, - homogeneous_wind_speed, ref_rotor_diameter, x_inertial_start, y_inertial_start, reference_height ) + self.reinitialize( + wind_directions=wind_directions_copy, + wind_speeds=wind_speeds_copy, + wind_shear=wind_shear_copy + ) + + print(velocity_deficit_profiles) # Placeholder -# velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( -# direction, -# downstream_dists, -# profile_range, -# resolution, -# ref_turbine_diameter, -# x_inertial_start, -# y_inertial_start -# ) -# # return velocity_deficit_profiles @property From b74f0edcdd0d5a18d9773b81d1f0a9902bca7148 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 29 Jun 2023 17:26:46 +0200 Subject: [PATCH 13/46] Velocity profiles returned successfully --- examples/29_plot_velocity_deficit_profiles.py | 10 ++++- floris/simulation/floris.py | 38 ++++++++++++++++++- floris/tools/floris_interface.py | 11 +++--- 3 files changed, 51 insertions(+), 8 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index f89deabe4..ace10ed6e 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -27,7 +27,15 @@ fi = FlorisInterface("inputs/gch.yaml") fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) -fi.sample_velocity_deficit_profiles(direction = 'y', resolution=10, homogeneous_wind_speed=7.0) +velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( + direction = 'y', + resolution=10, + homogeneous_wind_speed=7.0 +) + +print(velocity_deficit_profiles[0]) +print(velocity_deficit_profiles[-1]) + #horizontal_plane = fi.calculate_horizontal_plane( # x_resolution=200, diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 381b50775..ef551bf0c 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -16,6 +16,8 @@ from pathlib import Path +import numpy as np +import pandas as pd import yaml from attrs import define, field @@ -321,13 +323,14 @@ def solve_for_velocity_deficit_profiles( downstream_dists, profile_range, resolution, + homogeneous_wind_speed, ref_rotor_diameter, x_inertial_start, y_inertial_start, reference_height ): - VelocityProfileGrid( + field_grid = VelocityProfileGrid( direction=direction, downstream_dists=downstream_dists, profile_range=profile_range, @@ -346,7 +349,40 @@ def solve_for_velocity_deficit_profiles( time_series=self.flow_field.time_series ) + self.flow_field.initialize_velocity_field(field_grid) + + vel_model = self.wake.model_strings["velocity_model"] + + if vel_model in ("cc", "turbopark"): + raise NotImplementedError( + "solve_for_velocity_deficit_profiles is currently only available with the " + "gauss, jensen, and empirical_guass models." + ) + elif vel_model == "empirical_gauss": + full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) + else: + full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) + + x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (-1, resolution)) + y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (-1, resolution)) + z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (-1, resolution)) + u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (-1, resolution)) + velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed + + velocity_deficit_profiles = [] + + for i in range(len(downstream_dists)): + df = pd.DataFrame( + { + "x/D": x[i]/ref_rotor_diameter, + "y/D": y[i]/ref_rotor_diameter, + "z/D": z[i]/ref_rotor_diameter, + "velocity_deficit": velocity_deficit[i] + } + ) + velocity_deficit_profiles.append(df) + return velocity_deficit_profiles def finalize(self): # Once the wake calculation is finished, unsort the values to match diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index fc4814ea0..307859d78 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -989,7 +989,7 @@ def sample_velocity_deficit_profiles( ref_rotor_diameter = rotor_diameters[0] else: raise ValueError( - "'ref_rotor_diameter' needs to be defined manually since there are multiple" + "'ref_rotor_diameter' needs to be defined manually since there are multiple " "turbine types." ) @@ -1008,7 +1008,7 @@ def sample_velocity_deficit_profiles( wind_direction = wind_directions_copy[0] else: raise ValueError( - "Multiple wind directions detected. Specify a single 'wind_direction' for" + "Multiple wind directions detected. Specify a single 'wind_direction' for " "which to sample the velocity profiles." ) @@ -1018,7 +1018,7 @@ def sample_velocity_deficit_profiles( # Maybe add info msg that wind_speed is homogeneous else: raise ValueError( - "Multiple wind speeds detected. Specify a single 'homogeneous_wind_speed' for" + "Multiple wind speeds detected. Specify a single 'homogeneous_wind_speed' for " "which to sample the velocity profiles." ) @@ -1058,6 +1058,7 @@ def sample_velocity_deficit_profiles( downstream_dists, profile_range, resolution, + homogeneous_wind_speed, ref_rotor_diameter, x_inertial_start, y_inertial_start, @@ -1070,9 +1071,7 @@ def sample_velocity_deficit_profiles( wind_shear=wind_shear_copy ) - print(velocity_deficit_profiles) # Placeholder - -# return velocity_deficit_profiles + return velocity_deficit_profiles @property def layout_x(self): From 51079f7947be4e68698a0668838f0206f0d4091c Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 30 Jun 2023 18:30:19 +0200 Subject: [PATCH 14/46] Small fixes --- examples/29_plot_velocity_deficit_profiles.py | 12 ++++++++-- floris/tools/floris_interface.py | 24 ++++--------------- 2 files changed, 14 insertions(+), 22 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index ace10ed6e..81ec2eda6 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -33,9 +33,17 @@ homogeneous_wind_speed=7.0 ) -print(velocity_deficit_profiles[0]) -print(velocity_deficit_profiles[-1]) +for df in velocity_deficit_profiles: + print(df) +velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( + direction = 'z', + resolution=10, + homogeneous_wind_speed=7.0 +) + +for df in velocity_deficit_profiles: + print(df) #horizontal_plane = fi.calculate_horizontal_plane( # x_resolution=200, diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 307859d78..409414af2 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -989,8 +989,8 @@ def sample_velocity_deficit_profiles( ref_rotor_diameter = rotor_diameters[0] else: raise ValueError( - "'ref_rotor_diameter' needs to be defined manually since there are multiple " - "turbine types." + f"Multiple rotor_diameters detected: {np.unique(rotor_diameters)}. " + "Provide a single 'ref_rotor_diameter'." ) if downstream_dists is None: @@ -1008,7 +1008,7 @@ def sample_velocity_deficit_profiles( wind_direction = wind_directions_copy[0] else: raise ValueError( - "Multiple wind directions detected. Specify a single 'wind_direction' for " + "Multiple wind directions detected. Provide a single 'wind_direction' for " "which to sample the velocity profiles." ) @@ -1018,7 +1018,7 @@ def sample_velocity_deficit_profiles( # Maybe add info msg that wind_speed is homogeneous else: raise ValueError( - "Multiple wind speeds detected. Specify a single 'homogeneous_wind_speed' for " + "Multiple wind speeds detected. Provide a single 'homogeneous_wind_speed' for " "which to sample the velocity profiles." ) @@ -1037,22 +1037,6 @@ def sample_velocity_deficit_profiles( wind_shear=0.0 ) -# print( -# direction, -# downstream_dists, -# profile_range, -# resolution, -# wind_direction, -# homogeneous_wind_speed, -# ref_rotor_diameter, -# x_inertial_start, -# y_inertial_start, -# reference_height, -# self.floris.flow_field.wind_directions, -# self.floris.flow_field.wind_speeds, -# self.floris.flow_field.wind_shear -# ) - velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( direction, downstream_dists, From 5c40659ab45f76091d538a4a19868ae39e44a1d1 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 2 Jul 2023 16:55:25 +0200 Subject: [PATCH 15/46] Clarifications --- floris/simulation/floris.py | 32 +++++++++++++++++--------------- floris/simulation/grid.py | 1 + floris/tools/floris_interface.py | 13 +++++++++---- 3 files changed, 27 insertions(+), 19 deletions(-) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index ef551bf0c..fb9950b24 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -44,6 +44,7 @@ VelocityProfileGrid, WakeModelManager, ) +from floris.type_dec import NDArrayFloat from floris.utilities import load_yaml @@ -319,16 +320,16 @@ def solve_for_points(self, x, y, z): def solve_for_velocity_deficit_profiles( self, - direction, - downstream_dists, - profile_range, - resolution, - homogeneous_wind_speed, - ref_rotor_diameter, - x_inertial_start, - y_inertial_start, - reference_height - ): + direction: str, + downstream_dists: NDArrayFloat | list, + profile_range: NDArrayFloat | list, + resolution: int, + homogeneous_wind_speed: float, + ref_rotor_diameter: float, + x_inertial_start: float, + y_inertial_start: float, + reference_height: float + ) -> list[pd.DataFrame]: field_grid = VelocityProfileGrid( direction=direction, @@ -363,15 +364,16 @@ def solve_for_velocity_deficit_profiles( else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) - x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (-1, resolution)) - y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (-1, resolution)) - z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (-1, resolution)) - u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (-1, resolution)) + nProfiles = len(downstream_dists) + x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (nProfiles, resolution)) + y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (nProfiles, resolution)) + z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (nProfiles, resolution)) + u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (nProfiles, resolution)) velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed velocity_deficit_profiles = [] - for i in range(len(downstream_dists)): + for i in range(nProfiles): df = pd.DataFrame( { "x/D": x[i]/ref_rotor_diameter, diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 3a93d8746..8fd4f3d83 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -771,6 +771,7 @@ def set_grid(self) -> None: """ res = self.resolution nProfiles = len(self.downstream_dists) + # downsteam_dists are defined from the following starting point coordinates_inertial_start = np.array( [[self.x_inertial_start, self.y_inertial_start, self.reference_height]] diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 409414af2..fb2bf7f07 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -952,17 +952,17 @@ def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloa def sample_velocity_deficit_profiles( self, - direction: NDArrayFloat = 'y', + direction: str = 'y', downstream_dists: NDArrayFloat | list = None, profile_range: NDArrayFloat | list = None, - resolution=100, + resolution: int = 100, wind_direction: float = None, homogeneous_wind_speed: float = None, ref_rotor_diameter: float = None, x_inertial_start: float = None, y_inertial_start: float = None, reference_height: float = None - ): + ) -> list[pd.DataFrame]: """ Extract velocity deficit profiles @@ -1015,7 +1015,12 @@ def sample_velocity_deficit_profiles( if homogeneous_wind_speed is None: if len(wind_speeds_copy) == 1: homogeneous_wind_speed = wind_speeds_copy[0] - # Maybe add info msg that wind_speed is homogeneous + self.logger.warning( + "'homogeneous_wind_speed' not provided. Setting it to the single wind speed " + "found in 'wind_speeds'. Note that the inflow is always homogeneous when " + "calculateing the velocity deficit profiles. This is done by temporarily " + "setting 'wind_shear' to 0.0" + ) else: raise ValueError( "Multiple wind speeds detected. Provide a single 'homogeneous_wind_speed' for " From 0a7ba7b54e178ee1683200c0507f3568031cbc32 Mon Sep 17 00:00:00 2001 From: vallbog Date: Tue, 4 Jul 2023 17:35:05 +0200 Subject: [PATCH 16/46] Draft of init for VelocityProfilesFigure class --- examples/29_plot_velocity_deficit_profiles.py | 37 +++++++++---- floris/simulation/floris.py | 12 ++--- floris/simulation/grid.py | 13 ++--- floris/tools/floris_interface.py | 8 +-- floris/tools/visualization.py | 54 +++++++++++++++++++ 5 files changed, 99 insertions(+), 25 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index 81ec2eda6..1d8bfaf87 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -18,32 +18,51 @@ import floris.tools.visualization as wakeviz from floris.tools import FlorisInterface +from floris.tools.visualization import VelocityProfilesFigure """ docstr """ - +D = 126.0 fi = FlorisInterface("inputs/gch.yaml") fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) +downstream_dists = D * np.array([3, 5, 7]) + velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( direction = 'y', + downstream_dists=downstream_dists, resolution=10, homogeneous_wind_speed=7.0 ) -for df in velocity_deficit_profiles: - print(df) +#for df in velocity_deficit_profiles: +# print(df) -velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( - direction = 'z', - resolution=10, - homogeneous_wind_speed=7.0 +profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y', 'z'] +) + +downstream_dists = D * np.array([1, 3, 5, 7, 9]) +profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['z'] ) -for df in velocity_deficit_profiles: - print(df) +plt.show() + + +#velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( +# direction = 'z', + +# resolution=10, +# homogeneous_wind_speed=7.0 +#) +# +#for df in velocity_deficit_profiles: +# print(df) #horizontal_plane = fi.calculate_horizontal_plane( # x_resolution=200, diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index fb9950b24..16240dd8e 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -364,16 +364,16 @@ def solve_for_velocity_deficit_profiles( else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) - nProfiles = len(downstream_dists) - x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (nProfiles, resolution)) - y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (nProfiles, resolution)) - z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (nProfiles, resolution)) - u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (nProfiles, resolution)) + nprofiles = len(downstream_dists) + x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (nprofiles, resolution)) + y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (nprofiles, resolution)) + z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (nprofiles, resolution)) + u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (nprofiles, resolution)) velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed velocity_deficit_profiles = [] - for i in range(nProfiles): + for i in range(nprofiles): df = pd.DataFrame( { "x/D": x[i]/ref_rotor_diameter, diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 8fd4f3d83..c2e5fdd5c 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -770,13 +770,14 @@ def set_grid(self) -> None: Set points for calculation based on a series of user-supplied coordinates. """ res = self.resolution - nProfiles = len(self.downstream_dists) + nprofiles = len(self.downstream_dists) # downsteam_dists are defined from the following starting point coordinates_inertial_start = np.array( [[self.x_inertial_start, self.y_inertial_start, self.reference_height]] ) + # Starting point in rotated coordinates x_start, y_start, _, _, _ = rotate_coordinates_rel_west( self.wind_directions, coordinates_inertial_start, @@ -786,7 +787,7 @@ def set_grid(self) -> None: x_start, y_start = x_start[0,0,0], y_start[0,0,0] downstream_dists_transpose = np.atleast_2d(self.downstream_dists).T - x = (x_start + downstream_dists_transpose) * np.ones((nProfiles, res)) + x = (x_start + downstream_dists_transpose) * np.ones((nprofiles, res)) if self.direction == 'y': y_single_profile = np.linspace( @@ -794,8 +795,8 @@ def set_grid(self) -> None: y_start + self.profile_range[1], res ) - y = y_single_profile * np.ones((nProfiles, res)) - z = self.reference_height * np.ones((nProfiles, res)) + y = y_single_profile * np.ones((nprofiles, res)) + z = self.reference_height * np.ones((nprofiles, res)) elif self.direction == 'z': z_min_profile = self.reference_height + self.profile_range[0] if z_min_profile <= 0.0: @@ -806,8 +807,8 @@ def set_grid(self) -> None: self.reference_height + self.profile_range[1], res ) - z = z_single_profile * np.ones((nProfiles, res)) - y = y_start * np.ones((nProfiles, res)) + z = z_single_profile * np.ones((nprofiles, res)) + y = y_start * np.ones((nprofiles, res)) self.x_sorted = x.flatten()[None,None,:,None,None] self.y_sorted = y.flatten()[None,None,:,None,None] diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index fb2bf7f07..759742a7a 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -984,12 +984,12 @@ def sample_velocity_deficit_profiles( raise ValueError("'direction' must be either y or z") if ref_rotor_diameter is None: - rotor_diameters = self.floris.farm.rotor_diameters - if np.all(rotor_diameters == rotor_diameters[0]): - ref_rotor_diameter = rotor_diameters[0] + unique_rotor_diameters = np.unique(self.floris.farm.rotor_diameters) + if len(unique_rotor_diameters) == 1: + ref_rotor_diameter = unique_rotor_diameters[0] else: raise ValueError( - f"Multiple rotor_diameters detected: {np.unique(rotor_diameters)}. " + f"Multiple rotor_diameters detected: {unique_rotor_diameters}. " "Provide a single 'ref_rotor_diameter'." ) diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index 18efebdc6..a6034e09e 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -17,17 +17,23 @@ import warnings from typing import Union +import attrs import matplotlib as mpl import matplotlib.colors as mplcolors import matplotlib.pyplot as plt import numpy as np import pandas as pd +from attrs import define, field from matplotlib import rcParams from scipy.spatial import ConvexHull from floris.simulation import Floris from floris.tools.cut_plane import CutPlane from floris.tools.floris_interface import FlorisInterface +from floris.type_dec import ( + floris_array_converter, + NDArrayFloat, +) from floris.utilities import rotate_coordinates_rel_west, wind_delta @@ -695,3 +701,51 @@ def calculate_horizontal_plane_with_turbines( horizontal_plane = CutPlane(df, x_resolution, y_resolution, "z") return horizontal_plane + +@define +class VelocityProfilesFigure(): + """ + docstr + """ + downstream_dists_D: NDArrayFloat = field(converter=floris_array_converter) + layout: list[str] = field(default=['y']) + + nrows: int = field(init=False) + ncols: int = field(init=False) + fig: plt.Figure = field(init=False) + axs: np.ndarray | plt.Axes = field(init=False) + + def __attrs_post_init__(self): + self.nrows = len(self.layout) + self.ncols = len(self.downstream_dists_D) + width_per_col = 6.4 / 3 + height_per_row = 6.4 / 2 + figsize = [width_per_col * self.ncols, height_per_row * self.nrows] + self.fig, axs = plt.subplots( + self.nrows, + self.ncols, + figsize=figsize, + layout='tight', + sharex='all', + sharey='row' + ) + + @layout.validator + def layout_validator(self, instance : attrs.Attribute, value : list[str]) -> None: + allowed_layouts = [['y'], ['z'], ['y', 'z'], ['z', 'y']] + if value not in allowed_layouts: + raise ValueError(f"'layout' must be one of the following: {allowed_layouts}") + + +#def find_velocity_profile_direction(df): +# unique_y = np.unique(df['y/D']) +# unique_z = np.unique(df['z/D']) +# if len(unique_y) == 1: +# direction = 'y' +# elif len(unique_z) == 1: +# direction = 'z' +# else: +# raise ValueError( +# "Velocity profile at" +# ) +# From 621556aa20bbac5728dd6eb251636bd601fb56a7 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 6 Jul 2023 11:35:57 +0200 Subject: [PATCH 17/46] Small addition to class init --- examples/29_plot_velocity_deficit_profiles.py | 7 +++++++ floris/tools/visualization.py | 15 ++++++++++++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index 1d8bfaf87..f7c746880 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -51,6 +51,13 @@ layout=['z'] ) + +downstream_dists = D * np.array([3]) +profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y'] +) + plt.show() diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index a6034e09e..de85c0c9c 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -719,16 +719,25 @@ def __attrs_post_init__(self): self.nrows = len(self.layout) self.ncols = len(self.downstream_dists_D) width_per_col = 6.4 / 3 - height_per_row = 6.4 / 2 - figsize = [width_per_col * self.ncols, height_per_row * self.nrows] + height_per_row = 7.0 / 2 + figsize = [0.5 + width_per_col * self.ncols, height_per_row * self.nrows] self.fig, axs = plt.subplots( self.nrows, self.ncols, figsize=figsize, layout='tight', - sharex='all', + sharex='col', sharey='row' ) + self.axs = np.atleast_2d(axs) + + for ax in self.axs[-1]: + ax.set_xlabel(r'$\Delta U / U_\infty$', fontsize=14) + ax.tick_params('x', labelsize=14) + + for profile_direction, ax in zip(self.layout, self.axs[:,0]): + ax.set_ylabel(f'${profile_direction}/D$', fontsize=14) + ax.tick_params('y', labelsize=14) @layout.validator def layout_validator(self, instance : attrs.Attribute, value : list[str]) -> None: From ae29d21ee8d80ec1526f9adc6151f0fa5d599d8a Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 21 Jul 2023 00:56:31 +0200 Subject: [PATCH 18/46] Methods for adding profiles and reference lines --- examples/29_plot_velocity_deficit_profiles.py | 60 ++++-------- floris/tools/floris_interface.py | 2 +- floris/tools/visualization.py | 91 ++++++++++++++++--- 3 files changed, 100 insertions(+), 53 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index f7c746880..c1eb9bf49 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -29,54 +29,34 @@ fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) downstream_dists = D * np.array([3, 5, 7]) +homogeneous_wind_speed = 8.0 -velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( - direction = 'y', +velocity_deficit_profiles_y = fi.sample_velocity_deficit_profiles( + direction='y', downstream_dists=downstream_dists, - resolution=10, - homogeneous_wind_speed=7.0 + homogeneous_wind_speed=homogeneous_wind_speed +) +# Same range as y-profiles for an easier comparison +profile_range_z = D * np.array([0, 4]) - 90.0 +velocity_deficit_profiles_z = fi.sample_velocity_deficit_profiles( + direction='z', + downstream_dists=downstream_dists, + profile_range=profile_range_z, + homogeneous_wind_speed=homogeneous_wind_speed ) - -#for df in velocity_deficit_profiles: -# print(df) profiles_fig = VelocityProfilesFigure( downstream_dists_D=downstream_dists / D, layout=['y', 'z'] ) - -downstream_dists = D * np.array([1, 3, 5, 7, 9]) -profiles_fig = VelocityProfilesFigure( - downstream_dists_D=downstream_dists / D, - layout=['z'] -) - - -downstream_dists = D * np.array([3]) -profiles_fig = VelocityProfilesFigure( - downstream_dists_D=downstream_dists / D, - layout=['y'] +profiles_fig.add_profiles( + velocity_deficit_profiles_y + velocity_deficit_profiles_z, + color='k' ) +margin = 0.05 +profiles_fig.set_xlim([0.0 - margin, 0.6 + margin]) +profiles_fig.add_ref_lines_y_D([-0.5, 0.5]) +ref_lines_z_D = 90.0 / D + np.array([-0.5, 0.5]) +profiles_fig.add_ref_lines_z_D(ref_lines_z_D) plt.show() - - -#velocity_deficit_profiles = fi.sample_velocity_deficit_profiles( -# direction = 'z', - -# resolution=10, -# homogeneous_wind_speed=7.0 -#) -# -#for df in velocity_deficit_profiles: -# print(df) - -#horizontal_plane = fi.calculate_horizontal_plane( -# x_resolution=200, -# y_resolution=100, -# height=90.0, -#) -# -#wakeviz.visualize_cut_plane(horizontal_plane, title="Horizontal plane") -# -#wakeviz.show_plots() diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 759742a7a..619aacea8 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -1018,7 +1018,7 @@ def sample_velocity_deficit_profiles( self.logger.warning( "'homogeneous_wind_speed' not provided. Setting it to the single wind speed " "found in 'wind_speeds'. Note that the inflow is always homogeneous when " - "calculateing the velocity deficit profiles. This is done by temporarily " + "calculating the velocity deficit profiles. This is done by temporarily " "setting 'wind_shear' to 0.0" ) else: diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index de85c0c9c..354364c00 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -735,6 +735,9 @@ def __attrs_post_init__(self): ax.set_xlabel(r'$\Delta U / U_\infty$', fontsize=14) ax.tick_params('x', labelsize=14) + for ax, x_D in zip(self.axs[0], self.downstream_dists_D): + ax.set_title(f'$x/D = {x_D:.1f}$', fontsize=14) + for profile_direction, ax in zip(self.layout, self.axs[:,0]): ax.set_ylabel(f'${profile_direction}/D$', fontsize=14) ax.tick_params('y', labelsize=14) @@ -745,16 +748,80 @@ def layout_validator(self, instance : attrs.Attribute, value : list[str]) -> Non if value not in allowed_layouts: raise ValueError(f"'layout' must be one of the following: {allowed_layouts}") + def add_profiles(self, velocity_deficit_profiles, **kwargs): + for df in velocity_deficit_profiles: + ax, profile_direction = self.match_profile_to_axes(df) + profile_direction_D = f'{profile_direction}/D' + ax.plot(df['velocity_deficit'], df[profile_direction_D], **kwargs) + + def match_profile_to_axes(self, df): + x_D = np.unique(df['x/D']) + if len(x_D) > 1: + raise ValueError( + "The streamwise location x/D must be constant for each velocity profile." + ) + + unique_y = np.unique(df['y/D']) + unique_z = np.unique(df['z/D']) + if len(unique_y) == 1: + profile_direction = 'z' + elif len(unique_z) == 1: + profile_direction = 'y' + else: + raise ValueError( + "Velocity deficit profile at x/D = {x_D} is neither in the cross-stream (y)" + "nor the vertical (z) direction." + ) + row = self.layout.index(profile_direction) + + col = None + for i in range(self.ncols): + if np.abs(x_D - self.downstream_dists_D[i]) < 0.001: + col = i + break + if col is None: + raise ValueError( + "Could not add a velocity deficit profile at downstream distance " + f"x/D = {x_D}. The downstream distance must be one of the following " + "values with which the instance of the VelocityProfilesFigure was initialized: " + f"{self.downstream_dists_D}" + ) + return self.axs[row][col], profile_direction -#def find_velocity_profile_direction(df): -# unique_y = np.unique(df['y/D']) -# unique_z = np.unique(df['z/D']) -# if len(unique_y) == 1: -# direction = 'y' -# elif len(unique_z) == 1: -# direction = 'z' -# else: -# raise ValueError( -# "Velocity profile at" -# ) -# + def set_xlim(self, xlim): + for ax in self.axs[-1]: + ax.set_xlim(xlim) + + def add_ref_lines_y_D(self, ref_lines, **kwargs): + try: + row_y = self.layout.index('y') + except Exception: + print( + "Could not add reference lines to cross-stream (y) velocity profiles. No " + "such profiles exist in the figure." + ) + self.add_ref_lines(ref_lines, row_y, **kwargs) + + def add_ref_lines_z_D(self, ref_lines, **kwargs): + try: + row_z = self.layout.index('z') + except Exception: + print( + "Could not add reference lines to vertical (z) velocity profiles. No " + "such profiles exist in the figure." + ) + self.add_ref_lines(ref_lines, row_z, **kwargs) + + def add_ref_lines(self, ref_lines, row, **kwargs): + default_params = { + 'linestyle': (0, (4, 2)), + 'color': 'k', + 'linewidth': 1.2 + } + for key in default_params: + if key not in kwargs: + kwargs[key] = default_params[key] + + for ax in self.axs[row]: + for coordinate in ref_lines: + ax.plot([0.0, 1.0], [coordinate, coordinate], **kwargs) From 3b847a168ffea8b7ffd5343645ccb9040aae9650 Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 21 Jul 2023 18:13:22 +0200 Subject: [PATCH 19/46] Improved example --- examples/29_plot_velocity_deficit_profiles.py | 111 ++++++++++++------ floris/tools/visualization.py | 26 ++-- 2 files changed, 88 insertions(+), 49 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index c1eb9bf49..21cb433e4 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -24,39 +24,78 @@ """ docstr """ -D = 126.0 -fi = FlorisInterface("inputs/gch.yaml") -fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) - -downstream_dists = D * np.array([3, 5, 7]) -homogeneous_wind_speed = 8.0 - -velocity_deficit_profiles_y = fi.sample_velocity_deficit_profiles( - direction='y', - downstream_dists=downstream_dists, - homogeneous_wind_speed=homogeneous_wind_speed -) -# Same range as y-profiles for an easier comparison -profile_range_z = D * np.array([0, 4]) - 90.0 -velocity_deficit_profiles_z = fi.sample_velocity_deficit_profiles( - direction='z', - downstream_dists=downstream_dists, - profile_range=profile_range_z, - homogeneous_wind_speed=homogeneous_wind_speed -) - -profiles_fig = VelocityProfilesFigure( - downstream_dists_D=downstream_dists / D, - layout=['y', 'z'] -) -profiles_fig.add_profiles( - velocity_deficit_profiles_y + velocity_deficit_profiles_z, - color='k' -) -margin = 0.05 -profiles_fig.set_xlim([0.0 - margin, 0.6 + margin]) -profiles_fig.add_ref_lines_y_D([-0.5, 0.5]) -ref_lines_z_D = 90.0 / D + np.array([-0.5, 0.5]) -profiles_fig.add_ref_lines_z_D(ref_lines_z_D) - -plt.show() + +def get_profiles(direction, profile_range=None, resolution=100): + return fi.sample_velocity_deficit_profiles( + direction=direction, + downstream_dists=downstream_dists, + profile_range=profile_range, + resolution=resolution, + homogeneous_wind_speed=homogeneous_wind_speed + ) + +if __name__ == '__main__': + D = 126.0 + hub_height = 90.0 + fi = FlorisInterface("inputs/gch.yaml") + fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) + + downstream_dists = D * np.array([3, 5, 7]) + homogeneous_wind_speed = 8.0 + + # Same range as y-profiles for an easier comparison + profile_range_z = D * np.array([0, 4]) - hub_height + profiles = get_profiles('y') + get_profiles('z', profile_range_z) + + profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y', 'z'], + ) + profiles_fig.add_profiles(profiles, color='k') + + # Change velocity model to jensen, get the velocity deficit profiles, + # and plot them + floris_dict = fi.floris.as_dict() + floris_dict['wake']['model_strings']['velocity_model'] = 'jensen' + fi = FlorisInterface(floris_dict) + profiles_y = get_profiles('y', resolution=400) + profiles_z = get_profiles('z', profile_range_z, resolution=400) + profiles_fig.add_profiles(profiles_y + profiles_z, color='r') + + margin = 0.05 + profiles_fig.set_xlim([0.0 - margin, 0.6 + margin]) + profiles_fig.add_ref_lines_y_D([-0.5, 0.5]) + ref_lines_z_D = hub_height / D + np.array([-0.5, 0.5]) + profiles_fig.add_ref_lines_z_D(ref_lines_z_D) + + profiles_fig.axs[0,0].legend(['gauss', 'jensen'], fontsize=11) + profiles_fig.fig.suptitle( + 'Velocity decifit profiles from different velocity models', + fontsize=14 + ) + + # Profiles are independent of wind direction for a single turbine. + # This is because x is always in the streamwise direction + fi = FlorisInterface("inputs/gch.yaml") + fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) + downstream_dists = D * np.array([3]) + profiles = get_profiles('y', resolution=400) + profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y'] + ) + profiles_fig.add_profiles(profiles, color='k') + + fi.reinitialize(wind_directions=[315.0]) + profiles = get_profiles('y', resolution=21) + profiles_fig.add_profiles( + profiles, + linestyle='None', + marker='o', + color='b', + markerfacecolor='None', + markeredgewidth=1.3 + ) + profiles_fig.axs[0,0].legend(['WD = 270 deg', 'WD = 315 deg']) + + plt.show() diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index 354364c00..bed335577 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -709,6 +709,8 @@ class VelocityProfilesFigure(): """ downstream_dists_D: NDArrayFloat = field(converter=floris_array_converter) layout: list[str] = field(default=['y']) + width_per_col: float = field(default=2.07) + height_per_row: float = field(default=3.0) nrows: int = field(init=False) ncols: int = field(init=False) @@ -718,9 +720,7 @@ class VelocityProfilesFigure(): def __attrs_post_init__(self): self.nrows = len(self.layout) self.ncols = len(self.downstream_dists_D) - width_per_col = 6.4 / 3 - height_per_row = 7.0 / 2 - figsize = [0.5 + width_per_col * self.ncols, height_per_row * self.nrows] + figsize = [0.7 + self.width_per_col * self.ncols, 1.0 + self.height_per_row * self.nrows] self.fig, axs = plt.subplots( self.nrows, self.ncols, @@ -756,7 +756,9 @@ def add_profiles(self, velocity_deficit_profiles, **kwargs): def match_profile_to_axes(self, df): x_D = np.unique(df['x/D']) - if len(x_D) > 1: + if len(x_D) == 1: + x_D = x_D[0] + else: raise ValueError( "The streamwise location x/D must be constant for each velocity profile." ) @@ -769,7 +771,7 @@ def match_profile_to_axes(self, df): profile_direction = 'y' else: raise ValueError( - "Velocity deficit profile at x/D = {x_D} is neither in the cross-stream (y)" + f"Velocity deficit profile at x/D = {x_D} is neither in the cross-stream (y) " "nor the vertical (z) direction." ) row = self.layout.index(profile_direction) @@ -793,23 +795,21 @@ def set_xlim(self, xlim): ax.set_xlim(xlim) def add_ref_lines_y_D(self, ref_lines, **kwargs): - try: - row_y = self.layout.index('y') - except Exception: - print( + if 'y' not in self.layout: + raise Exception( "Could not add reference lines to cross-stream (y) velocity profiles. No " "such profiles exist in the figure." ) + row_y = self.layout.index('y') self.add_ref_lines(ref_lines, row_y, **kwargs) def add_ref_lines_z_D(self, ref_lines, **kwargs): - try: - row_z = self.layout.index('z') - except Exception: - print( + if 'z' not in self.layout: + raise Exception( "Could not add reference lines to vertical (z) velocity profiles. No " "such profiles exist in the figure." ) + row_z = self.layout.index('z') self.add_ref_lines(ref_lines, row_z, **kwargs) def add_ref_lines(self, ref_lines, row, **kwargs): From a31cb878f2cb7e002ebe3435234ab7db4fa4c3a3 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sat, 22 Jul 2023 13:19:10 +0200 Subject: [PATCH 20/46] Always construct vawt_blade_lengts. Set to 0.0 for HAWTs --- floris/simulation/farm.py | 6 ++++-- floris/simulation/solver.py | 3 +++ floris/simulation/turbine.py | 2 +- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index a7b5187d5..956056de7 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -211,11 +211,13 @@ def construct_rotor_diameters(self): ]) def construct_vawt_blade_lengths(self): + vawt_blade_lengths = [] for turb in self.turbine_definitions: try: - self.vawt_blade_lengths = turb['vawt_blade_length'] + vawt_blade_lengths.append(turb['vawt_blade_length']) except KeyError: - self.vawt_blade_lengths = 0.0 + vawt_blade_lengths.append(0.0) + self.vawt_blade_lengths = np.array(vawt_blade_lengths) def construct_turbine_TSRs(self): self.TSRs = np.array([turb['TSR'] for turb in self.turbine_definitions]) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index b2a53c0ec..08b87987a 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -672,6 +672,7 @@ def full_flow_cc_solver( turbine_grid_farm.construct_turbine_power_interps() turbine_grid_farm.construct_hub_heights() turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_vawt_blade_lengths() turbine_grid_farm.construct_turbine_TSRs() turbine_grid_farm.construct_turbine_pPs() turbine_grid_farm.construct_turbine_pTs() @@ -1108,6 +1109,7 @@ def full_flow_turbopark_solver( # turbine_grid_farm.construct_turbine_power_interps() # turbine_grid_farm.construct_hub_heights() # turbine_grid_farm.construct_rotor_diameters() + # turbine_grid_farm.construct_vawt_blade_lengths() # turbine_grid_farm.construct_turbine_TSRs() # turbine_grid_farm.construc_turbine_pPs() # turbine_grid_farm.construct_coordinates() @@ -1347,6 +1349,7 @@ def full_flow_empirical_gauss_solver( turbine_grid_farm.construct_turbine_power_interps() turbine_grid_farm.construct_hub_heights() turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_vawt_blade_lengths() turbine_grid_farm.construct_turbine_TSRs() turbine_grid_farm.construct_turbine_pPs() turbine_grid_farm.construct_turbine_pTs() diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index 604cb6ef8..e7a312917 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -627,7 +627,7 @@ class Turbine(BaseClass): power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) floating_tilt_table = field(default=None) floating_correct_cp_ct_for_tilt = field(default=None) - vawt_blade_length = field(default=None) + vawt_blade_length: float = field(default=0.0) # rloc: float = float_attrib() # TODO: goes here or on the Grid? # use_points_on_perimeter: bool = bool_attrib() From e68007d18fd652ed924e35e10feea1d9b3e7da72 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 23 Jul 2023 14:51:13 +0200 Subject: [PATCH 21/46] Able to use vawt solver for velocity profiles --- floris/simulation/floris.py | 4 +++- floris/simulation/wake_velocity/super_gaussian_vawt.py | 4 ++-- floris/tools/visualization.py | 2 +- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index d86e0ae17..3889671b8 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -371,10 +371,12 @@ def solve_for_velocity_deficit_profiles( if vel_model in ("cc", "turbopark"): raise NotImplementedError( "solve_for_velocity_deficit_profiles is currently only available with the " - "gauss, jensen, and empirical_guass models." + "gauss, jensen, empirical_guass, and super_gaussian_vawt models." ) elif vel_model == "empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) + elif vel_model == "super_gaussian_vawt": + full_flow_vawt_solver(self.farm, self.flow_field, field_grid, self.wake) else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index f879e20aa..6da4188e7 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -60,8 +60,8 @@ class SuperGaussianVAWTVelocityDeficit(BaseModel): :style: unsrt :filter: docname in docnames """ - wake_expansion_coeff_y: float = field(default=0.50) - wake_expansion_coeff_z: float = field(default=0.50) + wake_expansion_coeff_y: float = field(default=0.45) + wake_expansion_coeff_z: float = field(default=0.45) ay: float = field(default=0.95) az: float = field(default=4.5) by: float = field(default=0.35) diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index bed335577..f348b27a4 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -816,7 +816,7 @@ def add_ref_lines(self, ref_lines, row, **kwargs): default_params = { 'linestyle': (0, (4, 2)), 'color': 'k', - 'linewidth': 1.2 + 'linewidth': 1.1 } for key in default_params: if key not in kwargs: From 158cb995a437a35704fdbddd385ade0a35652fdd Mon Sep 17 00:00:00 2001 From: vallbog Date: Wed, 26 Jul 2023 08:43:40 +0200 Subject: [PATCH 22/46] Corrected self.axs --- floris/tools/visualization.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index bed335577..3ae663dbd 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -721,15 +721,15 @@ def __attrs_post_init__(self): self.nrows = len(self.layout) self.ncols = len(self.downstream_dists_D) figsize = [0.7 + self.width_per_col * self.ncols, 1.0 + self.height_per_row * self.nrows] - self.fig, axs = plt.subplots( + self.fig, self.axs = plt.subplots( self.nrows, self.ncols, figsize=figsize, layout='tight', sharex='col', - sharey='row' + sharey='row', + squeeze=False ) - self.axs = np.atleast_2d(axs) for ax in self.axs[-1]: ax.set_xlabel(r'$\Delta U / U_\infty$', fontsize=14) @@ -788,7 +788,7 @@ def match_profile_to_axes(self, df): "values with which the instance of the VelocityProfilesFigure was initialized: " f"{self.downstream_dists_D}" ) - return self.axs[row][col], profile_direction + return self.axs[row,col], profile_direction def set_xlim(self, xlim): for ax in self.axs[-1]: @@ -816,7 +816,7 @@ def add_ref_lines(self, ref_lines, row, **kwargs): default_params = { 'linestyle': (0, (4, 2)), 'color': 'k', - 'linewidth': 1.2 + 'linewidth': 1.1 } for key in default_params: if key not in kwargs: From bb62a153e62d605094195abc20494d75d6d513ce Mon Sep 17 00:00:00 2001 From: vallbog Date: Tue, 1 Aug 2023 10:23:43 +0200 Subject: [PATCH 23/46] Small detail --- floris/simulation/solver.py | 3 ++- floris/simulation/wake_velocity/super_gaussian_vawt.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 08b87987a..c2ac6818e 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -1513,7 +1513,8 @@ def vawt_solver( # Integrate this into the main data structure. # Move on to the next turbine. # - # Note that no turbulence model has been implmented in this solver. + # Note that no turbulence model has been implemented in this solver. This means that all + # turbines percieve a fixed turbulence intensity of flow_field.turbulence_intensity # <> deficit_model_args = model_manager.velocity_model.prepare_function(grid, flow_field) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 6da4188e7..0f798e0e1 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -171,7 +171,7 @@ def function( # At any streamwise location x behind the turbine, the velocity deficit in the # y-z plane is given by multiplying the maximum deficit C = C(x) with two Gaussian # shape functions fy = fy(y) and fz = fz(z) - eta = 1 / ny + 1 / nz + eta = ne.evaluate("1 / ny + 1 / nz") ny_inv = ne.evaluate("1 / ny") nz_inv = ne.evaluate("1 / nz") gamma_y, gamma_z = gamma(ny_inv), gamma(nz_inv) From 77d0e42f33f0baf75de1ecea2d3deec178b4c83d Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 3 Aug 2023 10:12:00 +0200 Subject: [PATCH 24/46] Comments, style check, and updated output coords --- examples/29_plot_velocity_deficit_profiles.py | 41 ++++--- floris/simulation/__init__.py | 2 +- floris/simulation/floris.py | 50 +++++--- floris/simulation/flow_field.py | 14 ++- floris/simulation/grid.py | 96 +++++++++------ floris/tools/floris_interface.py | 82 +++++++------ floris/tools/visualization.py | 110 +++++++++++++----- 7 files changed, 259 insertions(+), 136 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index 21cb433e4..cb347d6f9 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -25,13 +25,12 @@ docstr """ -def get_profiles(direction, profile_range=None, resolution=100): +def get_profiles(direction, resolution=100): return fi.sample_velocity_deficit_profiles( direction=direction, downstream_dists=downstream_dists, - profile_range=profile_range, resolution=resolution, - homogeneous_wind_speed=homogeneous_wind_speed + homogeneous_wind_speed=homogeneous_wind_speed, ) if __name__ == '__main__': @@ -43,9 +42,11 @@ def get_profiles(direction, profile_range=None, resolution=100): downstream_dists = D * np.array([3, 5, 7]) homogeneous_wind_speed = 8.0 - # Same range as y-profiles for an easier comparison - profile_range_z = D * np.array([0, 4]) - hub_height - profiles = get_profiles('y') + get_profiles('z', profile_range_z) + # Velocity deficit profiles can be obtained in either the cross-stream (y) + # or the vertical direction (z). The default profile_range is [-2 * D, 2 * D] with + # the default origin being (x, y, z) = (0.0, 0.0, reference_wind_height). + # TODO: Explain difference between starting point and coordinates seen in fig + profiles = get_profiles('y') + get_profiles('z') profiles_fig = VelocityProfilesFigure( downstream_dists_D=downstream_dists / D, @@ -54,19 +55,18 @@ def get_profiles(direction, profile_range=None, resolution=100): profiles_fig.add_profiles(profiles, color='k') # Change velocity model to jensen, get the velocity deficit profiles, - # and plot them + # and add them to the figure floris_dict = fi.floris.as_dict() floris_dict['wake']['model_strings']['velocity_model'] = 'jensen' fi = FlorisInterface(floris_dict) profiles_y = get_profiles('y', resolution=400) - profiles_z = get_profiles('z', profile_range_z, resolution=400) + profiles_z = get_profiles('z', resolution=400) profiles_fig.add_profiles(profiles_y + profiles_z, color='r') margin = 0.05 profiles_fig.set_xlim([0.0 - margin, 0.6 + margin]) - profiles_fig.add_ref_lines_y_D([-0.5, 0.5]) - ref_lines_z_D = hub_height / D + np.array([-0.5, 0.5]) - profiles_fig.add_ref_lines_z_D(ref_lines_z_D) + profiles_fig.add_ref_lines_y([-0.5, 0.5]) + profiles_fig.add_ref_lines_z([-0.5, 0.5]) profiles_fig.axs[0,0].legend(['gauss', 'jensen'], fontsize=11) profiles_fig.fig.suptitle( @@ -74,8 +74,10 @@ def get_profiles(direction, profile_range=None, resolution=100): fontsize=14 ) - # Profiles are independent of wind direction for a single turbine. - # This is because x is always in the streamwise direction + # Show that profiles are independent of the wind direction for a single turbine. + # This is because x is always in the streamwise direction. + # Also show that the coordinates x / D, y / D and z / D returned by + # sample_velocity_deficit_profiles are relative to the starting point mentioned above. fi = FlorisInterface("inputs/gch.yaml") fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) downstream_dists = D * np.array([3]) @@ -86,8 +88,15 @@ def get_profiles(direction, profile_range=None, resolution=100): ) profiles_fig.add_profiles(profiles, color='k') - fi.reinitialize(wind_directions=[315.0]) - profiles = get_profiles('y', resolution=21) + fi.reinitialize(wind_directions=[315.0], layout_x=[D], layout_y=[D]) + profiles = fi.sample_velocity_deficit_profiles( + direction='y', + downstream_dists=downstream_dists, + resolution=21, + homogeneous_wind_speed=homogeneous_wind_speed, + x_inertial_start=D, + y_inertial_start=D, + ) profiles_fig.add_profiles( profiles, linestyle='None', @@ -96,6 +105,6 @@ def get_profiles(direction, profile_range=None, resolution=100): markerfacecolor='None', markeredgewidth=1.3 ) - profiles_fig.axs[0,0].legend(['WD = 270 deg', 'WD = 315 deg']) + profiles_fig.axs[0,0].legend(['WD = 270 deg', 'WD = 315 deg,\nand turbine\nmoved']) plt.show() diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index c87f39d6e..3cf4e6887 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -44,7 +44,7 @@ FlowFieldPlanarGrid, Grid, PointsGrid, - VelocityProfileGrid, + LinesGrid, TurbineGrid, TurbineCubatureGrid ) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 16240dd8e..964b7112d 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -35,13 +35,13 @@ full_flow_sequential_solver, full_flow_turbopark_solver, Grid, + LinesGrid, PointsGrid, sequential_solver, State, TurbineCubatureGrid, TurbineGrid, turbopark_solver, - VelocityProfileGrid, WakeModelManager, ) from floris.type_dec import NDArrayFloat @@ -328,15 +328,20 @@ def solve_for_velocity_deficit_profiles( ref_rotor_diameter: float, x_inertial_start: float, y_inertial_start: float, - reference_height: float - ) -> list[pd.DataFrame]: + reference_height: float, + ) -> list[pd.DataFrame]: + """ + Extract velocity deficit profiles. See + :py:meth:`~floris.tools.floris_interface.FlorisInterface.sample_velocity_deficit_profiles` + for more details. + """ - field_grid = VelocityProfileGrid( + # Create a grid that contains coordinates for all the sample points in all profiles + field_grid = LinesGrid( direction=direction, downstream_dists=downstream_dists, - profile_range=profile_range, + line_range=profile_range, resolution=resolution, - ref_rotor_diameter=ref_rotor_diameter, x_inertial_start=x_inertial_start, y_inertial_start=y_inertial_start, reference_height=reference_height, @@ -347,14 +352,14 @@ def solve_for_velocity_deficit_profiles( grid_resolution=1, wind_directions=self.flow_field.wind_directions, wind_speeds=self.flow_field.wind_speeds, - time_series=self.flow_field.time_series + time_series=self.flow_field.time_series, ) self.flow_field.initialize_velocity_field(field_grid) vel_model = self.wake.model_strings["velocity_model"] - if vel_model in ("cc", "turbopark"): + if vel_model in ["cc", "turbopark"]: raise NotImplementedError( "solve_for_velocity_deficit_profiles is currently only available with the " "gauss, jensen, and empirical_guass models." @@ -364,22 +369,31 @@ def solve_for_velocity_deficit_profiles( else: full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) - nprofiles = len(downstream_dists) - x = np.reshape(field_grid.x_sorted[0,0,:,0,0], (nprofiles, resolution)) - y = np.reshape(field_grid.y_sorted[0,0,:,0,0], (nprofiles, resolution)) - z = np.reshape(field_grid.z_sorted[0,0,:,0,0], (nprofiles, resolution)) - u = np.reshape(self.flow_field.u_sorted[0,0,:,0,0], (nprofiles, resolution)) + n_profiles = len(downstream_dists) + x_relative_start = np.reshape( + field_grid.x_sorted[0, 0, :, 0, 0] - field_grid.x_start, + (n_profiles, resolution), + ) + y_relative_start = np.reshape( + field_grid.y_sorted[0, 0, :, 0, 0] - field_grid.y_start, + (n_profiles, resolution), + ) + z_relative_start = np.reshape( + field_grid.z_sorted[0, 0, :, 0, 0] - reference_height, + (n_profiles, resolution), + ) + u = np.reshape(self.flow_field.u_sorted[0, 0, :, 0, 0], (n_profiles, resolution)) velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed velocity_deficit_profiles = [] - for i in range(nprofiles): + for i in range(n_profiles): df = pd.DataFrame( { - "x/D": x[i]/ref_rotor_diameter, - "y/D": y[i]/ref_rotor_diameter, - "z/D": z[i]/ref_rotor_diameter, - "velocity_deficit": velocity_deficit[i] + "x/D": x_relative_start[i]/ref_rotor_diameter, + "y/D": y_relative_start[i]/ref_rotor_diameter, + "z/D": z_relative_start[i]/ref_rotor_diameter, + "velocity_deficit": velocity_deficit[i], } ) velocity_deficit_profiles.append(df) diff --git a/floris/simulation/flow_field.py b/floris/simulation/flow_field.py index a8dbf7393..a287f065f 100644 --- a/floris/simulation/flow_field.py +++ b/floris/simulation/flow_field.py @@ -127,11 +127,15 @@ def initialize_velocity_field(self, grid: Grid) -> None: # for height, using it here to apply the shear law makes that dimension store the vertical # wind profile. wind_profile_plane = (grid.z_sorted / self.reference_wind_height) ** self.wind_shear - dwind_profile_plane = ( - self.wind_shear - * (1 / self.reference_wind_height) ** self.wind_shear - * (grid.z_sorted) ** (self.wind_shear - 1) - ) + if self.wind_shear == 0.0: + # Avoid division by zero at z = 0.0 when wind_shear = 0.0 + dwind_profile_plane = np.zeros_like(grid.z_sorted) + else: + dwind_profile_plane = ( + self.wind_shear + * (1 / self.reference_wind_height) ** self.wind_shear + * (grid.z_sorted) ** (self.wind_shear - 1) + ) # If no hetergeneous inflow defined, then set all speeds ups to 1.0 if self.het_map is None: diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index c2e5fdd5c..4c9bd6954 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -118,7 +118,7 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter # TODO move this to the grid types and off of the base class """Check that grid resolution is given as int or Vec3 with int components.""" if isinstance(value, int) and \ - isinstance(self, (TurbineGrid, TurbineCubatureGrid, PointsGrid, VelocityProfileGrid)): + isinstance(self, (TurbineGrid, TurbineCubatureGrid, PointsGrid, LinesGrid)): return elif isinstance(value, Iterable) and isinstance(self, FlowFieldPlanarGrid): assert type(value[0]) is int @@ -746,31 +746,62 @@ def set_grid(self) -> None: self.z_sorted = z[:,:,:,None,None] @define -class VelocityProfileGrid(Grid): +class LinesGrid(Grid): """ - docstr + Create a grid of parallel lines. For each specified downstream distance of a + starting point, a line is created in either the cross-stream (y) or the vertical + direction (z). + + Args: + direction: Direction of the lines. Either `y` or `z`. + downstream_dists: A list/array of streamwise locations for where to create the lines. + line_range: Determines the extent of the lines. The range is defined about a point + which lies some distance directly downstream of the starting point. `line_range` + is a two-element list/array. For example, if `direction` is y and `line_range` is + [-2 * D, 2 * D], then for every downstream distance, a line will be created + at some reference height for y_start - 2 * D <= y <= y_start + 2 * D, where D is + the turbine diameter and `y_start` is the y-coordinate of the starting point in + rotated coordinates. + resolution: Number of points in each line. + x_inertial_start: x-coordinate of starting point in the inertial coordinate system. + y_inertial_start: y-coordinate of starting point in the inertial coordinate system. + reference_height: If `direction` is y, then `reference_height` defines the height of the + xy-plane in which the lines are created. If `direction` is z, then each line is + created in the vertical direction with the `line_range` being relative to the + `reference_height`. + x_center_of_rotation: The x-coordinate of the center of rotation around which the + starting point is rotated to account for wind direction changes. If not supplied, + the center of rotation will be the same as the starting point. + y_center_of_rotation: The y-coordinate of the center of rotation around which the + starting point is rotated to account for wind direction changes. If not supplied, + the center of rotation will be the same as the starting point. + See :py:class:`floris.simulation.grid.Grid` for super arguments. """ direction: str downstream_dists: NDArrayFloat = field(converter=floris_array_converter) - profile_range: NDArrayFloat = field(converter=floris_array_converter) + line_range: NDArrayFloat = field(converter=floris_array_converter) resolution: int - ref_rotor_diameter: float x_inertial_start: float y_inertial_start: float reference_height: float - x_center_of_rotation: float | None = field(default=None) - y_center_of_rotation: float | None = field(default=None) + x_center_of_rotation: float = field(default=None) + y_center_of_rotation: float = field(default=None) + + x_start: float = field(init=False) + y_start: float = field(init=False) def __attrs_post_init__(self) -> None: super().__attrs_post_init__() + if len(self.wind_directions) > 1: + raise NotImplementedError( + "LinesGrid currently only works for a single wind direction." + ) + self.set_grid() def set_grid(self) -> None: - """ - Set points for calculation based on a series of user-supplied coordinates. - """ res = self.resolution - nprofiles = len(self.downstream_dists) + n_lines = len(self.downstream_dists) # downsteam_dists are defined from the following starting point coordinates_inertial_start = np.array( @@ -778,38 +809,35 @@ def set_grid(self) -> None: ) # Starting point in rotated coordinates - x_start, y_start, _, _, _ = rotate_coordinates_rel_west( + self.x_start, self.y_start, _, _, _ = rotate_coordinates_rel_west( self.wind_directions, coordinates_inertial_start, x_center_of_rotation=self.x_center_of_rotation, - y_center_of_rotation=self.y_center_of_rotation + y_center_of_rotation=self.y_center_of_rotation, ) - x_start, y_start = x_start[0,0,0], y_start[0,0,0] + self.x_start, self.y_start = self.x_start[0, 0, 0], self.y_start[0, 0, 0] downstream_dists_transpose = np.atleast_2d(self.downstream_dists).T - x = (x_start + downstream_dists_transpose) * np.ones((nprofiles, res)) + # The x-coordinate is fixed for every line (every row in `x`) + x = (self.x_start + downstream_dists_transpose) * np.ones((n_lines, res)) if self.direction == 'y': - y_single_profile = np.linspace( - y_start + self.profile_range[0], - y_start + self.profile_range[1], - res + y_single_line = np.linspace( + self.y_start + self.line_range[0], + self.y_start + self.line_range[1], + res, ) - y = y_single_profile * np.ones((nprofiles, res)) - z = self.reference_height * np.ones((nprofiles, res)) + y = y_single_line * np.ones((n_lines, res)) + z = self.reference_height * np.ones((n_lines, res)) elif self.direction == 'z': - z_min_profile = self.reference_height + self.profile_range[0] - if z_min_profile <= 0.0: - # Arbitrary small value to avoid possible errors at z = 0.0 - z_min_profile = 0.1 - z_single_profile = np.linspace( - z_min_profile, - self.reference_height + self.profile_range[1], - res + z_single_line = np.linspace( + self.reference_height + self.line_range[0], + self.reference_height + self.line_range[1], + res, ) - z = z_single_profile * np.ones((nprofiles, res)) - y = y_start * np.ones((nprofiles, res)) + z = z_single_line * np.ones((n_lines, res)) + y = self.y_start * np.ones((n_lines, res)) - self.x_sorted = x.flatten()[None,None,:,None,None] - self.y_sorted = y.flatten()[None,None,:,None,None] - self.z_sorted = z.flatten()[None,None,:,None,None] + self.x_sorted = x.flatten()[None, None, :, None, None] + self.y_sorted = y.flatten()[None, None, :, None, None] + self.z_sorted = z.flatten()[None, None, :, None, None] diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 619aacea8..b625a9786 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -959,29 +959,44 @@ def sample_velocity_deficit_profiles( wind_direction: float = None, homogeneous_wind_speed: float = None, ref_rotor_diameter: float = None, - x_inertial_start: float = None, - y_inertial_start: float = None, - reference_height: float = None + x_inertial_start: float = 0.0, + y_inertial_start: float = 0.0, + reference_height: float = None, ) -> list[pd.DataFrame]: """ - Extract velocity deficit profiles + Extract velocity deficit profiles at a set of downstream distances from a starting point + (usually a turbine location). For each downstream distance, a profile is sampled along + a line in either the cross-stream direction (y) or the vertical direction (z). + Velocity deficit is here defined as (homogeneous_wind_speed - u)/homogeneous_wind_speed, + where u is the wake velocity obtained when wind_shear = 0.0. Args: - downstream_dists: Streamwise location for each velocity profile. Default starting point - is at (x, y) = (0, 0). - - - reference_height: If 'direction' is y, then 'reference_height' defines the height of the - xy-plane in which the velocity profiles are sampled. If 'direction' is z, then the - velocity is sampled along the vertical direction with the 'profile_range' being - relative to the 'reference_height'. + direction: At each downstream location, this is the direction in which to sample the + profile. Either `y` or `z`. + downstream_dists: A list/array of streamwise locations for where to sample the profiles. + Default starting point in inertial coordinates is (0.0, 0.0, reference_height). + profile_range: Determines the extent of the line along which the profiles are sampled. + The range is defined about a point which lies some distance directly downstream of + the starting point. + resolution: Number of sample points in each profile. + wind_direction: A single wind direction. + homogeneous_wind_speed: A single wind speed. It is called homogeneous since 'wind_shear' + is temporarily set to 0.0 in this method. + ref_rotor_diameter: A reference rotor diameter which is used to normalize the + coordinates. + x_inertial_start: x-coordinate of starting point in the inertial coordinate system. + y_inertial_start: y-coordinate of starting point in the inertial coordinate system. + reference_height: If `direction` is y, then `reference_height` defines the height of the + xy-plane in which the velocity profiles are sampled. If `direction` is z, then the + velocity is sampled along the vertical direction with the `profile_range` being + relative to the `reference_height`. Returns: - 3DArrayFloat containing wind speed with dimensions - (# of wind directions, # of wind speeds, # of sample points) + A list of pandas DataFrame objects where each DataFrame represents one velocity deficit + profile. """ - if direction not in ('y', 'z'): - raise ValueError("'direction' must be either y or z") + if direction not in ['y', 'z']: + raise ValueError("`direction` must be either y or z.") if ref_rotor_diameter is None: unique_rotor_diameters = np.unique(self.floris.farm.rotor_diameters) @@ -989,8 +1004,10 @@ def sample_velocity_deficit_profiles( ref_rotor_diameter = unique_rotor_diameters[0] else: raise ValueError( - f"Multiple rotor_diameters detected: {unique_rotor_diameters}. " - "Provide a single 'ref_rotor_diameter'." + "Please provide a `ref_rotor_diameter`. This is needed to normalize the " + "coordinates. Could not select a value automatically since the number of " + "unique rotor diameters in the turbine layout is not 1. " + f"Found the following rotor diameters: {unique_rotor_diameters}." ) if downstream_dists is None: @@ -1008,38 +1025,33 @@ def sample_velocity_deficit_profiles( wind_direction = wind_directions_copy[0] else: raise ValueError( - "Multiple wind directions detected. Provide a single 'wind_direction' for " - "which to sample the velocity profiles." + "Could not determine a wind direction for which to sample the velocity " + "profiles. Either provide a single `wind_direction` as an argument to this " + "method, or initialize the Floris object with a single wind direction." ) if homogeneous_wind_speed is None: if len(wind_speeds_copy) == 1: homogeneous_wind_speed = wind_speeds_copy[0] self.logger.warning( - "'homogeneous_wind_speed' not provided. Setting it to the single wind speed " - "found in 'wind_speeds'. Note that the inflow is always homogeneous when " - "calculating the velocity deficit profiles. This is done by temporarily " - "setting 'wind_shear' to 0.0" + "`homogeneous_wind_speed` not provided. Setting it to the following wind speed " + f"found in the current flow field: {wind_speeds_copy[0]} m/s. Note that the " + "inflow is always homogeneous when calculating the velocity deficit profiles. " + "This is done by temporarily setting `wind_shear` to 0.0" ) else: raise ValueError( - "Multiple wind speeds detected. Provide a single 'homogeneous_wind_speed' for " - "which to sample the velocity profiles." + "Could not determine a wind speed for which to sample the velocity " + "profiles. Provide a single `homogeneous_wind_speed` to this method." ) - if x_inertial_start is None: - x_inertial_start = 0.0 - - if y_inertial_start is None: - y_inertial_start = 0.0 - if reference_height is None: reference_height = self.floris.flow_field.reference_wind_height self.reinitialize( wind_directions=[wind_direction], wind_speeds=[homogeneous_wind_speed], - wind_shear=0.0 + wind_shear=0.0, ) velocity_deficit_profiles = self.floris.solve_for_velocity_deficit_profiles( @@ -1051,13 +1063,13 @@ def sample_velocity_deficit_profiles( ref_rotor_diameter, x_inertial_start, y_inertial_start, - reference_height + reference_height, ) self.reinitialize( wind_directions=wind_directions_copy, wind_speeds=wind_speeds_copy, - wind_shear=wind_shear_copy + wind_shear=wind_shear_copy, ) return velocity_deficit_profiles diff --git a/floris/tools/visualization.py b/floris/tools/visualization.py index 3ae663dbd..55ecb6793 100644 --- a/floris/tools/visualization.py +++ b/floris/tools/visualization.py @@ -705,30 +705,43 @@ def calculate_horizontal_plane_with_turbines( @define class VelocityProfilesFigure(): """ - docstr + Create a figure which displays velocity deficit profiles at several downstream + locations of a turbine. + + Args: + downstream_dists_D: A list/array of streamwise locations at which the velocity deficit + profiles have been sampled. The locations should be normalized by the turbine + diameter D. + layout: A one- or two-element list defining the direction of the profiles and in which + order the directions are plotted. For example, ['y', 'z'] initializes a figure where + y-profiles are expected on the top row of Axes in the figure, and z-profiles are + expected on the bottom row. + ax_width: Roughly the width of each Axes. + ax_height: Roughly the height of each Axes. + """ downstream_dists_D: NDArrayFloat = field(converter=floris_array_converter) layout: list[str] = field(default=['y']) - width_per_col: float = field(default=2.07) - height_per_row: float = field(default=3.0) + ax_width: float = field(default=2.07) + ax_height: float = field(default=3.0) - nrows: int = field(init=False) - ncols: int = field(init=False) + n_rows: int = field(init=False) + n_cols: int = field(init=False) fig: plt.Figure = field(init=False) - axs: np.ndarray | plt.Axes = field(init=False) + axs: np.ndarray = field(init=False) - def __attrs_post_init__(self): - self.nrows = len(self.layout) - self.ncols = len(self.downstream_dists_D) - figsize = [0.7 + self.width_per_col * self.ncols, 1.0 + self.height_per_row * self.nrows] + def __attrs_post_init__(self) -> None: + self.n_rows = len(self.layout) + self.n_cols = len(self.downstream_dists_D) + figsize = [0.7 + self.ax_width * self.n_cols, 1.0 + self.ax_height * self.n_rows] self.fig, self.axs = plt.subplots( - self.nrows, - self.ncols, + self.n_rows, + self.n_cols, figsize=figsize, layout='tight', sharex='col', sharey='row', - squeeze=False + squeeze=False, ) for ax in self.axs[-1]: @@ -738,7 +751,7 @@ def __attrs_post_init__(self): for ax, x_D in zip(self.axs[0], self.downstream_dists_D): ax.set_title(f'$x/D = {x_D:.1f}$', fontsize=14) - for profile_direction, ax in zip(self.layout, self.axs[:,0]): + for ax, profile_direction in zip(self.axs[:,0], self.layout): ax.set_ylabel(f'${profile_direction}/D$', fontsize=14) ax.tick_params('y', labelsize=14) @@ -746,15 +759,26 @@ def __attrs_post_init__(self): def layout_validator(self, instance : attrs.Attribute, value : list[str]) -> None: allowed_layouts = [['y'], ['z'], ['y', 'z'], ['z', 'y']] if value not in allowed_layouts: - raise ValueError(f"'layout' must be one of the following: {allowed_layouts}") + raise ValueError(f"'layout' must be one of the following: {allowed_layouts}.") - def add_profiles(self, velocity_deficit_profiles, **kwargs): + def add_profiles( + self, + velocity_deficit_profiles: list[pd.DataFrame], + **kwargs + ) -> None: + """ + Add a list of velocity deficit profiles to the figure. Each profile is represented + as a pandas DataFrame. `kwargs` are passed to `ax.plot`. + """ for df in velocity_deficit_profiles: ax, profile_direction = self.match_profile_to_axes(df) profile_direction_D = f'{profile_direction}/D' ax.plot(df['velocity_deficit'], df[profile_direction_D], **kwargs) - def match_profile_to_axes(self, df): + def match_profile_to_axes( + self, + df: pd.DataFrame, + ) -> tuple[plt.Axes, str]: x_D = np.unique(df['x/D']) if len(x_D) == 1: x_D = x_D[0] @@ -777,7 +801,7 @@ def match_profile_to_axes(self, df): row = self.layout.index(profile_direction) col = None - for i in range(self.ncols): + for i in range(self.n_cols): if np.abs(x_D - self.downstream_dists_D[i]) < 0.001: col = i break @@ -785,34 +809,66 @@ def match_profile_to_axes(self, df): raise ValueError( "Could not add a velocity deficit profile at downstream distance " f"x/D = {x_D}. The downstream distance must be one of the following " - "values with which the instance of the VelocityProfilesFigure was initialized: " - f"{self.downstream_dists_D}" + "values with which this VelocityProfilesFigure object was initialized: " + f"{self.downstream_dists_D}." ) return self.axs[row,col], profile_direction - def set_xlim(self, xlim): + def set_xlim( + self, + xlim: list[float] | NDArrayFloat, + ) -> None: for ax in self.axs[-1]: ax.set_xlim(xlim) - def add_ref_lines_y_D(self, ref_lines, **kwargs): + def add_ref_lines_y( + self, + ref_lines_y_D: list[float] | NDArrayFloat, + **kwargs + ) -> None: + """ + Add reference lines to the VelocityProfilesFigure which go along the XAxis. + Commonly used to show the extent of the turbine. + Args: + ref_lines_y_D: A list of y-coordinates normalized by the turbine diameter D. + One coordinate per reference line. + **kwargs: Additional parameters to pass to `ax.plot`. + """ if 'y' not in self.layout: raise Exception( "Could not add reference lines to cross-stream (y) velocity profiles. No " "such profiles exist in the figure." ) row_y = self.layout.index('y') - self.add_ref_lines(ref_lines, row_y, **kwargs) + self.add_ref_lines(ref_lines_y_D, row_y, **kwargs) - def add_ref_lines_z_D(self, ref_lines, **kwargs): + def add_ref_lines_z( + self, + ref_lines_z_D: list[float] | NDArrayFloat, + **kwargs + ) -> None: + """ + Add reference lines to the VelocityProfilesFigure which go along the XAxis. + Commonly used to show the extent of the turbine. + Args: + ref_lines_z_D: A list of z-coordinates normalized by the turbine diameter D. + One coordinate per reference line. + **kwargs: Additional parameters to pass to `ax.plot`. + """ if 'z' not in self.layout: raise Exception( "Could not add reference lines to vertical (z) velocity profiles. No " "such profiles exist in the figure." ) row_z = self.layout.index('z') - self.add_ref_lines(ref_lines, row_z, **kwargs) + self.add_ref_lines(ref_lines_z_D, row_z, **kwargs) - def add_ref_lines(self, ref_lines, row, **kwargs): + def add_ref_lines( + self, + ref_lines_D: list[float] | NDArrayFloat, + row: int, + **kwargs + ) -> None: default_params = { 'linestyle': (0, (4, 2)), 'color': 'k', @@ -823,5 +879,5 @@ def add_ref_lines(self, ref_lines, row, **kwargs): kwargs[key] = default_params[key] for ax in self.axs[row]: - for coordinate in ref_lines: + for coordinate in ref_lines_D: ax.plot([0.0, 1.0], [coordinate, coordinate], **kwargs) From 2990e74777490e04408fd1b9bcb2ca4a274b7221 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 3 Aug 2023 15:44:56 +0200 Subject: [PATCH 25/46] Change default coefficient --- floris/simulation/wake_velocity/super_gaussian_vawt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 0f798e0e1..84f46223e 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -60,8 +60,8 @@ class SuperGaussianVAWTVelocityDeficit(BaseModel): :style: unsrt :filter: docname in docnames """ - wake_expansion_coeff_y: float = field(default=0.45) - wake_expansion_coeff_z: float = field(default=0.45) + wake_expansion_coeff_y: float = field(default=0.50) + wake_expansion_coeff_z: float = field(default=0.50) ay: float = field(default=0.95) az: float = field(default=4.5) by: float = field(default=0.35) From c75c517c9007a39bc3f5f313b0964d0a5c3ca03a Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 4 Aug 2023 14:04:41 +0200 Subject: [PATCH 26/46] Details in example --- examples/29_plot_velocity_deficit_profiles.py | 123 ++++++++++++------ 1 file changed, 84 insertions(+), 39 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index cb347d6f9..1bcbc51e0 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -22,7 +22,14 @@ """ -docstr +The first part of this example illustrates how to plot velocity deficit profiles at +several location downstream of a turbine. Here we use the following definition: + velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed + , where u is the wake velocity obtained when the incoming wind speed is the + same at all heights and equal to `homogeneous_wind_speed`. +The second part of the example shows a special case of how the profiles are affected +by a change in wind direction as well as a change in turbine location and sampling +starting point. """ def get_profiles(direction, resolution=100): @@ -34,7 +41,7 @@ def get_profiles(direction, resolution=100): ) if __name__ == '__main__': - D = 126.0 + D = 126.0 # Turbine diameter hub_height = 90.0 fi = FlorisInterface("inputs/gch.yaml") fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) @@ -42,20 +49,27 @@ def get_profiles(direction, resolution=100): downstream_dists = D * np.array([3, 5, 7]) homogeneous_wind_speed = 8.0 - # Velocity deficit profiles can be obtained in either the cross-stream (y) - # or the vertical direction (z). The default profile_range is [-2 * D, 2 * D] with - # the default origin being (x, y, z) = (0.0, 0.0, reference_wind_height). - # TODO: Explain difference between starting point and coordinates seen in fig + # Below, `get_profiles('y')` returns three velocity deficit profiles sampled along + # three lines that are all parallel to the y-axis (cross-stream direction). The + # streamwise location of each line in given in `downstream_dists`. + # Similarly, `get_profiles('z')` samples profiles in the vertical direction (z) at the + # same streamwise locations. profiles = get_profiles('y') + get_profiles('z') + # Initialize a VelocityProfilesFigure. The workflow is similar to a matplotlib Figure: + # Initialize it, plot data, and then customize it further if needed. + # The provided value of `layout` puts y-profiles on the top row of the figure and + # z-profiles on the bottom row. profiles_fig = VelocityProfilesFigure( downstream_dists_D=downstream_dists / D, layout=['y', 'z'], ) + # Add profiles to the figure. This method automatically determines the direction and + # streamwise location of each profile from the profile coordinates. profiles_fig.add_profiles(profiles, color='k') # Change velocity model to jensen, get the velocity deficit profiles, - # and add them to the figure + # and add them to the figure. floris_dict = fi.floris.as_dict() floris_dict['wake']['model_strings']['velocity_model'] = 'jensen' fi = FlorisInterface(floris_dict) @@ -71,40 +85,71 @@ def get_profiles(direction, resolution=100): profiles_fig.axs[0,0].legend(['gauss', 'jensen'], fontsize=11) profiles_fig.fig.suptitle( 'Velocity decifit profiles from different velocity models', - fontsize=14 + fontsize=14, ) - # Show that profiles are independent of the wind direction for a single turbine. - # This is because x is always in the streamwise direction. - # Also show that the coordinates x / D, y / D and z / D returned by - # sample_velocity_deficit_profiles are relative to the starting point mentioned above. - fi = FlorisInterface("inputs/gch.yaml") - fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) + # Second part of example: + # Case 1: Show that, if we have a single turbine, the profiles are independent of + # the wind direction. This is because x is defined to be in the streamwise direction + # in sample_velocity_deficit_profiles and VelocityProfilesFigure. + # Case 2: Show that the coordinates x/D, y/D and z/D returned by + # sample_velocity_deficit_profiles are relative to the sampling starting point. downstream_dists = D * np.array([3]) - profiles = get_profiles('y', resolution=400) - profiles_fig = VelocityProfilesFigure( - downstream_dists_D=downstream_dists / D, - layout=['y'] - ) - profiles_fig.add_profiles(profiles, color='k') - - fi.reinitialize(wind_directions=[315.0], layout_x=[D], layout_y=[D]) - profiles = fi.sample_velocity_deficit_profiles( - direction='y', - downstream_dists=downstream_dists, - resolution=21, - homogeneous_wind_speed=homogeneous_wind_speed, - x_inertial_start=D, - y_inertial_start=D, - ) - profiles_fig.add_profiles( - profiles, - linestyle='None', - marker='o', - color='b', - markerfacecolor='None', - markeredgewidth=1.3 - ) - profiles_fig.axs[0,0].legend(['WD = 270 deg', 'WD = 315 deg,\nand turbine\nmoved']) + for case in [1, 2]: + # Reference + fi = FlorisInterface("inputs/gch.yaml") + fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) + profiles = get_profiles('y', resolution=400) + profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y'], + ax_width=3.5, + ax_height=3.5, + ) + profiles_fig.add_profiles(profiles, color='k') + + if case == 1: + # Change wind direction compared to reference + wind_directions = [315.0] + layout_x, layout_y = [0.0], [0.0] + # This is the default starting point in sample_velocity_deficit_profiles + x_inertial_start, y_inertial_start = 0.0, 0.0 + elif case == 2: + # Change turbine location compared to reference. Then, set the sampling starting + # point to the new turbine location using the arguments + # `x_inertial_start` and `y_inertial_start`. + wind_directions = [270.0] + layout_x, layout_y = [D], [D] + x_inertial_start, y_inertial_start = D, D + + fi.reinitialize(wind_directions=wind_directions, layout_x=layout_x, layout_y=layout_y) + profiles = fi.sample_velocity_deficit_profiles( + direction='y', + downstream_dists=downstream_dists, + resolution=21, + homogeneous_wind_speed=homogeneous_wind_speed, + x_inertial_start=x_inertial_start, + y_inertial_start=y_inertial_start, + ) + profiles_fig.add_profiles( + profiles, + linestyle='None', + marker='o', + color='b', + markerfacecolor='None', + markeredgewidth=1.3, + ) + + if case == 1: + profiles_fig.axs[0,0].legend(['WD = 270 deg', 'WD = 315 deg']) + elif case == 2: + profiles_fig.fig.suptitle( + 'Legend (x, y) locations in inertial coordinates.\n' + 'x/D and y/D relative to sampling start point', + ) + profiles_fig.axs[0,0].legend([ + 'turbine location: (0, 0)\nsampling start point: (0, 0)', + 'turbine location: (D, D)\nsampling start point: (D, D)', + ]) plt.show() From 18270c596ae2ffa685be84fe371a00ad79cd76f7 Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 7 Aug 2023 12:28:14 +0200 Subject: [PATCH 27/46] Started to add some comments --- docs/references.bib | 11 ++++ .../wake_velocity/super_gaussian_vawt.py | 64 +++---------------- 2 files changed, 20 insertions(+), 55 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index f08ec7cd0..7bb0085df 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -271,3 +271,14 @@ @Article{bay_2022 URL = {https://wes.copernicus.org/preprints/wes-2022-17/}, DOI = {10.5194/wes-2022-17} } + +@article{ouro2021theoretical, + author = {Pablo Ouro and Maxime Lazennec}, + journal = {Flow}, + title = {Theoretical modelling of the three-dimensional wake of vertical axis turbines}, + year = {2021}, + pages = {E3}, + volume = {1}, + doi = {10.1017/flo.2021.4}, + publisher = {Cambridge University Press}, +} diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 84f46223e..ce9a8d53f 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -29,31 +29,14 @@ @define class SuperGaussianVAWTVelocityDeficit(BaseModel): """ - The Empirical Gauss velocity model has a Gaussian profile - (see :cite:`bastankhah2016experimental` and - :cite:`King2019Controls`) throughout and expands in a (smoothed) - piecewise linear fashion. - - parameter_dictionary (dict): Model-specific parameters. - Default values are used when a parameter is not included - in `parameter_dictionary`. Possible key-value pairs include: - - - **wake_expansion_rates** (*list*): List of expansion - rates for the Gaussian wake width. Must be of length 1 - or greater. - - **breakpoints_D** (*list*): List of downstream - locations, specified in terms of rotor diameters, where - the expansion rates go into effect. Must be one element - shorter than wake_expansion_rates. May be empty. - - **sigma_0_D** (*float*): Initial width of the Gaussian - wake at the turbine location, specified as a multiplier - of the rotor diameter. - - **smoothing_length_D** (*float*): Distance over which - the corners in the piece-wise linear wake expansion rate - are smoothed (specified as a multiplier of the rotor - diameter). - - **mixing_gain_deflection** (*float*): Gain to set the - increase in wake expansion due to wake-induced mixing. + This is a super-Gaussian wake velocity model for Vertical Axis Wind Turbines (VAWTs). + The model is based on :cite:`ouro2021theoretical` and allows the wake to have + different characteristics in the cross-stream and vertical direction. The initial + wake shape is closely related to the turbine cross section, which is: + turbine diameter * length of the vertical turbine blades. + + Parameters: + wake_expansion_coeff_y: References: .. bibliography:: /references.bib @@ -101,42 +84,13 @@ def function( z: np.ndarray, wind_veer: float ) -> None: - """ - Calculates the velocity deficits in the wake. - - Args: - x_i (np.array): Streamwise direction grid coordinates of - the ith turbine (m). - y_i (np.array): Cross stream direction grid coordinates of - the ith turbine (m). - z_i (np.array): Vertical direction grid coordinates of - the ith turbine (m) [not used]. - axial_induction_i (np.array): Axial induction factor of the - ith turbine (-) [not used]. - yaw_angle_i (np.array): Yaw angle of the ith turbine (deg). - ct_i (np.array): Thrust coefficient for the ith turbine (-). - hub_height_i (float): Hub height for the ith turbine (m). - rotor_diameter_i (np.array): Rotor diameter for the ith - turbine (m). - - x (np.array): Streamwise direction grid coordinates of the - flow field domain (m). - y (np.array): Cross stream direction grid coordinates of the - flow field domain (m). - z (np.array): Vertical direction grid coordinates of the - flow field domain (m). - wind_veer (np.array): Wind veer (deg). - - Returns: - np.array: Velocity deficits (-). - """ # Model parameters need to be unpacked for use in Numexpr ay, by, cy = self.ay, self.by, self.cy az, bz, cz = self.az, self.bz, self.cz # No specific near or far wakes in this model - downstream_mask = np.array(x > x_i + 0.1) + downstream_mask = (x > x_i + 0.1) # Nondimensional coordinates. Is z_i always equal to hub_height_i? x_tilde = ne.evaluate("downstream_mask * (x - x_i) / rotor_diameter_i") From 8165bcf218ac696ca0ddfbdb1aa0472415a543f3 Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 7 Aug 2023 13:09:01 +0200 Subject: [PATCH 28/46] Updated comments --- examples/29_plot_velocity_deficit_profiles.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/29_plot_velocity_deficit_profiles.py b/examples/29_plot_velocity_deficit_profiles.py index 1bcbc51e0..3e8a14e5c 100644 --- a/examples/29_plot_velocity_deficit_profiles.py +++ b/examples/29_plot_velocity_deficit_profiles.py @@ -16,7 +16,6 @@ import matplotlib.pyplot as plt import numpy as np -import floris.tools.visualization as wakeviz from floris.tools import FlorisInterface from floris.tools.visualization import VelocityProfilesFigure @@ -49,9 +48,9 @@ def get_profiles(direction, resolution=100): downstream_dists = D * np.array([3, 5, 7]) homogeneous_wind_speed = 8.0 - # Below, `get_profiles('y')` returns three velocity deficit profiles sampled along - # three lines that are all parallel to the y-axis (cross-stream direction). The - # streamwise location of each line in given in `downstream_dists`. + # Below, `get_profiles('y')` returns three velocity deficit profiles. These are sampled along + # three corresponding lines that are all parallel to the y-axis (cross-stream direction). + # The streamwise location of each line is given in `downstream_dists`. # Similarly, `get_profiles('z')` samples profiles in the vertical direction (z) at the # same streamwise locations. profiles = get_profiles('y') + get_profiles('z') @@ -89,14 +88,16 @@ def get_profiles(direction, resolution=100): ) # Second part of example: - # Case 1: Show that, if we have a single turbine, the profiles are independent of + # Case 1: Show that, if we have a single turbine, then the profiles are independent of # the wind direction. This is because x is defined to be in the streamwise direction # in sample_velocity_deficit_profiles and VelocityProfilesFigure. # Case 2: Show that the coordinates x/D, y/D and z/D returned by # sample_velocity_deficit_profiles are relative to the sampling starting point. + # By default, this starting point is at (0.0, 0.0, fi.floris.flow_field.reference_wind_height) + # in inertial coordinates. downstream_dists = D * np.array([3]) for case in [1, 2]: - # Reference + # The first added profile is a reference fi = FlorisInterface("inputs/gch.yaml") fi.reinitialize(layout_x=[0.0], layout_y=[0.0]) profiles = get_profiles('y', resolution=400) @@ -112,7 +113,7 @@ def get_profiles(direction, resolution=100): # Change wind direction compared to reference wind_directions = [315.0] layout_x, layout_y = [0.0], [0.0] - # This is the default starting point in sample_velocity_deficit_profiles + # Same as the default starting point but specified for completeness x_inertial_start, y_inertial_start = 0.0, 0.0 elif case == 2: # Change turbine location compared to reference. Then, set the sampling starting @@ -122,6 +123,7 @@ def get_profiles(direction, resolution=100): layout_x, layout_y = [D], [D] x_inertial_start, y_inertial_start = D, D + # Plot a second profile to show that it is equivalent to the reference fi.reinitialize(wind_directions=wind_directions, layout_x=layout_x, layout_y=layout_y) profiles = fi.sample_velocity_deficit_profiles( direction='y', From 4c8bf1e5909c9e62dcbb8902760822a9d49cd65e Mon Sep 17 00:00:00 2001 From: vallbog Date: Wed, 9 Aug 2023 00:33:33 +0200 Subject: [PATCH 29/46] Wake model comments --- .../wake_velocity/super_gaussian_vawt.py | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index ce9a8d53f..f7a777d9b 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -31,13 +31,19 @@ class SuperGaussianVAWTVelocityDeficit(BaseModel): """ This is a super-Gaussian wake velocity model for Vertical Axis Wind Turbines (VAWTs). The model is based on :cite:`ouro2021theoretical` and allows the wake to have - different characteristics in the cross-stream and vertical direction. The initial + different characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: turbine diameter * length of the vertical turbine blades. Parameters: - wake_expansion_coeff_y: - + wake_expansion_coeff_y: The wake expands linearly in y with a rate of + wake_expansion_coeff_y * turbulence_intensity_i. + wake_expansion_coeff_z: The wake expands linearly in z with a rate of + wake_expansion_coeff_z * turbulence_intensity_i. + ay, by, cy: Parameters that control how the shape function exponent `ny` evolves + in the streamwise direction. + az, bz, cz: Parameters that control how the shape function exponent `nz` evolves + in the streamwise direction. References: .. bibliography:: /references.bib :style: unsrt @@ -102,28 +108,27 @@ def function( ky_star = self.wake_expansion_coeff_y * turbulence_intensity_i kz_star = self.wake_expansion_coeff_z * turbulence_intensity_i - # Relative initial wake expansion + # `beta` is the initial wake expansion relative to the turbine cross-section fac = np.sqrt(1 - ct_i) beta = 0.5 * (1 + fac) / fac - # Initial wake width + # Characteristic nondimensional wake width at x - x_i = 0 ny_0 = ay + cy nz_0 = az + cz eta_0 = 1/ny_0 + 1/nz_0 epsilon = beta * ny_0 * nz_0 / (2 ** (2 * eta_0 + 2) * gamma(1 / ny_0) * gamma(1 / nz_0)) epsilon **= 1 / (2 * eta_0) - epsilon_y, epsilon_z = epsilon, epsilon - # Characteristic wake widths grow linearly in the streamwise direction - sigma_y_tilde = ne.evaluate("ky_star * x_tilde + epsilon_y") - sigma_z_tilde = ne.evaluate("kz_star * x_tilde_H + epsilon_z") + # Characteristic nondimensional wake widths grow linearly in the streamwise direction + sigma_y_tilde = ne.evaluate("ky_star * x_tilde + epsilon") + sigma_z_tilde = ne.evaluate("kz_star * x_tilde_H + epsilon") - # Gaussian exponents which determine the wake shape + # Eponents which determine the wake shape ny = ne.evaluate("ay * exp(-by * x_tilde) + cy") nz = ne.evaluate("az * exp(-bz * x_tilde_H) + cz") # At any streamwise location x behind the turbine, the velocity deficit in the - # y-z plane is given by multiplying the maximum deficit C = C(x) with two Gaussian + # yz-plane is given by multiplying the maximum deficit C = C(x) with two super-Gaussian # shape functions fy = fy(y) and fz = fz(z) eta = ne.evaluate("1 / ny + 1 / nz") ny_inv = ne.evaluate("1 / ny") From aab7c978a102ff7930c9f24a3a574e84c4818995 Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 11 Aug 2023 13:17:39 +0200 Subject: [PATCH 30/46] Added VAWT example --- examples/30_vertical_axis_wind_turbine.py | 132 ++++++++++++++++++ examples/inputs/super_gaussian_vawt.yaml | 53 +++++++ .../turbine_files/t1_long_blades_vawt.yaml | 28 ++++ .../inputs/turbine_files/t1_vawt_200kW.yaml | 28 ++++ .../wake_velocity/super_gaussian_vawt.py | 2 +- 5 files changed, 242 insertions(+), 1 deletion(-) create mode 100644 examples/30_vertical_axis_wind_turbine.py create mode 100644 examples/inputs/super_gaussian_vawt.yaml create mode 100644 examples/inputs/turbine_files/t1_long_blades_vawt.yaml create mode 100644 examples/inputs/turbine_files/t1_vawt_200kW.yaml diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py new file mode 100644 index 000000000..8d7881c2e --- /dev/null +++ b/examples/30_vertical_axis_wind_turbine.py @@ -0,0 +1,132 @@ +# 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 pandas as pd +import yaml +from matplotlib.ticker import AutoMinorLocator, MultipleLocator + +from floris.tools import FlorisInterface +from floris.tools.visualization import VelocityProfilesFigure, visualize_cut_plane + + +""" +This example shows a characteristic wake of a vertical-axis wind turbine (VAWT). The super-Gaussian +velocity model with default coefficients is used, which allows the wake to have different +characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape +is closely related to the turbine cross section, which is: + rotor diameter * length of the vertical turbine blades. +When plotting the velocity deficit profiles, we use the following definition: + velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed + , where u is the wake velocity obtained when the incoming wind speed is the + same at all heights and equal to `homogeneous_wind_speed`. +See example 29 for more details about how to sample and plot these kinds of profiles. +""" + +D = 26.0 # Rotor diameter +vawt_blade_length = 48.0 # Length of vertical turbine blades +hub_height = 40.0 +# Streamwise location of each velocity deficit profile +downstream_dists = D * np.array([1, 4, 7, 10]) +homogeneous_wind_speed = 7.0 + +fi = FlorisInterface('inputs/super_gaussian_vawt.yaml') + +# Velocity field in a horizontal plane +horizontal_plane = fi.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=hub_height, +) +horizontal_plane.df['x1'] /= D +horizontal_plane.df['x2'] /= D +fig, ax = plt.subplots() +visualize_cut_plane( + horizontal_plane, + ax, + title='Streamwise velocity [m/s] in a horizontal plane at hub height', +) +# The figure's XAxis and YAxis are in inertial coordinates which are not affected by +# the wind direction +ax.set_xlabel('$x_{inertial} / D$', fontsize=12) +ax.set_ylabel('$y_{inertial} / D$', fontsize=12) +ax.tick_params('both', labelsize=12) +fig.set_size_inches(8, 4) + +# Velocity field in a vertical plane +y_plane = fi.calculate_y_plane( + x_resolution=200, + z_resolution=100, + crossstream_dist=0.0, + z_bounds=np.array([-2, 2]) * D + hub_height, +) +y_plane.df['x1'] /= D +y_plane.df['x2'] /= D +fig, ax = plt.subplots() +visualize_cut_plane( + y_plane, + ax, + title='Streamwise velocity [m/s] in a xz-plane going through\n' + 'the center of the turbine', +) +# The figure's XAxis and YAxis are in inertial coordinates which are not affected by +# the wind direction +ax.set_xlabel('$x_{inertial} / D$', fontsize=12) +ax.set_ylabel('$z_{inertial} / D$', fontsize=12) +ax.tick_params('both', labelsize=12) +fig.set_size_inches(8, 4) + +# Velocity deficit profiles. +# The coordinates x/D, y/D and z/D returned by sample_velocity_deficit_profiles +# (and seen in the figure) are relative to the sampling starting point. +# Here, the following default starting point, in inertial coordinates, is used: +# (0.0, 0.0, fi.floris.flow_field.reference_wind_height) +profiles_y = fi.sample_velocity_deficit_profiles( + direction='y', + downstream_dists=downstream_dists, + homogeneous_wind_speed=homogeneous_wind_speed, +) +profiles_z = fi.sample_velocity_deficit_profiles( + direction='z', + downstream_dists=downstream_dists, + homogeneous_wind_speed=homogeneous_wind_speed, +) + +profiles_fig = VelocityProfilesFigure( + downstream_dists_D=downstream_dists / D, + layout=['y', 'z'], + ax_width=1.8, +) +profiles_fig.add_profiles(profiles_y + profiles_z, color='k') + +# Add dashed reference lines that show the extent of the turbine. +# Each line is defined by a coordinate that is normalized by `D`. +profiles_fig.add_ref_lines_y([-0.5, 0.5]) +H_half = vawt_blade_length / 2 +ref_lines_z_D = [-H_half / D, H_half / D] +profiles_fig.add_ref_lines_z(ref_lines_z_D) + +profiles_fig.set_xlim([0.0 - 0.05, 0.6 + 0.05]) +for ax in profiles_fig.axs[:,0]: + ax.set_ylim([-2 - 0.05, 2 + 0.05]) + +for ax in profiles_fig.axs[-1]: + ax.xaxis.set_major_locator(MultipleLocator(0.4)) + ax.xaxis.set_minor_locator(AutoMinorLocator(2)) + +profiles_fig.fig.suptitle('Velocity deficit profiles', fontsize=14) + +plt.show() diff --git a/examples/inputs/super_gaussian_vawt.yaml b/examples/inputs/super_gaussian_vawt.yaml new file mode 100644 index 000000000..0ce2f74f1 --- /dev/null +++ b/examples/inputs/super_gaussian_vawt.yaml @@ -0,0 +1,53 @@ +name: Super-Gaussian VAWT +description: Single vertical-axis wind turbine using super-Gaussian model +floris_version: v3.4 +logging: + console: + enable: true + level: WARNING + file: + enable: false + level: WARNING +solver: + type: turbine_grid + turbine_grid_points: 3 +farm: + layout_x: + - 0.0 + layout_y: + - 0.0 + turbine_type: + - !include 'turbine_files/t1_long_blades_vawt.yaml' +flow_field: + air_density: 1.225 + reference_wind_height: -1 + turbulence_intensity: 0.091 + wind_directions: + - 270.0 + wind_shear: 0.0 + wind_speeds: + - 7.0 + wind_veer: 0.0 +wake: + model_strings: + combination_model: fls + deflection_model: none + turbulence_model: none + velocity_model: super_gaussian_vawt + enable_secondary_steering: false + enable_yaw_added_recovery: false + enable_transverse_velocities: false + wake_deflection_parameters: + none: + wake_velocity_parameters: + super_gaussian_vawt: + wake_expansion_coeff_y: 0.50 + wake_expansion_coeff_z: 0.50 + ay: 0.95 + az: 4.5 + by: 0.35 + bz: 0.70 + cy: 2.4 + cz: 2.4 + wake_turbulence_parameters: + none: diff --git a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml new file mode 100644 index 000000000..ed5be569f --- /dev/null +++ b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml @@ -0,0 +1,28 @@ +# This is a vertical-axis wind turbine (VAWT). +turbine_type: 't1_long_blades_vawt' +generator_efficiency: 1.0 +hub_height: 40.0 +# Cosine exponent for power loss due to yaw misalignment. Assume that this value +# doesn't matter for VAWTs since they cannot have yaw. +pP: 1.88 +# Cosine exponent for power loss due to tilt. Placeholder value. +pT: 1.88 +# The overall geometry of a VAWT is characterized by the rotor diameter and +# the length of the vertical turbine blades. +rotor_diameter: 26.0 +vawt_blade_length: 48.0 +TSR: 3.8 +ref_density_cp_ct: 1.225 +ref_tilt_cp_ct: 0.0 +# Ct is assumed to be constant. This may be valid if the inflow velocity is close +# to the turbine's rated wind speed. The Cp value is a placeholder. +power_thrust_table: + power: + - 0.33 + - 0.33 + thrust: + - 0.64 + - 0.64 + wind_speed: + - 0.0 + - 50.0 diff --git a/examples/inputs/turbine_files/t1_vawt_200kW.yaml b/examples/inputs/turbine_files/t1_vawt_200kW.yaml new file mode 100644 index 000000000..01e66e265 --- /dev/null +++ b/examples/inputs/turbine_files/t1_vawt_200kW.yaml @@ -0,0 +1,28 @@ +# This is a vertical-axis wind turbine (VAWT). +turbine_type: 't1_vawt_200kW' +generator_efficiency: 1.0 +hub_height: 40.0 +# Cosine exponent for power loss due to yaw misalignment. Assume that this value +# doesn't matter for VAWTs since they cannot have yaw. +pP: 1.88 +# Cosine exponent for power loss due to tilt. Placeholder value. +pT: 1.88 +# The overall geometry of a VAWT is characterized by the rotor diameter and +# the length of the vertical turbine blades. +rotor_diameter: 26.0 +vawt_blade_length: 24.0 +TSR: 3.8 +ref_density_cp_ct: 1.225 +ref_tilt_cp_ct: 0.0 +# Cp and Ct are assumed to be constant. This may be valid if the inflow velocity is close +# to the turbine's rated wind speed. +power_thrust_table: + power: + - 0.33 + - 0.33 + thrust: + - 0.64 + - 0.64 + wind_speed: + - 0.0 + - 50.0 diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index f7a777d9b..8f74d678d 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -33,7 +33,7 @@ class SuperGaussianVAWTVelocityDeficit(BaseModel): The model is based on :cite:`ouro2021theoretical` and allows the wake to have different characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: - turbine diameter * length of the vertical turbine blades. + rotor diameter * length of the vertical turbine blades. Parameters: wake_expansion_coeff_y: The wake expands linearly in y with a rate of From ecb7fecc406ad0735965f125d87b7bdb909b8ee7 Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 13 Aug 2023 14:54:39 +0200 Subject: [PATCH 31/46] Details in example --- docs/references.bib | 11 ++ examples/30_vertical_axis_wind_turbine.py | 14 +- .../turbine_files/t1_long_blades_vawt.yaml | 151 +++++++++++++++++- .../inputs/turbine_files/t1_vawt_200kW.yaml | 28 ---- 4 files changed, 170 insertions(+), 34 deletions(-) delete mode 100644 examples/inputs/turbine_files/t1_vawt_200kW.yaml diff --git a/docs/references.bib b/docs/references.bib index 7bb0085df..b0c116117 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -282,3 +282,14 @@ @article{ouro2021theoretical doi = {10.1017/flo.2021.4}, publisher = {Cambridge University Press}, } + +@article{abkar2019theoretical, + author = {Mahdi Abkar}, + journal = {Energies}, + title = {Theoretical modeling of vertical-axis wind turbine wakes}, + year = {2019}, + number = {1}, + volume = {12}, + article-number = {10}, + doi = {10.3390/en12010010}, +} diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py index 8d7881c2e..4492880ef 100644 --- a/examples/30_vertical_axis_wind_turbine.py +++ b/examples/30_vertical_axis_wind_turbine.py @@ -24,16 +24,22 @@ """ -This example shows a characteristic wake of a vertical-axis wind turbine (VAWT). The super-Gaussian -velocity model with default coefficients is used, which allows the wake to have different -characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape -is closely related to the turbine cross section, which is: +This example shows a characteristic wake of a vertical-axis wind turbine (VAWT). It is based +on case 3 in :cite:``abkar2019theoretical. The super-Gaussian velocity model with default +coefficients is used, which allows the wake to have differenti characteristics in the +cross-stream (y) and vertical direction (z). The initial wake shape is closely related to +the turbine cross section, which is: rotor diameter * length of the vertical turbine blades. When plotting the velocity deficit profiles, we use the following definition: velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed , where u is the wake velocity obtained when the incoming wind speed is the same at all heights and equal to `homogeneous_wind_speed`. See example 29 for more details about how to sample and plot these kinds of profiles. + +References: + .. bibliography:: /references.bib + :style: unsrt + :filter: docname in docnames """ D = 26.0 # Rotor diameter diff --git a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml index ed5be569f..4ad7053fb 100644 --- a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml +++ b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml @@ -14,15 +14,162 @@ vawt_blade_length: 48.0 TSR: 3.8 ref_density_cp_ct: 1.225 ref_tilt_cp_ct: 0.0 -# Ct is assumed to be constant. This may be valid if the inflow velocity is close -# to the turbine's rated wind speed. The Cp value is a placeholder. +# Ct and Cp are assumed to be constant. This may be valid if the inflow velocity is close +# to the turbine's rated wind speed. The Cp value is a first guess. power_thrust_table: power: - 0.33 - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 + - 0.33 thrust: - 0.64 - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 + - 0.64 wind_speed: - 0.0 + - 2.0 + - 2.5 + - 3.0 + - 3.5 + - 4.0 + - 4.5 + - 5.0 + - 5.5 + - 6.0 + - 6.5 + - 7.0 + - 7.5 + - 8.0 + - 8.5 + - 9.0 + - 9.5 + - 10.0 + - 10.5 + - 11.0 + - 11.5 + - 12.0 + - 12.5 + - 13.0 + - 13.5 + - 14.0 + - 14.5 + - 15.0 + - 15.5 + - 16.0 + - 16.5 + - 17.0 + - 17.5 + - 18.0 + - 18.5 + - 19.0 + - 19.5 + - 20.0 + - 20.5 + - 21.0 + - 21.5 + - 22.0 + - 22.5 + - 23.0 + - 23.5 + - 24.0 + - 24.5 + - 25.0 + - 25.01 + - 25.02 - 50.0 diff --git a/examples/inputs/turbine_files/t1_vawt_200kW.yaml b/examples/inputs/turbine_files/t1_vawt_200kW.yaml deleted file mode 100644 index 01e66e265..000000000 --- a/examples/inputs/turbine_files/t1_vawt_200kW.yaml +++ /dev/null @@ -1,28 +0,0 @@ -# This is a vertical-axis wind turbine (VAWT). -turbine_type: 't1_vawt_200kW' -generator_efficiency: 1.0 -hub_height: 40.0 -# Cosine exponent for power loss due to yaw misalignment. Assume that this value -# doesn't matter for VAWTs since they cannot have yaw. -pP: 1.88 -# Cosine exponent for power loss due to tilt. Placeholder value. -pT: 1.88 -# The overall geometry of a VAWT is characterized by the rotor diameter and -# the length of the vertical turbine blades. -rotor_diameter: 26.0 -vawt_blade_length: 24.0 -TSR: 3.8 -ref_density_cp_ct: 1.225 -ref_tilt_cp_ct: 0.0 -# Cp and Ct are assumed to be constant. This may be valid if the inflow velocity is close -# to the turbine's rated wind speed. -power_thrust_table: - power: - - 0.33 - - 0.33 - thrust: - - 0.64 - - 0.64 - wind_speed: - - 0.0 - - 50.0 From dd19720d5fd1e58f4783303c59c27d3df564eae7 Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 14 Aug 2023 16:54:24 +0200 Subject: [PATCH 32/46] Draft of VerticalAxisTurbine class --- floris/simulation/__init__.py | 10 +- floris/simulation/farm.py | 13 ++- floris/simulation/turbine.py | 168 ++++++++++++++++++++++++++++++++++ 3 files changed, 189 insertions(+), 2 deletions(-) diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 07aca3cf9..dffea6fb7 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -37,7 +37,15 @@ import floris.logging_manager from .base import BaseClass, BaseModel, State -from .turbine import average_velocity, axial_induction, Ct, power, rotor_effective_velocity, Turbine +from .turbine import ( + average_velocity, + axial_induction, + Ct, + power, + rotor_effective_velocity, + Turbine, + VerticalAxisTurbine, +) from .farm import Farm from .grid import ( FlowFieldGrid, diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index 956056de7..e109a2ff9 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -24,6 +24,7 @@ BaseClass, State, Turbine, + VerticalAxisTurbine, ) from floris.simulation.turbine import compute_tilt_angles_for_floating_turbines from floris.type_dec import ( @@ -244,7 +245,17 @@ def construct_turbine_correct_cp_ct_for_tilt(self): ) def construct_turbine_map(self): - self.turbine_map = [Turbine.from_dict(turb) for turb in self.turbine_definitions] + self.turbine_map = [] + for turb in self.turbine_definitions: + try: + is_vertical_axis_turbine = turb['is_vertical_axis_turbine'] + except KeyError: + is_vertical_axis_turbine = False + + if is_vertical_axis_turbine: + self.turbine_map.append(VerticalAxisTurbine.from_dict(turb)) + else: + self.turbine_map.append(Turbine.from_dict(turb)) def construct_turbine_fCts(self): self.turbine_fCts = { diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index e7a312917..cc40527de 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -627,6 +627,9 @@ class Turbine(BaseClass): power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) floating_tilt_table = field(default=None) floating_correct_cp_ct_for_tilt = field(default=None) + # Placeholder attributes for a vertical-axis wind turbine. If using such a turbine, + # please instantiate a VerticalAxisTurbine instead of a Turbine. + is_vertical_axis_turbine = field(default=False) vawt_blade_length: float = field(default=0.0) # rloc: float = float_attrib() # TODO: goes here or on the Grid? @@ -760,3 +763,168 @@ def check_for_cp_ct_correct_flag_if_floating( "If a floating tilt/wind_speed table is defined, the boolean flag" "floating_correct_cp_ct_for_tilt must also be defined." ) + +@define +class VerticalAxisTurbine(BaseClass): + """ + Turbine is a class containing objects pertaining to the individual + turbines. + + Turbine is a model class representing a particular wind turbine. It + is largely a container of data and parameters, but also contains + methods to probe properties for output. + + Parameters: + rotor_diameter (:py:obj: float): The rotor diameter (m). + hub_height (:py:obj: float): The hub height (m). + pP (:py:obj: float): The cosine exponent relating the yaw + misalignment angle to power. + pT (:py:obj: float): The cosine exponent relating the rotor + tilt angle to power. + generator_efficiency (:py:obj: float): The generator + efficiency factor used to scale the power production. + ref_density_cp_ct (:py:obj: float): The density at which the provided + cp and ct is defined + power_thrust_table (PowerThrustTable): A dictionary containing the + following key-value pairs: + + power (:py:obj: List[float]): The coefficient of power at + different wind speeds. + thrust (:py:obj: List[float]): The coefficient of thrust + at different wind speeds. + wind_speed (:py:obj: List[float]): The wind speeds for + which the power and thrust values are provided (m/s). + """ + + turbine_type: str = field() + rotor_diameter: float = field() + vawt_blade_length: float = field() + hub_height: float = field() + pP: float = field() # Remove if possible + pT: float = field() + TSR: float = field() + generator_efficiency: float = field() + ref_density_cp_ct: float = field() + ref_tilt_cp_ct: float = field() + power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) + floating_tilt_table = field(default=None) + floating_correct_cp_ct_for_tilt = field(default=None) + is_vertical_axis_turbine = field(default=True) + + # Initialized in the post_init function + rotor_radius: float = field(init=False) + turbine_yz_cross_section_area: float = field(init=False) + fCp_interp: interp1d = field(init=False) + fCt_interp: interp1d = field(init=False) + power_interp: interp1d = field(init=False) + tilt_interp: interp1d = field(init=False) + + def __attrs_post_init__(self) -> None: + + # Post-init initialization for the power curve interpolation functions + wind_speeds = self.power_thrust_table.wind_speed + self.fCp_interp = interp1d( + wind_speeds, + self.power_thrust_table.power, + fill_value=(0.0, 1.0), + bounds_error=False, + ) + inner_power = ( + 0.5 * self.turbine_yz_cross_section_area + * self.fCp_interp(wind_speeds) + * self.generator_efficiency + * wind_speeds ** 3 + ) + self.power_interp = interp1d( + wind_speeds, + inner_power + ) + + """ + Given an array of wind speeds, this function returns an array of the + interpolated thrust coefficients from the power / thrust table used + to define the Turbine. The values are bound by the range of the input + values. Any requested wind speeds outside of the range of input wind + speeds are assigned Ct of 0.0001 or 0.9999. + + The fill_value arguments sets (upper, lower) bounds for any values + outside of the input range. + """ + self.fCt_interp = interp1d( + wind_speeds, + self.power_thrust_table.thrust, + fill_value=(0.0001, 0.9999), + bounds_error=False, + ) + + # If defined, create a tilt interpolation function for floating turbines. + # fill_value currently set to apply the min or max tilt angles if outside + # of the interpolation range. + if self.floating_tilt_table is not None: + self.floating_tilt_table = TiltTable.from_dict(self.floating_tilt_table) + self.fTilt_interp = interp1d( + self.floating_tilt_table.wind_speeds, + self.floating_tilt_table.tilt, + fill_value=(0.0, self.floating_tilt_table.tilt[-1]), + bounds_error=False, + ) + self.tilt_interp = self.fTilt_interp + self.correct_cp_ct_for_tilt = self.floating_correct_cp_ct_for_tilt + else: + self.fTilt_interp = None + self.tilt_interp = None + self.correct_cp_ct_for_tilt = False + + @rotor_diameter.validator + def reset_rotor_diameter_dependencies(self, instance: attrs.Attribute, value: float) -> None: + """Resets the `turbine_yz_cross_section_area` and `rotor_radius` attributes.""" + self.turbine_yz_cross_section_area = self.rotor_diameter * self.vawt_blade_length + # Temporarily turn off validators to avoid infinite recursion + with attrs.validators.disabled(): + # Reset the value + self.rotor_radius = value / 2.0 + + @rotor_radius.validator + def reset_rotor_radius_dependencies(self, instance: attrs.Attribute, value: float) -> None: + """ + Resets the `rotor_diameter` attribute to trigger the recalculation of + `turbine_yz_cross_section_area`. + """ + self.rotor_diameter = value * 2.0 + + @vawt_blade_length.validator + def reset_vawt_blade_length_dependencies(self, instance: attrs.Attribute, value: float) -> None: + """Resets the `turbine_yz_cross_section_area` attribute.""" + self.turbine_yz_cross_section_area = self.rotor_diameter * self.vawt_blade_length + + @floating_tilt_table.validator + def check_floating_tilt_table(self, instance: attrs.Attribute, value: Any) -> None: + """ + Check that if the tile/wind_speed table is defined, that the tilt and + wind_speed arrays are the same length so that the interpolation will work. + """ + if self.floating_tilt_table is not None: + if ( + len(self.floating_tilt_table["tilt"]) + != len(self.floating_tilt_table["wind_speeds"]) + ): + raise ValueError( + "tilt and wind_speeds must be the same length for the interpolation to work." + ) + + @floating_correct_cp_ct_for_tilt.validator + def check_for_cp_ct_correct_flag_if_floating( + self, + instance: attrs.Attribute, + value: Any + ) -> None: + """ + Check that the boolean flag exists for correcting Cp/Ct for tilt + if a tile/wind_speed table is also defined. + """ + if self.floating_tilt_table is not None: + if self.floating_correct_cp_ct_for_tilt is None: + raise ValueError( + "If a floating tilt/wind_speed table is defined, the boolean flag" + "floating_correct_cp_ct_for_tilt must also be defined." + ) From 9b50f3dd664ae94f611f2b208c38805ad67ff1c8 Mon Sep 17 00:00:00 2001 From: vallbog Date: Tue, 15 Aug 2023 15:02:31 +0200 Subject: [PATCH 33/46] Made TurbineGrid work with mix of HAWTs and VAWTs --- floris/simulation/farm.py | 9 ++++ floris/simulation/floris.py | 3 ++ floris/simulation/grid.py | 83 ++++++++++++++++++++++++-------- floris/tools/floris_interface.py | 30 +++++++++--- 4 files changed, 97 insertions(+), 28 deletions(-) diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index e109a2ff9..fd6af012b 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -84,6 +84,7 @@ class Farm(BaseClass): correct_cp_ct_for_tilt: NDArrayFloat = field(init=False, default=[]) correct_cp_ct_for_tilt_sorted: NDArrayFloat = field(init=False, default=[]) turbine_fTilts: list = field(init=False, default=[]) + is_vertical_axis_turbine: list[bool] = field(init=False, default=[]) def __attrs_post_init__(self) -> None: # Turbine definitions can be supplied in three ways: @@ -257,6 +258,14 @@ def construct_turbine_map(self): else: self.turbine_map.append(Turbine.from_dict(turb)) + def construct_is_vertical_axis_turbine(self): + self.is_vertical_axis_turbine = [] + for turb in self.turbine_definitions: + try: + self.is_vertical_axis_turbine.append(turb['is_vertical_axis_turbine']) + except KeyError: + self.is_vertical_axis_turbine.append(False) + def construct_turbine_fCts(self): self.turbine_fCts = { turb.turbine_type: turb.fCt_interp for turb in self.turbine_map diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 8a20cbdfa..b7051b683 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -91,6 +91,7 @@ def __attrs_post_init__(self) -> None: self.farm.construct_turbine_fTilts() self.farm.construct_turbine_correct_cp_ct_for_tilt() self.farm.construct_coordinates() + self.farm.construct_is_vertical_axis_turbine() self.farm.set_yaw_angles(self.flow_field.n_wind_directions, self.flow_field.n_wind_speeds) self.farm.set_tilt_to_ref_tilt( self.flow_field.n_wind_directions, @@ -105,6 +106,8 @@ def __attrs_post_init__(self) -> None: wind_speeds=self.flow_field.wind_speeds, grid_resolution=self.solver["turbine_grid_points"], time_series=self.flow_field.time_series, + is_vertical_axis_turbine=self.farm.is_vertical_axis_turbine, + vawt_blade_lengths=self.farm.vawt_blade_lengths, ) elif self.solver["type"] == "turbine_cubature_grid": self.grid = TurbineCubatureGrid( diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 4c9bd6954..a99dc24fb 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -118,7 +118,7 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter # TODO move this to the grid types and off of the base class """Check that grid resolution is given as int or Vec3 with int components.""" if isinstance(value, int) and \ - isinstance(self, (TurbineGrid, TurbineCubatureGrid, PointsGrid, LinesGrid)): + isinstance(self, (TurbineCubatureGrid, PointsGrid, LinesGrid)): return elif isinstance(value, Iterable) and isinstance(self, FlowFieldPlanarGrid): assert type(value[0]) is int @@ -127,6 +127,22 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter assert type(value[0]) is int assert type(value[1]) is int assert type(value[2]) is int + elif isinstance(self, TurbineGrid): + # Always store `grid_resolution` as an Iterable(int,) of length 2 when using the + # TurbineGrid class. This is the resolution in the cross-stream and vertical direction + # respectively. + if isinstance(value, int): + self.grid_resolution = [value, value] + elif isinstance(value, Iterable) and len(value) == 1: + assert type(value[0]) is int + self.grid_resolution = [value[0], value[0]] + elif isinstance(value, Iterable) and len(value) == 2: + assert type(value[0]) is int + assert type(value[1]) is int + else: + raise ValueError( + '`grid_resolution` must be of type int or Iterable(int,) of length 1 or 2.' + ) else: raise TypeError("`grid_resolution` must be of type int or Iterable(int,)") @@ -150,6 +166,8 @@ class TurbineGrid(Grid): series. """ # TODO: describe these and the differences between `sorted_indices` and `sorted_coord_indices` + is_vertical_axis_turbine: np.ndarray = field(converter=np.array) + vawt_blade_lengths: NDArrayFloat = field() sorted_indices: NDArrayInt = field(init=False) sorted_coord_indices: NDArrayInt = field(init=False) unsorted_indices: NDArrayInt = field(init=False) @@ -225,43 +243,68 @@ def set_grid(self) -> None: # the rotor radius. # Defaults to 0.5. - # Create the data for the turbine grids - radius_ratio = 0.5 - disc_area_radius = radius_ratio * self.reference_turbine_diameter / 2 + width_ratio = 0.5 + height_ratio = 0.5 + + # Use a square area for traditional horisontal-axis turbines and a rectangular area + # for vertical-axis turbines + width = width_ratio * self.reference_turbine_diameter + height = self.vawt_blade_lengths * self.is_vertical_axis_turbine + height += self.reference_turbine_diameter * np.invert(self.is_vertical_axis_turbine) + height *= height_ratio + template_grid = np.ones( ( self.n_wind_directions, self.n_wind_speeds, self.n_turbines, - self.grid_resolution, - self.grid_resolution, + self.grid_resolution[0], + self.grid_resolution[1], ), dtype=floris_float_type ) - # Calculate the radial distance from the center of the turbine rotor. - # If a grid resolution of 1 is selected, create a disc_grid of zeros, as - # np.linspace would just return the starting value of -1 * disc_area_radius - # which would place the point below the center of the rotor. - if self.grid_resolution == 1: - disc_grid = np.zeros((np.shape(disc_area_radius)[0], 1 )) + + # Conceptually, each turbine's grid can be seen as a meshgrid of a coordinate vector in + # `cross_stream_vecs` and a coordinate vector in `vertical_vecs`. These variables store one + # vector per turbine, since the turbine diameters or vawt blade lenths may be different. + if self.grid_resolution[0] == 1: + # If the first grid resolution is 1, then let each coordinate vector in `cross_stream_vecs` + # be [0.0], as np.linspace would just return the starting value of -width / 2 which would + # place the grid point(s) to the side of the center of the turbine. + cross_stream_vecs = np.zeros((self.n_turbines, 1)) + else: + cross_stream_vecs = np.linspace( + -width / 2, + width / 2, + self.grid_resolution[0], + dtype=floris_float_type, + axis=1, + ) + + if self.grid_resolution[1] == 1: + # If the second grid resolution is 1, then let each coordinate vector in `vertical_vecs` + # be [0.0], as np.linspace would just return the starting value of -height / 2 which would + # place the grid point(s) below the center of the turbine. + vertical_vecs = np.zeros((self.n_turbines, 1)) else: - disc_grid = np.linspace( - -1 * disc_area_radius, - disc_area_radius, - self.grid_resolution, + vertical_vecs = np.linspace( + -height / 2, + height / 2, + self.grid_resolution[1], dtype=floris_float_type, - axis=1 + axis=1, ) + # Construct the turbine grids # Here, they are already rotated to the correct orientation for each wind direction _x = x[:, :, :, None, None] * template_grid ones_grid = np.ones( - (self.n_turbines, self.grid_resolution, self.grid_resolution), + (self.n_turbines, self.grid_resolution[0], self.grid_resolution[1]), dtype=floris_float_type ) - _y = y[:, :, :, None, None] + template_grid * ( disc_grid[None, None, :, :, None]) - _z = z[:, :, :, None, None] + template_grid * ( disc_grid[:, None, :] * ones_grid ) + _y = y[:, :, :, None, None] + template_grid * ( cross_stream_vecs[None, None, :, :, None]) + _z = z[:, :, :, None, None] + template_grid * ( vertical_vecs[:, None, :] * ones_grid ) # Sort the turbines at each wind direction diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index b625a9786..ab117612e 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -15,6 +15,7 @@ from __future__ import annotations from pathlib import Path +from typing import Iterable import numpy as np import pandas as pd @@ -79,15 +80,28 @@ def __init__(self, configuration: dict | str | Path): ) self.logger.warning(err_msg, stack_info=True) - # Check the turbine_grid_points is reasonable + # Check that turbine_grid_points is reasonable if self.floris.solver["type"] == "turbine_grid": - if self.floris.solver["turbine_grid_points"] > 3: - self.logger.error( - f"turbine_grid_points value is {self.floris.solver['turbine_grid_points']} " - "which is larger than the recommended value of less than or equal to 3. " - "High amounts of turbine grid points reduce the computational performance " - "but have a small change on accuracy." - ) + ngrid = self.floris.solver["turbine_grid_points"] + if np.any(np.array(ngrid) > 3): + if isinstance(ngrid, int) or \ + (isinstance(ngrid, Iterable) and len(np.array(ngrid)) == 1): + self.logger.error( + f"`turbine_grid_points` is {ngrid} (which is an alias for grid resolution " + f"{self.floris.grid.grid_resolution} in the cross-stream and vertical " + "direction respectively). This is larger then the recommended value of " + "less than or equal to 3 in each direction. High amounts of turbine grid " + "points reduce the computational performance but have a small change on " + "accuracy." + ) + else: + self.logger.error( + f"`turbine_grid_points` is {ngrid}. This is larger then the recommended " + "value of less than or equal to 3 in each direction. High amounts of " + "turbine grid points reduce the computational performance but have a small " + "change on accuracy." + ) + raise ValueError("turbine_grid_points must be less than or equal to 3.") def assign_hub_height_to_ref_height(self): From 36b95d799beced383544cd2925ed3cfb37d83788 Mon Sep 17 00:00:00 2001 From: vallbog Date: Tue, 15 Aug 2023 21:18:03 +0200 Subject: [PATCH 34/46] Power calculation works for VAWTs --- examples/30_vertical_axis_wind_turbine.py | 44 +++++++++++++++++-- .../turbine_files/t1_long_blades_vawt.yaml | 2 +- floris/simulation/floris.py | 1 + floris/simulation/grid.py | 5 +++ floris/simulation/solver.py | 16 +++++++ floris/tools/floris_interface.py | 6 +++ 6 files changed, 69 insertions(+), 5 deletions(-) diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py index 4492880ef..ee7037c32 100644 --- a/examples/30_vertical_axis_wind_turbine.py +++ b/examples/30_vertical_axis_wind_turbine.py @@ -20,13 +20,17 @@ from matplotlib.ticker import AutoMinorLocator, MultipleLocator from floris.tools import FlorisInterface -from floris.tools.visualization import VelocityProfilesFigure, visualize_cut_plane +from floris.tools.visualization import ( + plot_rotor_values, + VelocityProfilesFigure, + visualize_cut_plane, +) """ -This example shows a characteristic wake of a vertical-axis wind turbine (VAWT). It is based -on case 3 in :cite:``abkar2019theoretical. The super-Gaussian velocity model with default -coefficients is used, which allows the wake to have differenti characteristics in the +The first part of this example shows a characteristic wake of a vertical-axis wind turbine (VAWT). +It is based on case 3 in :cite:``abkar2019theoretical. The super-Gaussian velocity model with +default coefficients is used, which allows the wake to have differenti characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: rotor diameter * length of the vertical turbine blades. @@ -36,6 +40,9 @@ same at all heights and equal to `homogeneous_wind_speed`. See example 29 for more details about how to sample and plot these kinds of profiles. +The second part of the example shows how turbine velocities and turbine powers are obtained in a +layout with two VAWTs. + References: .. bibliography:: /references.bib :style: unsrt @@ -135,4 +142,33 @@ profiles_fig.fig.suptitle('Velocity deficit profiles', fontsize=14) +# Switch to a two-turbine layout. Then, calculate the velocities at each turbine in a +# yz-grid centered on the hub of the turbine. The grids are based on the VAWTs' rectangular +# cross-section and have a resolution of 3x5 to make them reasonably equidistant. +solver_settings = { + 'type': 'turbine_grid', + 'turbine_grid_points': [3, 5] +} +fi.reinitialize(layout_x=[0.0, 78.0], layout_y=[0.0, 0.0], solver_settings=solver_settings) +fi.calculate_wake() + +# Plot the velocities in each grid +fig, axes, _, _ = plot_rotor_values( + fi.floris.flow_field.u, + wd_index=0, + ws_index=0, + n_rows=1, + n_cols=2, + return_fig_objects=True, +) +fig.suptitle( + 'Streamwise velocity [m/s] in yz-grids centered on\n' + 'the hub of each turbine, in a two-turbine layout' +) + +# Calculate the power generated by each turbine. The average velocity in each yz-grid +# is used in the calculation. +turbine_powers = fi.get_turbine_powers() / 1000 +print(f'Turbine powers for a two-turbine layout [kW] =\n{turbine_powers}') + plt.show() diff --git a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml index 4ad7053fb..917896254 100644 --- a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml +++ b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml @@ -1,5 +1,5 @@ -# This is a vertical-axis wind turbine (VAWT). turbine_type: 't1_long_blades_vawt' +is_vertical_axis_turbine: True generator_efficiency: 1.0 hub_height: 40.0 # Cosine exponent for power loss due to yaw misalignment. Assume that this value diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index b7051b683..a1870625f 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -117,6 +117,7 @@ def __attrs_post_init__(self) -> None: wind_speeds=self.flow_field.wind_speeds, time_series=self.flow_field.time_series, grid_resolution=self.solver["turbine_grid_points"], + is_vertical_axis_turbine=self.farm.is_vertical_axis_turbine, ) elif self.solver["type"] == "flow_field_grid": self.grid = FlowFieldGrid( diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index a99dc24fb..c5ed6903b 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -353,6 +353,7 @@ class TurbineCubatureGrid(Grid): time_series (:py:obj:`bool`): Flag to indicate whether the supplied wind data is a time series. """ + is_vertical_axis_turbine: np.ndarray = field(converter=np.array) sorted_indices: NDArrayInt = field(init=False) sorted_coord_indices: NDArrayInt = field(init=False) unsorted_indices: NDArrayInt = field(init=False) @@ -362,6 +363,10 @@ class TurbineCubatureGrid(Grid): def __attrs_post_init__(self) -> None: super().__attrs_post_init__() + if np.any(self.is_vertical_axis_turbine): + raise NotImplementedError( + '`TurbineCubatureGrid` does not support vertical-axis turbines.' + ) self.set_grid() def set_grid(self) -> None: diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index c2ac6818e..c83549438 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -279,6 +279,7 @@ def full_flow_sequential_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -288,6 +289,8 @@ def full_flow_sequential_solver( wind_speeds=turbine_grid_flow_field.wind_speeds, grid_resolution=3, time_series=turbine_grid_flow_field.time_series, + is_vertical_axis_turbine=turbine_grid_farm.is_vertical_axis_turbine, + vawt_blade_lengths=turbine_grid_farm.vawt_blade_lengths, ) turbine_grid_farm.expand_farm_properties( turbine_grid_flow_field.n_wind_directions, @@ -681,6 +684,7 @@ def full_flow_cc_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -690,6 +694,8 @@ def full_flow_cc_solver( wind_speeds=turbine_grid_flow_field.wind_speeds, grid_resolution=3, time_series=turbine_grid_flow_field.time_series, + is_vertical_axis_turbine=turbine_grid_farm.is_vertical_axis_turbine, + vawt_blade_lengths=turbine_grid_farm.vawt_blade_lengths, ) turbine_grid_farm.expand_farm_properties( turbine_grid_flow_field.n_wind_directions, @@ -1113,6 +1119,7 @@ def full_flow_turbopark_solver( # turbine_grid_farm.construct_turbine_TSRs() # turbine_grid_farm.construc_turbine_pPs() # turbine_grid_farm.construct_coordinates() + # turbine_grid_farm.construct_is_vertical_axis_turbine() # turbine_grid = TurbineGrid( @@ -1121,6 +1128,9 @@ def full_flow_turbopark_solver( # wind_directions=turbine_grid_flow_field.wind_directions, # wind_speeds=turbine_grid_flow_field.wind_speeds, # grid_resolution=11, + # is_vertical_axis_turbine=turbine_grid_farm.is_vertical_axis_turbine, + # vawt_blade_lengths=turbine_grid_farm.vawt_blade_lengths, + # ) # turbine_grid_farm.expand_farm_properties( # turbine_grid_flow_field.n_wind_directions, @@ -1358,6 +1368,7 @@ def full_flow_empirical_gauss_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -1367,6 +1378,8 @@ def full_flow_empirical_gauss_solver( wind_speeds=turbine_grid_flow_field.wind_speeds, grid_resolution=3, time_series=turbine_grid_flow_field.time_series, + is_vertical_axis_turbine=turbine_grid_farm.is_vertical_axis_turbine, + vawt_blade_lengths=turbine_grid_farm.vawt_blade_lengths, ) turbine_grid_farm.expand_farm_properties( turbine_grid_flow_field.n_wind_directions, @@ -1625,6 +1638,7 @@ def full_flow_vawt_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -1634,6 +1648,8 @@ def full_flow_vawt_solver( wind_speeds=turbine_grid_flow_field.wind_speeds, grid_resolution=3, time_series=turbine_grid_flow_field.time_series, + is_vertical_axis_turbine=turbine_grid_farm.is_vertical_axis_turbine, + vawt_blade_lengths=turbine_grid_farm.vawt_blade_lengths, ) turbine_grid_farm.expand_farm_properties( turbine_grid_flow_field.n_wind_directions, diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index ab117612e..e308d8016 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -104,6 +104,12 @@ def __init__(self, configuration: dict | str | Path): raise ValueError("turbine_grid_points must be less than or equal to 3.") + if len(np.unique(self.floris.farm.is_vertical_axis_turbine)) == 2: + raise NotImplementedError( + 'There is currently no solver that supports a mix of horisontal-axis turbines ' + 'and vertical-axis turbines in the same farm.' + ) + def assign_hub_height_to_ref_height(self): # Confirm can do this operation From db44bb692492a3f1abb0d4914c5e295d0778cc26 Mon Sep 17 00:00:00 2001 From: vallbog Date: Tue, 15 Aug 2023 22:41:17 +0200 Subject: [PATCH 35/46] Bug fix related to new def of grid_resolution, and updated test conf --- floris/simulation/solver.py | 8 ++++---- floris/simulation/turbine.py | 2 +- tests/conftest.py | 6 +++++- tests/flow_field_unit_test.py | 4 ++-- 4 files changed, 12 insertions(+), 8 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index c83549438..d73ef4665 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -223,7 +223,7 @@ def sequential_solver( # Calculate wake overlap for wake-added turbulence (WAT) area_overlap = ( np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4)) - / (grid.grid_resolution * grid.grid_resolution) + / (grid.grid_resolution[0] * grid.grid_resolution[1]) ) area_overlap = area_overlap[:, :, :, None, None] @@ -630,7 +630,7 @@ def cc_solver( # Calculate wake overlap for wake-added turbulence (WAT) area_overlap = 1 - ( np.sum(turb_u_wake <= 0.05, axis=(3, 4)) - / (grid.grid_resolution * grid.grid_resolution) + / (grid.grid_resolution[0] * grid.grid_resolution[1]) ) area_overlap = area_overlap[:, :, :, None, None] @@ -1063,7 +1063,7 @@ def turbopark_solver( # Calculate wake overlap for wake-added turbulence (WAT) area_overlap = ( np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4)) - / (grid.grid_resolution * grid.grid_resolution) + / (grid.grid_resolution[0] * grid.grid_resolution[1]) ) area_overlap = area_overlap[:, :, :, None, None] @@ -1320,7 +1320,7 @@ def empirical_gauss_solver( # Calculate wake overlap for wake-added turbulence (WAT) area_overlap = np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4))\ - / (grid.grid_resolution * grid.grid_resolution) + / (grid.grid_resolution[0] * grid.grid_resolution[1]) # Compute wake induced mixing factor mixing_factor[:,:,:,i] += \ diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index cc40527de..6b78e72d6 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -800,7 +800,7 @@ class VerticalAxisTurbine(BaseClass): rotor_diameter: float = field() vawt_blade_length: float = field() hub_height: float = field() - pP: float = field() # Remove if possible + pP: float = field() pT: float = field() TSR: float = field() generator_efficiency: float = field() diff --git a/tests/conftest.py b/tests/conftest.py index efa4fd13c..6cc541564 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -104,6 +104,8 @@ def print_test_values( ROTOR_DIAMETER = 126.0 TURBINE_GRID_RESOLUTION = 2 TIME_SERIES = False +IS_VERTICAL_AXIS_TURBINE = np.full(N_TURBINES, False) +VAWT_BLADE_LENGTHS = np.zeros(N_TURBINES) ## Unit test fixtures @@ -132,7 +134,9 @@ def turbine_grid_fixture(sample_inputs_fixture) -> TurbineGrid: wind_directions=np.array(WIND_DIRECTIONS), wind_speeds=np.array(WIND_SPEEDS), grid_resolution=TURBINE_GRID_RESOLUTION, - time_series=TIME_SERIES + time_series=TIME_SERIES, + is_vertical_axis_turbine=IS_VERTICAL_AXIS_TURBINE, + vawt_blade_lengths=VAWT_BLADE_LENGTHS, ) @pytest.fixture diff --git a/tests/flow_field_unit_test.py b/tests/flow_field_unit_test.py index 5b84403c7..da3514734 100644 --- a/tests/flow_field_unit_test.py +++ b/tests/flow_field_unit_test.py @@ -36,8 +36,8 @@ def test_initialize_velocity_field(flow_field_fixture, turbine_grid_fixture: Tur assert np.shape(flow_field_fixture.u_sorted)[0] == flow_field_fixture.n_wind_directions assert np.shape(flow_field_fixture.u_sorted)[1] == flow_field_fixture.n_wind_speeds assert np.shape(flow_field_fixture.u_sorted)[2] == N_TURBINES - assert np.shape(flow_field_fixture.u_sorted)[3] == turbine_grid_fixture.grid_resolution - assert np.shape(flow_field_fixture.u_sorted)[4] == turbine_grid_fixture.grid_resolution + assert np.shape(flow_field_fixture.u_sorted)[3] == turbine_grid_fixture.grid_resolution[0] + assert np.shape(flow_field_fixture.u_sorted)[4] == turbine_grid_fixture.grid_resolution[1] # Check that the wind speed profile was created correctly. By setting the shear # exponent to 1.0 above, the shear profile is a linear function of height and From fc6097c4635adc2d035ed994149d84699cd8d00d Mon Sep 17 00:00:00 2001 From: vallbog Date: Wed, 16 Aug 2023 23:38:45 +0200 Subject: [PATCH 36/46] WIP: Add tests for super_gaussian_vawt velocity model --- tests/conftest.py | 237 ++++++++++++++ .../super_gaussian_vawt_regression_test.py | 304 ++++++++++++++++++ 2 files changed, 541 insertions(+) create mode 100644 tests/reg_tests/super_gaussian_vawt_regression_test.py diff --git a/tests/conftest.py b/tests/conftest.py index 6cc541564..50de4ad19 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -107,6 +107,21 @@ def print_test_values( IS_VERTICAL_AXIS_TURBINE = np.full(N_TURBINES, False) VAWT_BLADE_LENGTHS = np.zeros(N_TURBINES) +X_COORDS_VAWT = [ + 0.0, + 5 * 26.0, + 10 * 26.0, +] +Y_COORDS_VAWT = [ + 0.0, + 0.0, + 0.0, +] +Z_COORDS_VAWT = [ + 40.0, + 40.0, + 40.0, +] ## Unit test fixtures @@ -365,12 +380,193 @@ def __init__(self): } self.turbine_floating["floating_correct_cp_ct_for_tilt"] = True + self.vertical_axis_turbine = { + "turbine_type": "t1_long_blades_vawt", + "is_vertical_axis_turbine": True, + "rotor_diameter": 26.0, + "vawt_blade_length": 48.0, + "hub_height": 40.0, + "pP": 1.88, + "pT": 1.88, + "generator_efficiency": 1.0, + "ref_density_cp_ct": 1.225, + "ref_tilt_cp_ct": 0.0, + "power_thrust_table": { + "power": [ + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + 0.33, + ], + "thrust": [ + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + 0.64, + ], + "wind_speed": [ + 0.0, + 2.0, + 2.5, + 3.0, + 3.5, + 4.0, + 4.5, + 5.0, + 5.5, + 6.0, + 6.5, + 7.0, + 7.5, + 8.0, + 8.5, + 9.0, + 9.5, + 10.0, + 10.5, + 11.0, + 11.5, + 12.0, + 12.5, + 13.0, + 13.5, + 14.0, + 14.5, + 15.0, + 15.5, + 16.0, + 16.5, + 17.0, + 17.5, + 18.0, + 18.5, + 19.0, + 19.5, + 20.0, + 20.5, + 21.0, + 21.5, + 22.0, + 22.5, + 23.0, + 23.5, + 24.0, + 24.5, + 25.0, + 25.01, + 25.02, + 50.0, + ], + }, + "TSR": 3.8, + } + self.farm = { "layout_x": X_COORDS, "layout_y": Y_COORDS, "turbine_type": [self.turbine] } + self.farm_vawt = { + "layout_x": X_COORDS_VAWT, + "layout_y": Y_COORDS_VAWT, + "turbine_type": [self.vertical_axis_turbine], + } + self.flow_field = { "wind_speeds": WIND_SPEEDS, "wind_directions": WIND_DIRECTIONS, @@ -381,6 +577,9 @@ def __init__(self): "reference_wind_height": self.turbine["hub_height"], } + self.flow_field_vawt = copy.deepcopy(self.flow_field) + self.flow_field_vawt["reference_wind_height"] = self.vertical_axis_turbine["hub_height"] + self.wake = { "model_strings": { "velocity_model": "jensen", @@ -442,6 +641,16 @@ def __init__(self): "smoothing_length_D": 2.0, "mixing_gain_velocity": 2.0 }, + "super_gaussian_vawt": { + "wake_expansion_coeff_y": 0.50, + "wake_expansion_coeff_z": 0.50, + "ay": 0.95, + "az": 4.5, + "by": 0.35, + "bz": 0.70, + "cy": 2.4, + "cz": 2.4, + }, }, "wake_turbulence_parameters": { "crespo_hernandez": { @@ -459,6 +668,15 @@ def __init__(self): "enable_transverse_velocities": False, } + self.wake_vawt = copy.deepcopy(self.wake) + model_strings_vawt = { + "velocity_model": "super_gaussian_vawt", + "deflection_model": "none", + "combination_model": "sosfs", + "turbulence_model": "none", + } + self.wake_vawt["model_strings"] = model_strings_vawt + self.floris = { "farm": self.farm, "flow_field": self.flow_field, @@ -475,3 +693,22 @@ def __init__(self): "description": "Inputs used for testing", "floris_version": "v3.0.0", } + + # Create a vertical-axis wind turbine (VAWT) configuration + # that is similar to `self.floris` above + self.floris_vawt = { + "farm": self.farm_vawt, + "flow_field": self.flow_field_vawt, + "wake": self.wake_vawt, + "solver": { + "type": "turbine_grid", + "turbine_grid_points": 3, + }, + "logging": { + "console": {"enable": True, "level": 1}, + "file": {"enable": False, "level": 1}, + }, + "name": "conftest", + "description": "Inputs used for testing", + "floris_version": "v3.0.0", + } diff --git a/tests/reg_tests/super_gaussian_vawt_regression_test.py b/tests/reg_tests/super_gaussian_vawt_regression_test.py new file mode 100644 index 000000000..4ab4cd9f8 --- /dev/null +++ b/tests/reg_tests/super_gaussian_vawt_regression_test.py @@ -0,0 +1,304 @@ +# Copyright 2022 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 numpy as np + +from floris.simulation import ( + average_velocity, + axial_induction, + Ct, + Floris, + power, + rotor_effective_velocity, +) +from tests.conftest import ( + assert_results_arrays, + N_TURBINES, + N_WIND_DIRECTIONS, + N_WIND_SPEEDS, + print_test_values, +) + + +DEBUG = False +VELOCITY_MODEL = "super_gaussian_vawt" +DEFLECTION_MODEL = "none" +TURBULENCE_MODEL = "none" + + +baseline = np.array( + [ + # 8 m/s + [ + [7.9736330, 0.7636044, 1691326.6483808, 0.2568973], + [5.1827276, 0.8807411, 441118.3637433, 0.3273306], + [4.9925898, 0.8926413, 385869.8808447, 0.3361718], + ], + # 9m/s + [ + [8.9703371, 0.7625570, 2407841.6718785, 0.2563594], + [5.8355012, 0.8438407, 650343.4078478, 0.3024150], + [5.6871296, 0.8514332, 598874.9374620, 0.3072782], + ], + # 10 m/s + [ + [9.9670412, 0.7529384, 3298067.1555604, 0.2514735], + [6.5341306, 0.8110034, 925882.5592972, 0.2826313], + [6.4005794, 0.8169593, 869713.2904634, 0.2860837], + ], + # 11 m/s + [ + [10.9637454, 0.7306256, 4363191.9880631, 0.2404936], + [7.3150380, 0.7819182, 1309551.0796815, 0.2665039], + [7.1452486, 0.7874908, 1219637.5477980, 0.2695064], + ], + ] +) + + +def test_regression_tandem(sample_inputs_fixture): + """ + Tandem turbines + """ + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["deflection_model"] = \ + DEFLECTION_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["turbulence_model"] = \ + TURBULENCE_MODEL + + floris = Floris.from_dict(sample_inputs_fixture.floris_vawt) + exit() # pytest -sk 'test_regression_tandem' super_gaussian_vawt_regression_test.py + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + n_turbines = floris.farm.n_turbines + n_wind_speeds = floris.flow_field.n_wind_speeds + n_wind_directions = floris.flow_field.n_wind_directions + + velocities = floris.flow_field.u + yaw_angles = floris.farm.yaw_angles + tilt_angles = floris.farm.tilt_angles + ref_tilt_cp_cts = ( + np.ones((n_wind_directions, n_wind_speeds, n_turbines)) + * floris.farm.ref_tilt_cp_cts + ) + test_results = np.zeros((n_wind_directions, n_wind_speeds, n_turbines, 4)) + + farm_avg_velocities = average_velocity( + velocities, + ) + farm_eff_velocities = rotor_effective_velocity( + floris.flow_field.air_density, + floris.farm.ref_density_cp_cts, + velocities, + yaw_angles, + tilt_angles, + ref_tilt_cp_cts, + floris.farm.pPs, + floris.farm.pTs, + floris.farm.turbine_fTilts, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + ) + farm_cts = Ct( + velocities, + yaw_angles, + tilt_angles, + ref_tilt_cp_cts, + floris.farm.turbine_fCts, + floris.farm.turbine_fTilts, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + ) + farm_powers = power( + floris.farm.ref_density_cp_cts, + farm_eff_velocities, + floris.farm.turbine_power_interps, + floris.farm.turbine_type_map, + ) + farm_axial_inductions = axial_induction( + velocities, + yaw_angles, + tilt_angles, + ref_tilt_cp_cts, + floris.farm.turbine_fCts, + floris.farm.turbine_fTilts, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + ) + for i in range(n_wind_directions): + for j in range(n_wind_speeds): + for k in range(n_turbines): + test_results[i, j, k, 0] = farm_avg_velocities[i, j, k] + test_results[i, j, k, 1] = farm_cts[i, j, k] + test_results[i, j, k, 2] = farm_powers[i, j, k] + test_results[i, j, k, 3] = farm_axial_inductions[i, j, k] + + if DEBUG: + print_test_values( + farm_avg_velocities, + farm_cts, + farm_powers, + farm_axial_inductions, + ) + + assert_results_arrays(test_results[0], baseline) + + +def test_regression_rotation(sample_inputs_fixture): + """ + Turbines in tandem and rotated. + The result from 270 degrees should match the results from 360 degrees. + + Wind from the West (Left) + + ^ + | + y + + 1|1 3 + | + | + | + 0|0 2 + |----------| + 0 1 x-> + + + Wind from the North (Top), rotated + + ^ + | + y + + 1|3 2 + | + | + | + 0|1 0 + |----------| + 0 1 x-> + + In 270, turbines 2 and 3 are waked. In 360, turbines 0 and 2 are waked. + The test compares turbines 2 and 3 with 0 and 2 from 270 and 360. + """ + TURBINE_DIAMETER = 126.0 + + sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.floris["wake"]["model_strings"]["turbulence_model"] = TURBULENCE_MODEL + sample_inputs_fixture.floris["farm"]["layout_x"] = [ + 0.0, + 0.0, + 5 * TURBINE_DIAMETER, + 5 * TURBINE_DIAMETER, + ] + sample_inputs_fixture.floris["farm"]["layout_y"] = [ + 0.0, + 5 * TURBINE_DIAMETER, + 0.0, + 5 * TURBINE_DIAMETER + ] + sample_inputs_fixture.floris["flow_field"]["wind_directions"] = [270.0, 360.0] + sample_inputs_fixture.floris["flow_field"]["wind_speeds"] = [8.0] + + floris = Floris.from_dict(sample_inputs_fixture.floris) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + farm_avg_velocities = average_velocity(floris.flow_field.u) + + t0_270 = farm_avg_velocities[0, 0, 0] # upstream + t1_270 = farm_avg_velocities[0, 0, 1] # upstream + t2_270 = farm_avg_velocities[0, 0, 2] # waked + t3_270 = farm_avg_velocities[0, 0, 3] # waked + + t0_360 = farm_avg_velocities[1, 0, 0] # waked + t1_360 = farm_avg_velocities[1, 0, 1] # upstream + t2_360 = farm_avg_velocities[1, 0, 2] # waked + t3_360 = farm_avg_velocities[1, 0, 3] # upstream + + assert np.allclose(t0_270, t1_360) + assert np.allclose(t1_270, t3_360) + assert np.allclose(t2_270, t0_360) + assert np.allclose(t3_270, t2_360) + + +def test_regression_small_grid_rotation(sample_inputs_fixture): + """ + Where wake models are masked based on the x-location of a turbine, numerical precision + can cause masking to fail unexpectedly. For example, in the configuration here one of + the turbines has these delta x values; + + [[4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13] + [4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13 4.54747351e-13]] + + and therefore the masking statement is False when it should be True. This causes the current + turbine to be affected by its own wake. This test requires that at least in this particular + configuration the masking correctly filters grid points. + """ + sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL + sample_inputs_fixture.floris["wake"]["model_strings"]["turbulence_model"] = TURBULENCE_MODEL + X, Y = np.meshgrid( + 6.0 * 126.0 * np.arange(0, 5, 1), + 6.0 * 126.0 * np.arange(0, 5, 1) + ) + X = X.flatten() + Y = Y.flatten() + + sample_inputs_fixture.floris["farm"]["layout_x"] = X + sample_inputs_fixture.floris["farm"]["layout_y"] = Y + + floris = Floris.from_dict(sample_inputs_fixture.floris) + floris.initialize_domain() + floris.steady_state_atmospheric_condition() + + # farm_avg_velocities = average_velocity(floris.flow_field.u) + velocities = floris.flow_field.u + yaw_angles = floris.farm.yaw_angles + tilt_angles = floris.farm.tilt_angles + ref_tilt_cp_cts = np.ones((1, 1, len(X))) * floris.farm.ref_tilt_cp_cts + + farm_eff_velocities = rotor_effective_velocity( + floris.flow_field.air_density, + floris.farm.ref_density_cp_cts, + velocities, + yaw_angles, + tilt_angles, + ref_tilt_cp_cts, + floris.farm.pPs, + floris.farm.pTs, + floris.farm.turbine_fTilts, + floris.farm.correct_cp_ct_for_tilt, + floris.farm.turbine_type_map, + ) + farm_powers = power( + floris.farm.ref_density_cp_cts, + farm_eff_velocities, + floris.farm.turbine_power_interps, + floris.farm.turbine_type_map, + ) + + # A "column" is oriented parallel to the wind direction + # Columns 1 - 4 should have the same power profile + # Column 5 is completely unwaked in this model + assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,5:10]) + assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,10:15]) + assert np.allclose(farm_powers[2,0,0:5], farm_powers[2,0,15:20]) + assert np.allclose(farm_powers[2,0,20], farm_powers[2,0,0]) + assert np.allclose(farm_powers[2,0,21], farm_powers[2,0,21:25]) From 635d93b4507aa731b7382f8392e15d89d618f63c Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 17 Aug 2023 11:34:10 +0200 Subject: [PATCH 37/46] Added tests for super_gaussian_vawt velocity model --- .../super_gaussian_vawt_regression_test.py | 67 ++++++++++--------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/tests/reg_tests/super_gaussian_vawt_regression_test.py b/tests/reg_tests/super_gaussian_vawt_regression_test.py index 4ab4cd9f8..f8433147e 100644 --- a/tests/reg_tests/super_gaussian_vawt_regression_test.py +++ b/tests/reg_tests/super_gaussian_vawt_regression_test.py @@ -41,27 +41,27 @@ [ # 8 m/s [ - [7.9736330, 0.7636044, 1691326.6483808, 0.2568973], - [5.1827276, 0.8807411, 441118.3637433, 0.3273306], - [4.9925898, 0.8926413, 385869.8808447, 0.3361718], + [7.9808916, 0.6400000, 128284.1947176, 0.2000000], + [5.7072555, 0.6400000, 47157.2775396, 0.2000000], + [5.3899835, 0.6400000, 39671.9573422, 0.2000000], ], - # 9m/s + # 9 m/s [ - [8.9703371, 0.7625570, 2407841.6718785, 0.2563594], - [5.8355012, 0.8438407, 650343.4078478, 0.3024150], - [5.6871296, 0.8514332, 598874.9374620, 0.3072782], + [8.9785030, 0.6400000, 182645.8538054, 0.2000000], + [6.4206624, 0.6400000, 66928.1744359, 0.2000000], + [6.0637314, 0.6400000, 56371.3864073, 0.2000000], ], # 10 m/s [ - [9.9670412, 0.7529384, 3298067.1555604, 0.2514735], - [6.5341306, 0.8110034, 925882.5592972, 0.2826313], - [6.4005794, 0.8169593, 869713.2904634, 0.2860837], + [9.9761145, 0.6400000, 250533.3207157, 0.2000000], + [7.1340694, 0.6400000, 91857.4256469, 0.2000000], + [6.7374793, 0.6400000, 77466.6641405, 0.2000000], ], # 11 m/s [ - [10.9637454, 0.7306256, 4363191.9880631, 0.2404936], - [7.3150380, 0.7819182, 1309551.0796815, 0.2665039], - [7.1452486, 0.7874908, 1219637.5477980, 0.2695064], + [10.9737259, 0.6400000, 333449.2621454, 0.2000000], + [7.8474763, 0.6400000, 122218.0126134, 0.2000000], + [7.4112273, 0.6400000, 102886.3004900, 0.2000000], ], ] ) @@ -78,7 +78,6 @@ def test_regression_tandem(sample_inputs_fixture): TURBULENCE_MODEL floris = Floris.from_dict(sample_inputs_fixture.floris_vawt) - exit() # pytest -sk 'test_regression_tandem' super_gaussian_vawt_regression_test.py floris.initialize_domain() floris.steady_state_atmospheric_condition() @@ -193,27 +192,29 @@ def test_regression_rotation(sample_inputs_fixture): In 270, turbines 2 and 3 are waked. In 360, turbines 0 and 2 are waked. The test compares turbines 2 and 3 with 0 and 2 from 270 and 360. """ - TURBINE_DIAMETER = 126.0 + TURBINE_DIAMETER = 26.0 - sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL - sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL - sample_inputs_fixture.floris["wake"]["model_strings"]["turbulence_model"] = TURBULENCE_MODEL - sample_inputs_fixture.floris["farm"]["layout_x"] = [ + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["deflection_model"] = \ + DEFLECTION_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["turbulence_model"] = \ + TURBULENCE_MODEL + sample_inputs_fixture.floris_vawt["farm"]["layout_x"] = [ 0.0, 0.0, 5 * TURBINE_DIAMETER, 5 * TURBINE_DIAMETER, ] - sample_inputs_fixture.floris["farm"]["layout_y"] = [ + sample_inputs_fixture.floris_vawt["farm"]["layout_y"] = [ 0.0, 5 * TURBINE_DIAMETER, 0.0, - 5 * TURBINE_DIAMETER + 5 * TURBINE_DIAMETER, ] - sample_inputs_fixture.floris["flow_field"]["wind_directions"] = [270.0, 360.0] - sample_inputs_fixture.floris["flow_field"]["wind_speeds"] = [8.0] + sample_inputs_fixture.floris_vawt["flow_field"]["wind_directions"] = [270.0, 360.0] + sample_inputs_fixture.floris_vawt["flow_field"]["wind_speeds"] = [8.0] - floris = Floris.from_dict(sample_inputs_fixture.floris) + floris = Floris.from_dict(sample_inputs_fixture.floris_vawt) floris.initialize_domain() floris.steady_state_atmospheric_condition() @@ -251,20 +252,22 @@ def test_regression_small_grid_rotation(sample_inputs_fixture): turbine to be affected by its own wake. This test requires that at least in this particular configuration the masking correctly filters grid points. """ - sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL - sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL - sample_inputs_fixture.floris["wake"]["model_strings"]["turbulence_model"] = TURBULENCE_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["deflection_model"] = \ + DEFLECTION_MODEL + sample_inputs_fixture.floris_vawt["wake"]["model_strings"]["turbulence_model"] = \ + TURBULENCE_MODEL X, Y = np.meshgrid( - 6.0 * 126.0 * np.arange(0, 5, 1), - 6.0 * 126.0 * np.arange(0, 5, 1) + 6.0 * 26.0 * np.arange(0, 5, 1), + 6.0 * 26.0 * np.arange(0, 5, 1), ) X = X.flatten() Y = Y.flatten() - sample_inputs_fixture.floris["farm"]["layout_x"] = X - sample_inputs_fixture.floris["farm"]["layout_y"] = Y + sample_inputs_fixture.floris_vawt["farm"]["layout_x"] = X + sample_inputs_fixture.floris_vawt["farm"]["layout_y"] = Y - floris = Floris.from_dict(sample_inputs_fixture.floris) + floris = Floris.from_dict(sample_inputs_fixture.floris_vawt) floris.initialize_domain() floris.steady_state_atmospheric_condition() From 306a2ef4e33e972ad5559a16cd0db8a7d8be5245 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 17 Aug 2023 13:49:35 +0200 Subject: [PATCH 38/46] Use is_vertical_axis_turbine attribute to construct turbine_map --- floris/simulation/farm.py | 22 +++++++++------------- floris/simulation/floris.py | 2 +- floris/simulation/solver.py | 10 +++++----- 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index fd6af012b..f77761bce 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -245,27 +245,23 @@ def construct_turbine_correct_cp_ct_for_tilt(self): [turb.correct_cp_ct_for_tilt for turb in self.turbine_map] ) - def construct_turbine_map(self): - self.turbine_map = [] + def construct_is_vertical_axis_turbine(self): + is_vertical_axis_turbine = [] for turb in self.turbine_definitions: try: - is_vertical_axis_turbine = turb['is_vertical_axis_turbine'] + is_vertical_axis_turbine.append(turb['is_vertical_axis_turbine']) except KeyError: - is_vertical_axis_turbine = False + is_vertical_axis_turbine.append(False) + self.is_vertical_axis_turbine = np.array(is_vertical_axis_turbine) - if is_vertical_axis_turbine: + def construct_turbine_map(self): + self.turbine_map = [] + for turb, is_vawt in zip(self.turbine_definitions, self.is_vertical_axis_turbine): + if is_vawt: self.turbine_map.append(VerticalAxisTurbine.from_dict(turb)) else: self.turbine_map.append(Turbine.from_dict(turb)) - def construct_is_vertical_axis_turbine(self): - self.is_vertical_axis_turbine = [] - for turb in self.turbine_definitions: - try: - self.is_vertical_axis_turbine.append(turb['is_vertical_axis_turbine']) - except KeyError: - self.is_vertical_axis_turbine.append(False) - def construct_turbine_fCts(self): self.turbine_fCts = { turb.turbine_type: turb.fCt_interp for turb in self.turbine_map diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index a1870625f..15cfe00a7 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -77,6 +77,7 @@ def __attrs_post_init__(self) -> None: self.check_deprecated_inputs() # Initialize farm quanitities that depend on other objects + self.farm.construct_is_vertical_axis_turbine() self.farm.construct_turbine_map() self.farm.construct_turbine_fCts() self.farm.construct_turbine_power_interps() @@ -91,7 +92,6 @@ def __attrs_post_init__(self) -> None: self.farm.construct_turbine_fTilts() self.farm.construct_turbine_correct_cp_ct_for_tilt() self.farm.construct_coordinates() - self.farm.construct_is_vertical_axis_turbine() self.farm.set_yaw_angles(self.flow_field.n_wind_directions, self.flow_field.n_wind_speeds) self.farm.set_tilt_to_ref_tilt( self.flow_field.n_wind_directions, diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index d73ef4665..5648853d7 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -265,6 +265,7 @@ def full_flow_sequential_solver( turbine_grid_farm = copy.deepcopy(farm) turbine_grid_flow_field = copy.deepcopy(flow_field) + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.construct_turbine_map() turbine_grid_farm.construct_turbine_fCts() turbine_grid_farm.construct_turbine_power_interps() @@ -279,7 +280,6 @@ def full_flow_sequential_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() - turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -670,6 +670,7 @@ def full_flow_cc_solver( turbine_grid_farm = copy.deepcopy(farm) turbine_grid_flow_field = copy.deepcopy(flow_field) + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.construct_turbine_map() turbine_grid_farm.construct_turbine_fCts() turbine_grid_farm.construct_turbine_power_interps() @@ -684,7 +685,6 @@ def full_flow_cc_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() - turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -1110,6 +1110,7 @@ def full_flow_turbopark_solver( # turbine_grid_farm = copy.deepcopy(farm) # turbine_grid_flow_field = copy.deepcopy(flow_field) + # turbine_grid_farm.construct_is_vertical_axis_turbine() # turbine_grid_farm.construct_turbine_map() # turbine_grid_farm.construct_turbine_fCts() # turbine_grid_farm.construct_turbine_power_interps() @@ -1119,7 +1120,6 @@ def full_flow_turbopark_solver( # turbine_grid_farm.construct_turbine_TSRs() # turbine_grid_farm.construc_turbine_pPs() # turbine_grid_farm.construct_coordinates() - # turbine_grid_farm.construct_is_vertical_axis_turbine() # turbine_grid = TurbineGrid( @@ -1354,6 +1354,7 @@ def full_flow_empirical_gauss_solver( turbine_grid_farm = copy.deepcopy(farm) turbine_grid_flow_field = copy.deepcopy(flow_field) + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.construct_turbine_map() turbine_grid_farm.construct_turbine_fCts() turbine_grid_farm.construct_turbine_power_interps() @@ -1368,7 +1369,6 @@ def full_flow_empirical_gauss_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() - turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( @@ -1624,6 +1624,7 @@ def full_flow_vawt_solver( turbine_grid_farm = copy.deepcopy(farm) turbine_grid_flow_field = copy.deepcopy(flow_field) + turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.construct_turbine_map() turbine_grid_farm.construct_turbine_fCts() turbine_grid_farm.construct_turbine_power_interps() @@ -1638,7 +1639,6 @@ def full_flow_vawt_solver( turbine_grid_farm.construct_turbine_fTilts() turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() turbine_grid_farm.construct_coordinates() - turbine_grid_farm.construct_is_vertical_axis_turbine() turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( From 16ba53457d1243aea3dbaeca03549423272558bc Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 17 Aug 2023 17:46:22 +0200 Subject: [PATCH 39/46] Documentation --- .../turbine_files/t1_long_blades_vawt.yaml | 5 ++- floris/simulation/farm.py | 3 +- floris/simulation/grid.py | 25 ++++++++--- floris/simulation/turbine.py | 45 +++++++++++-------- floris/tools/floris_interface.py | 2 +- 5 files changed, 51 insertions(+), 29 deletions(-) diff --git a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml index 917896254..09ff5a3a0 100644 --- a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml +++ b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml @@ -1,9 +1,10 @@ +# This is a vertical-axis wind turbine (VAWT). turbine_type: 't1_long_blades_vawt' is_vertical_axis_turbine: True generator_efficiency: 1.0 hub_height: 40.0 -# Cosine exponent for power loss due to yaw misalignment. Assume that this value -# doesn't matter for VAWTs since they cannot have yaw. +# Cosine exponent for power loss due to yaw misalignment. This value shouldn't +# matter since the concept of yaw is not applicable to VAWTs. pP: 1.88 # Cosine exponent for power loss due to tilt. Placeholder value. pT: 1.88 diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index f77761bce..658755372 100644 --- a/floris/simulation/farm.py +++ b/floris/simulation/farm.py @@ -31,6 +31,7 @@ convert_to_path, floris_array_converter, iter_validator, + NDArrayBool, NDArrayFloat, NDArrayObject, ) @@ -84,7 +85,7 @@ class Farm(BaseClass): correct_cp_ct_for_tilt: NDArrayFloat = field(init=False, default=[]) correct_cp_ct_for_tilt_sorted: NDArrayFloat = field(init=False, default=[]) turbine_fTilts: list = field(init=False, default=[]) - is_vertical_axis_turbine: list[bool] = field(init=False, default=[]) + is_vertical_axis_turbine: NDArrayBool = field(init=False, default=[]) def __attrs_post_init__(self) -> None: # Turbine definitions can be supplied in three ways: diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index c5ed6903b..2a40e60a1 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -25,6 +25,7 @@ from floris.type_dec import ( floris_array_converter, floris_float_type, + NDArrayBool, NDArrayFloat, NDArrayInt, ) @@ -152,21 +153,31 @@ def set_grid(self) -> None: @define class TurbineGrid(Grid): - """See `Grid` for more details. + """ + See `Grid` for more details. Args: turbine_coordinates (`list[Vec3]`): The series of turbine coordinate (`Vec3`) objects. reference_turbine_diameter (:py:obj:`float`): A reference turbine's rotor diameter. grid_resolution (:py:obj:`int` | :py:obj:`Iterable(int,)`): The number of points in each - direction of the square grid on the rotor plane. For example, grid_resolution=3 - creates a 3x3 grid within the rotor swept area. + direction of a yz-grid centered on the hub of each turbine. For example, for a + horizontal-axis turbine, grid_resolution=[3, 3] creates a 3x3 square grid in the rotor + plane within the rotor swept area. Alias grid_resolution=3 or grid_resolution=[3] can + be used to create the same 3x3 grid (setting the same number of points in each + direction). wind_directions (:py:obj:`NDArrayFloat`): Wind directions supplied by the user. wind_speeds (:py:obj:`NDArrayFloat`): Wind speeds supplied by the user. time_series (:py:obj:`bool`): Flag to indicate whether the supplied wind data is a time series. + is_vertical_axis_turbine (:py:obj:`NDArrayBool`): Indicates whether a turbine is a + vertical_axis_turbine (value `True`) or a traditional horizontal-axis turbine + (value `False`). + vawt_blade_lengths (:py:obj:`NDArrayFloat`): Length of the vertical turbine blades for + vertical-axis turbines. The value doesn't matter for horizontal-axis turbines, but + zero is usually used as a placeholder. """ # TODO: describe these and the differences between `sorted_indices` and `sorted_coord_indices` - is_vertical_axis_turbine: np.ndarray = field(converter=np.array) + is_vertical_axis_turbine: NDArrayBool = field(converter=np.array) vawt_blade_lengths: NDArrayFloat = field() sorted_indices: NDArrayInt = field(init=False) sorted_coord_indices: NDArrayInt = field(init=False) @@ -183,6 +194,8 @@ def set_grid(self) -> None: """ Create grid points at each turbine for each wind direction and wind speed in the simulation. This creates the underlying data structure for the calculation. + Both horizontal-axis turbines and vertical-axis turbines are supported, including a mix + of the two in the same farm. arrays have shape (n wind directions, n wind speeds, n turbines, m grid spanwise, m grid vertically) @@ -246,7 +259,7 @@ def set_grid(self) -> None: width_ratio = 0.5 height_ratio = 0.5 - # Use a square area for traditional horisontal-axis turbines and a rectangular area + # Use a square area for horizontal-axis turbines and a rectangular area # for vertical-axis turbines width = width_ratio * self.reference_turbine_diameter height = self.vawt_blade_lengths * self.is_vertical_axis_turbine @@ -266,7 +279,7 @@ def set_grid(self) -> None: # Conceptually, each turbine's grid can be seen as a meshgrid of a coordinate vector in # `cross_stream_vecs` and a coordinate vector in `vertical_vecs`. These variables store one - # vector per turbine, since the turbine diameters or vawt blade lenths may be different. + # vector per turbine, since the turbine diameters or vawt blade lengths may be different. if self.grid_resolution[0] == 1: # If the first grid resolution is 1, then let each coordinate vector in `cross_stream_vecs` # be [0.0], as np.linspace would just return the starting value of -width / 2 which would diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index 6b78e72d6..75778f0c9 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -605,6 +605,12 @@ class Turbine(BaseClass): at different wind speeds. wind_speed (:py:obj: List[float]): The wind speeds for which the power and thrust values are provided (m/s). + is_vertical_axis_turbine (:py:obj bool, optional): Placeholder attribute + for a vertical-axis wind turbine. If using such a turbine, please + instantiate a `VerticalAxisTurbine` instead of a `Turbine`. + vawt_blade_length (:py:obj float, optional): Placeholder attribute + for a vertical-axis wind turbine. If using such a turbine, please + instantiate a `VerticalAxisTurbine` instead of a `Turbine`. ngrid (*int*, optional): The square root of the number of points to use on the turbine grid. This number will be squared so that the points can be evenly distributed. @@ -627,8 +633,6 @@ class Turbine(BaseClass): power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) floating_tilt_table = field(default=None) floating_correct_cp_ct_for_tilt = field(default=None) - # Placeholder attributes for a vertical-axis wind turbine. If using such a turbine, - # please instantiate a VerticalAxisTurbine instead of a Turbine. is_vertical_axis_turbine = field(default=False) vawt_blade_length: float = field(default=0.0) @@ -767,33 +771,36 @@ def check_for_cp_ct_correct_flag_if_floating( @define class VerticalAxisTurbine(BaseClass): """ - Turbine is a class containing objects pertaining to the individual - turbines. - - Turbine is a model class representing a particular wind turbine. It - is largely a container of data and parameters, but also contains - methods to probe properties for output. + VerticalAxisTurbine is a model class representing a particular + vertixal-axis wind turbine (VAWT). It is largely a container of data + and parameters, but also contains methods to probe properties for output. Parameters: - rotor_diameter (:py:obj: float): The rotor diameter (m). - hub_height (:py:obj: float): The hub height (m). - pP (:py:obj: float): The cosine exponent relating the yaw - misalignment angle to power. - pT (:py:obj: float): The cosine exponent relating the rotor + rotor_diameter (:py:obj:`float`): The rotor diameter (m). + vawt_blade_length (:py:obj:`float`): Length of the vertical turbine blades (m). + hub_height (:py:obj:`float`): The hub height (m). + pP (:py:obj:`float`): The cosine exponent relating the yaw + misalignment angle to power. This value shouldn't matter since the concept + of yaw is not applicable to VAWTs. + pT (:py:obj:`float`): The cosine exponent relating the rotor tilt angle to power. - generator_efficiency (:py:obj: float): The generator + TSR (:py:obj:`float`): Tip speed ratio defined as the tangential speed of the + blade tip normalized by the incoming wind speed. + generator_efficiency (:py:obj:`float`): The generator efficiency factor used to scale the power production. - ref_density_cp_ct (:py:obj: float): The density at which the provided - cp and ct is defined + ref_density_cp_ct (:py:obj:`float`): The density at which the provided + cp and ct is defined. power_thrust_table (PowerThrustTable): A dictionary containing the following key-value pairs: - power (:py:obj: List[float]): The coefficient of power at + power (:py:obj:`list[float]`): The coefficient of power at different wind speeds. - thrust (:py:obj: List[float]): The coefficient of thrust + thrust (:py:obj:`list[float]`): The coefficient of thrust at different wind speeds. - wind_speed (:py:obj: List[float]): The wind speeds for + wind_speed (:py:obj:`list[float]`): The wind speeds for which the power and thrust values are provided (m/s). + is_vertical_axis_turbine (:py:obj:`bool`, optional): Value `True` indicates that + this is a vertical-axis turbine and not a traditional horizontal-axis turbine. """ turbine_type: str = field() diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index e308d8016..dfccbb9b6 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -106,7 +106,7 @@ def __init__(self, configuration: dict | str | Path): if len(np.unique(self.floris.farm.is_vertical_axis_turbine)) == 2: raise NotImplementedError( - 'There is currently no solver that supports a mix of horisontal-axis turbines ' + 'There is currently no solver that supports a mix of horizontal-axis turbines ' 'and vertical-axis turbines in the same farm.' ) From 609328214d35a0ab7402603401ddf0ad5ae22988 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 17 Aug 2023 21:37:59 +0200 Subject: [PATCH 40/46] Clarification --- floris/tools/floris_interface.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index dfccbb9b6..a8abc370c 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -86,6 +86,7 @@ def __init__(self, configuration: dict | str | Path): if np.any(np.array(ngrid) > 3): if isinstance(ngrid, int) or \ (isinstance(ngrid, Iterable) and len(np.array(ngrid)) == 1): + # An alias (a single number) was used to set the grid resolution self.logger.error( f"`turbine_grid_points` is {ngrid} (which is an alias for grid resolution " f"{self.floris.grid.grid_resolution} in the cross-stream and vertical " @@ -94,7 +95,9 @@ def __init__(self, configuration: dict | str | Path): "points reduce the computational performance but have a small change on " "accuracy." ) - else: + elif isinstance(ngrid, Iterable) and len(np.array(ngrid)) == 2: + # The grid resolution was set in the same way as it is stored in the + # `TurbineGrid` class self.logger.error( f"`turbine_grid_points` is {ngrid}. This is larger then the recommended " "value of less than or equal to 3 in each direction. High amounts of " From 38d401744aa7955ff7986fdafd6f9d3a19851b08 Mon Sep 17 00:00:00 2001 From: vallbog Date: Thu, 17 Aug 2023 22:00:13 +0200 Subject: [PATCH 41/46] Fixed comments --- examples/30_vertical_axis_wind_turbine.py | 2 +- floris/simulation/wake_velocity/super_gaussian_vawt.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py index ee7037c32..0d73eacad 100644 --- a/examples/30_vertical_axis_wind_turbine.py +++ b/examples/30_vertical_axis_wind_turbine.py @@ -30,7 +30,7 @@ """ The first part of this example shows a characteristic wake of a vertical-axis wind turbine (VAWT). It is based on case 3 in :cite:``abkar2019theoretical. The super-Gaussian velocity model with -default coefficients is used, which allows the wake to have differenti characteristics in the +default coefficients is used, which allows the wake to have different characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: rotor diameter * length of the vertical turbine blades. diff --git a/floris/simulation/wake_velocity/super_gaussian_vawt.py b/floris/simulation/wake_velocity/super_gaussian_vawt.py index 8f74d678d..d9e7ac54a 100644 --- a/floris/simulation/wake_velocity/super_gaussian_vawt.py +++ b/floris/simulation/wake_velocity/super_gaussian_vawt.py @@ -29,7 +29,7 @@ @define class SuperGaussianVAWTVelocityDeficit(BaseModel): """ - This is a super-Gaussian wake velocity model for Vertical Axis Wind Turbines (VAWTs). + This is a super-Gaussian wake velocity model for vertical-axis wind turbines (VAWTs). The model is based on :cite:`ouro2021theoretical` and allows the wake to have different characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: @@ -112,7 +112,7 @@ def function( fac = np.sqrt(1 - ct_i) beta = 0.5 * (1 + fac) / fac - # Characteristic nondimensional wake width at x - x_i = 0 + # `epsilon` is the characteristic nondimensional wake width at x - x_i = 0 ny_0 = ay + cy nz_0 = az + cz eta_0 = 1/ny_0 + 1/nz_0 @@ -123,7 +123,7 @@ def function( sigma_y_tilde = ne.evaluate("ky_star * x_tilde + epsilon") sigma_z_tilde = ne.evaluate("kz_star * x_tilde_H + epsilon") - # Eponents which determine the wake shape + # Exponents which determine the wake shape ny = ne.evaluate("ay * exp(-by * x_tilde) + cy") nz = ne.evaluate("az * exp(-bz * x_tilde_H) + cz") From 554d77c96ab99647932908495c495f07b528b39d Mon Sep 17 00:00:00 2001 From: vallbog Date: Fri, 18 Aug 2023 10:08:39 +0200 Subject: [PATCH 42/46] Typo --- examples/30_vertical_axis_wind_turbine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py index 0d73eacad..385e3c907 100644 --- a/examples/30_vertical_axis_wind_turbine.py +++ b/examples/30_vertical_axis_wind_turbine.py @@ -29,7 +29,7 @@ """ The first part of this example shows a characteristic wake of a vertical-axis wind turbine (VAWT). -It is based on case 3 in :cite:``abkar2019theoretical. The super-Gaussian velocity model with +It is based on case 3 in :cite:`abkar2019theoretical`. The super-Gaussian velocity model with default coefficients is used, which allows the wake to have different characteristics in the cross-stream (y) and vertical direction (z). The initial wake shape is closely related to the turbine cross section, which is: From 7f48520aea44fd055a1d4a99b666140c1774c215 Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 13 Nov 2023 14:34:05 +0100 Subject: [PATCH 43/46] Add flexibility to width and height ratio in TurbineGrid --- floris/simulation/grid.py | 22 +++++++++++------ .../super_gaussian_vawt_regression_test.py | 24 +++++++++---------- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 2a40e60a1..d93fd617a 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -256,15 +256,23 @@ def set_grid(self) -> None: # the rotor radius. # Defaults to 0.5. + # For horisontal axis turbines, makes sure that the square grid is within the + # rotor swept area width_ratio = 0.5 height_ratio = 0.5 - - # Use a square area for horizontal-axis turbines and a rectangular area - # for vertical-axis turbines - width = width_ratio * self.reference_turbine_diameter - height = self.vawt_blade_lengths * self.is_vertical_axis_turbine - height += self.reference_turbine_diameter * np.invert(self.is_vertical_axis_turbine) - height *= height_ratio + # For vertical-axis turbines, let the rectangular grid coincide with the turbine + # cross-section area (as seen by the incoming flow) + width_ratio_vawt = 1.0 + height_ratio_vawt = 1.0 + + is_vawt = self.is_vertical_axis_turbine + is_hawt = np.invert(is_vawt) + width_ratios = width_ratio * is_hawt + width_ratio_vawt * is_vawt + height_ratios = height_ratio * is_hawt + height_ratio_vawt * is_vawt + + width = width_ratios * self.reference_turbine_diameter + height = self.reference_turbine_diameter * is_hawt + self.vawt_blade_lengths * is_vawt + height *= height_ratios template_grid = np.ones( ( diff --git a/tests/reg_tests/super_gaussian_vawt_regression_test.py b/tests/reg_tests/super_gaussian_vawt_regression_test.py index f8433147e..32abcd47f 100644 --- a/tests/reg_tests/super_gaussian_vawt_regression_test.py +++ b/tests/reg_tests/super_gaussian_vawt_regression_test.py @@ -41,27 +41,27 @@ [ # 8 m/s [ - [7.9808916, 0.6400000, 128284.1947176, 0.2000000], - [5.7072555, 0.6400000, 47157.2775396, 0.2000000], - [5.3899835, 0.6400000, 39671.9573422, 0.2000000], + [7.9131763, 0.6400000, 125205.2867467, 0.2000000], + [6.6384807, 0.6400000, 74051.6617894, 0.2000000], + [6.3656271, 0.6400000, 65300.4197119, 0.2000000], ], # 9 m/s [ - [8.9785030, 0.6400000, 182645.8538054, 0.2000000], - [6.4206624, 0.6400000, 66928.1744359, 0.2000000], - [6.0637314, 0.6400000, 56371.3864073, 0.2000000], + [8.9023233, 0.6400000, 178230.8663179, 0.2000000], + [7.4682908, 0.6400000, 105157.0165519, 0.2000000], + [7.1613305, 0.6400000, 92942.2211551, 0.2000000], ], # 10 m/s [ - [9.9761145, 0.6400000, 250533.3207157, 0.2000000], - [7.1340694, 0.6400000, 91857.4256469, 0.2000000], - [6.7374793, 0.6400000, 77466.6641405, 0.2000000], + [9.8914704, 0.6400000, 244442.7624721, 0.2000000], + [8.2981009, 0.6400000, 144511.9190004, 0.2000000], + [7.9570339, 0.6400000, 127199.4232151, 0.2000000], ], # 11 m/s [ - [10.9737259, 0.6400000, 333449.2621454, 0.2000000], - [7.8474763, 0.6400000, 122218.0126134, 0.2000000], - [7.4112273, 0.6400000, 102886.3004900, 0.2000000], + [10.8806174, 0.6400000, 325305.2089358, 0.2000000], + [9.1279110, 0.6400000, 192175.9529676, 0.2000000], + [8.7527373, 0.6400000, 169561.6231444, 0.2000000], ], ] ) From 01817070529fc4aa6a7a04ad7b2d2587ac540c1b Mon Sep 17 00:00:00 2001 From: vallbog Date: Sun, 19 Nov 2023 17:17:35 +0100 Subject: [PATCH 44/46] Change combination model in example --- examples/inputs/super_gaussian_vawt.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/inputs/super_gaussian_vawt.yaml b/examples/inputs/super_gaussian_vawt.yaml index 0ce2f74f1..452c439d4 100644 --- a/examples/inputs/super_gaussian_vawt.yaml +++ b/examples/inputs/super_gaussian_vawt.yaml @@ -30,7 +30,7 @@ flow_field: wind_veer: 0.0 wake: model_strings: - combination_model: fls + combination_model: sosfs deflection_model: none turbulence_model: none velocity_model: super_gaussian_vawt From a6f3a0a628dace8618d7d1daa42fda5192202bb5 Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 20 Nov 2023 09:50:15 +0100 Subject: [PATCH 45/46] Increase the number of grid points for accuracy --- examples/30_vertical_axis_wind_turbine.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/30_vertical_axis_wind_turbine.py b/examples/30_vertical_axis_wind_turbine.py index 385e3c907..7529e7d5c 100644 --- a/examples/30_vertical_axis_wind_turbine.py +++ b/examples/30_vertical_axis_wind_turbine.py @@ -144,10 +144,11 @@ # Switch to a two-turbine layout. Then, calculate the velocities at each turbine in a # yz-grid centered on the hub of the turbine. The grids are based on the VAWTs' rectangular -# cross-section and have a resolution of 3x5 to make them reasonably equidistant. +# cross-section and have a resolution of 12x22 such that the turbine powers can be calculated +# with reasonable accuracy. solver_settings = { 'type': 'turbine_grid', - 'turbine_grid_points': [3, 5] + 'turbine_grid_points': [12, 22], } fi.reinitialize(layout_x=[0.0, 78.0], layout_y=[0.0, 0.0], solver_settings=solver_settings) fi.calculate_wake() From 2ad39a0239d54aa4d1c7db08e2ffcade390d0cad Mon Sep 17 00:00:00 2001 From: vallbog Date: Mon, 20 Nov 2023 14:21:06 +0100 Subject: [PATCH 46/46] Small update of power_thrust_table --- examples/inputs/turbine_files/t1_long_blades_vawt.yaml | 6 +++--- tests/conftest.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml index 09ff5a3a0..7e7e03ce4 100644 --- a/examples/inputs/turbine_files/t1_long_blades_vawt.yaml +++ b/examples/inputs/turbine_files/t1_long_blades_vawt.yaml @@ -124,6 +124,9 @@ power_thrust_table: - 0.64 wind_speed: - 0.0 + - 0.5 + - 1.0 + - 1.5 - 2.0 - 2.5 - 3.0 @@ -171,6 +174,3 @@ power_thrust_table: - 24.0 - 24.5 - 25.0 - - 25.01 - - 25.02 - - 50.0 diff --git a/tests/conftest.py b/tests/conftest.py index 50de4ad19..80f8483af 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -500,6 +500,9 @@ def __init__(self): ], "wind_speed": [ 0.0, + 0.5, + 1.0, + 1.5, 2.0, 2.5, 3.0, @@ -547,9 +550,6 @@ def __init__(self): 24.0, 24.5, 25.0, - 25.01, - 25.02, - 50.0, ], }, "TSR": 3.8,