From c12755c559c9778c52d7cdb2b6c38249aaa0a9d3 Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 8 Jan 2025 22:02:00 -0700 Subject: [PATCH 1/3] Add load heuristic function --- floris/floris_model.py | 107 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 1 deletion(-) diff --git a/floris/floris_model.py b/floris/floris_model.py index f521b5faf..a73c0175f 100644 --- a/floris/floris_model.py +++ b/floris/floris_model.py @@ -36,6 +36,7 @@ nested_get, nested_set, print_nested_dict, + rotate_coordinates_rel_west, ) from floris.wind_data import ( TimeSeries, @@ -503,8 +504,112 @@ def run_no_wake(self) -> None: self.core.finalize() - ### Methods for extracting turbine performance after running + #TODO: This could probably be moved elsewhere + #TODO: load_h is a stand-in term for "load heuristic" until something better comes to mind + def _get_turbine_load_h(self) -> NDArrayFloat: + + #TODO: Probably these should be function input parameters + wake_slope = 0.3 + max_dist_D = 10.0 + weight_ws = 1.0 + weight_ct = 1.0 + weight_ti = 1.0 + + # This one I'm not sure, maybe it will be ambient ti, maybe a 1D TD load_ambient_ti + load_ambient_ti = 0.10 + + # Store for convenience + n_turbines = self.n_turbines + n_findex = self.n_findex + + # TODO for now assume the first one + D = self.core.farm.rotor_diameters[0,0] + + # Get the wind direction rotated farm layouts with wind flowing left to right + coord_rotate = rotate_coordinates_rel_west(self.wind_directions, + self.core.grid.turbine_coordinates) + x_rot = coord_rotate[0] + y_rot = coord_rotate[1] + + # Get all turbine thrust coefficients + cts = self.get_turbine_thrust_coefficients() + + # Get ws, was originally thinking to use the actual rotor wind speeds (commented first line) + # but think more consistent with load methods to use ambient here since those ambient flows + # will surround the wake somehow for load inducing purposes + # ws = np.mean(self.core.flow_field.u, axis=(2,3)) + ws = self.wind_speeds + # Make the same for every turbine + ws = np.tile(ws[:, np.newaxis], (1, n_turbines)) + + # Initialize the load_ti to the load_ambient_ti + load_ti = np.ones((n_findex, n_turbines)) * load_ambient_ti + + # For each findex and turbine, compute the "load turbulence" + for findex in range(n_findex): + + # t_i is the index of the "current" turbine whose "load_ti" we want to compute + for t_i in range(n_turbines): + nearest_t_dist = 1 + np.max(x_rot) - np.min(x_rot) # Initialize to the max + 1 + + # t_j is the index of the turbine who is proposed to be waking t_i + for t_j in range(n_turbines): + + # If the same turbine skip + if t_i == t_j: + continue + + # Compute the distance from proposed waking to turbine to the current turbine + # with the distance positive if the current turbine is downstream + x_diff = x_rot[findex, t_i] - x_rot[findex, t_j] + + # If the proposed waking turbine is downstream skip + if x_diff<=0: + continue + + # If the distance between the turbines is more than the max skip + if x_diff >= max_dist_D * D: + continue + # If this is not closer than previouly found wakes skip + if x_diff >= nearest_t_dist: + continue + + # Compute the cross-stream distance + y_diff_abs = np.abs(y_rot[findex, t_i] - y_rot[findex, t_j]) + + # If the t_i turbine is above the wake cone of the + # proposed waking turbine skip + if y_diff_abs > D + x_diff * wake_slope: + continue + + # Set nearest_t_dist to the current value + nearest_t_dist = x_diff + + # Get the c_t of the proposed waking turbine + c_t_j = cts[findex, t_j] + + # Get the ambient ws (note all turbines same value + # but might as well use the current turbine) + ambient_ws = ws[findex,t_i] + + # Compute the added turbulence + ti_add = ambient_ws / (1.5 + 0.8 * (x_diff/D) / np.sqrt(c_t_j) ) + + # Update the ti_term + load_ti[findex, t_i] = ( + + np.sqrt(ti_add**2 + (load_ambient_ti * ambient_ws)**2) / ambient_ws + + ) + + + + # Compute load_h + return weight_ws * ws + weight_ct * cts + weight_ti * load_ti + + + ### Methods for extracting turbine performance after running def _get_turbine_powers(self) -> NDArrayFloat: """Calculates the power at each turbine in the wind farm. From e9ba67d1bb366b64e44e5b2531bd7bfc8ceb2f63 Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 8 Jan 2025 22:02:13 -0700 Subject: [PATCH 2/3] Add simple example --- .../001b_derating_for_load.py | 90 +++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 examples/examples_control_types/001b_derating_for_load.py diff --git a/examples/examples_control_types/001b_derating_for_load.py b/examples/examples_control_types/001b_derating_for_load.py new file mode 100644 index 000000000..df2d4d81d --- /dev/null +++ b/examples/examples_control_types/001b_derating_for_load.py @@ -0,0 +1,90 @@ +"""Example of de-rating for load + +TBD. +""" + +import matplotlib.pyplot as plt +import numpy as np + +from floris import FlorisModel + + +fmodel = FlorisModel("../inputs/gch.yaml") + +# Change to the simple-derating model turbine +# (Note this could also be done with the mixed model) +fmodel.set_operation_model("simple-derating") + +# Convert to a simple two turbine layout with derating turbines +fmodel.set(layout_x=[0, 126.0 * 5], layout_y=[0.0, 0.0]) + + +# Set the wind directions and speeds to be constant over n_findex = N time steps +N = 50 +fmodel.set( + wind_directions=270 * np.ones(N), + wind_speeds=7.0 * np.ones(N), + turbulence_intensities=0.06 * np.ones(N), +) +fmodel.run() +turbine_powers_orig = fmodel.get_turbine_powers() + +# Get the nominal operating power of the upstream turbine +nom_pow = turbine_powers_orig[0,0] + +# # Derate the front turbine from full power to 75 % +power_setpoints_front = np.linspace(nom_pow,nom_pow *.75,N) +full_rating = np.ones_like(power_setpoints_front) * 5E6 + +# Only derate the front turbine +power_setpoints = np.column_stack([power_setpoints_front, full_rating]) +fmodel.set(power_setpoints=power_setpoints) +fmodel.run() +turbine_powers_derated = fmodel.get_turbine_powers() + +# Compute the mean power +power_mean = np.mean(turbine_powers_derated,axis=1) + +# Define de-rating as percent reduction +de_rating = 100 * (nom_pow - power_setpoints_front) / nom_pow + +# Grab the load heuristic +load_h = fmodel._get_turbine_load_h() + +# Compute mean load +load_h_mean = np.mean(load_h,axis=1) + +# Plot the results +fig, axarr = plt.subplots(2,3,sharex=True, sharey='row',figsize=(12,8)) + +ax = axarr[0,0] +ax.plot(de_rating, turbine_powers_derated[:,0]/turbine_powers_derated[0,0],'k') +ax.set_title('Turbine 0') +ax.set_ylabel('Power (ratio to nomimal)') + +ax = axarr[0,1] +ax.plot(de_rating, turbine_powers_derated[:,1]/turbine_powers_derated[0,1],'k') +ax.set_title('Turbine 1') + +ax = axarr[0,2] +ax.plot(de_rating, power_mean/power_mean[0],'k') +ax.set_title('Mean') + +ax = axarr[1,0] +ax.plot(de_rating, load_h[:,0]/load_h[0,0],'k') +ax.set_ylabel('Load H (ratio to nomimal)') +ax.set_xlabel('Front De-Rating %') + +ax = axarr[1,1] +ax.plot(de_rating, load_h[:,1]/load_h[0,1],'k') +ax.set_xlabel('Front De-Rating %') + +ax = axarr[1,2] +ax.plot(de_rating, load_h_mean/load_h_mean[0],'k') +ax.set_xlabel('Front De-Rating %') + +for ax in axarr.flatten(): + ax.grid(True) + + +plt.show() From 6127285757ff9ccba6d76d7ca0773e73e15b2a59 Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 8 Jan 2025 22:02:22 -0700 Subject: [PATCH 3/3] Add tests --- tests/floris_model_integration_test.py | 31 ++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/tests/floris_model_integration_test.py b/tests/floris_model_integration_test.py index eabe764db..16797f6f3 100644 --- a/tests/floris_model_integration_test.py +++ b/tests/floris_model_integration_test.py @@ -193,6 +193,37 @@ def test_run_no_wake(): power_no_wake = fmodel.get_turbine_powers() assert len(np.unique(power_no_wake)) == 1 +def test_get_turbine_load_h(): + + fmodel = FlorisModel(configuration=YAML_INPUT) + + wind_speeds = np.array([8.0, 8.0, 8.0]) + wind_directions = np.array([270.0, 270.0, 280.0]) + turbulence_intensities = np.array([0.06, 0.06, 0.06]) + n_findex = len(wind_directions) + + layout_x = np.array([0, 1000]) + layout_y = np.array([0, 1]) + n_turbines = len(layout_x) + + fmodel.set( + wind_speeds=wind_speeds, + wind_directions=wind_directions, + turbulence_intensities=turbulence_intensities, + layout_x=layout_x, + layout_y=layout_y, + ) + + fmodel.run() + + load_h = fmodel._get_turbine_load_h() + + assert load_h.shape[0] == n_findex + assert load_h.shape[1] == n_turbines + assert load_h[0, 0] == load_h[1, 0] + + + def test_get_turbine_powers(): # Get turbine powers should return n_findex x n_turbine powers # Apply the same wind speed and direction multiple times and confirm all equal