Add infrastructure for vertical-axis wind turbines

docs/references.bib

@@ -271,3 +271,25 @@ @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},
+}
+
+@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},
+}

examples/29_plot_velocity_deficit_profiles.py

@@ -0,0 +1,157 @@
+# 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
+
+from floris.tools import FlorisInterface
+from floris.tools.visualization import VelocityProfilesFigure
+
+
+"""
+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):
+    return fi.sample_velocity_deficit_profiles(
+        direction=direction,
+        downstream_dists=downstream_dists,
+        resolution=resolution,
+        homogeneous_wind_speed=homogeneous_wind_speed,
+    )
+
+if __name__ == '__main__':
+    D = 126.0 # Turbine diameter
+    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
+
+    # 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')
+
+    # 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.
+    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', 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([-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(
+        'Velocity decifit profiles from different velocity models',
+        fontsize=14,
+    )
+
+    # Second part of example:
+    # 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]:
+        # 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)
+        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]
+            # 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
+            # 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
+
+        # 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',
+            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()

examples/30_vertical_axis_wind_turbine.py

@@ -0,0 +1,175 @@
+# 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 (
+    plot_rotor_values,
+    VelocityProfilesFigure,
+    visualize_cut_plane,
+)
+
+
+"""
+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 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.
+
+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
+        :filter: docname in docnames
+"""
+
+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)
+
+# 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 12x22 such that the turbine powers can be calculated
+# with reasonable accuracy.
+solver_settings = {
+    'type': 'turbine_grid',
+    '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()
+
+# 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()