Bug fix: Add variable wind shear to GCH calculations (#679)

NREL · Jul 10, 2023 · 8908ab4 · 8908ab4
1 parent e6906fe
commit 8908ab4
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 5 deletions.
diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py
@@ -149,6 +149,7 @@ def sequential_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
             effective_yaw_i += added_yaw
 
@@ -178,6 +179,7 @@ def sequential_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
 
         if model_manager.enable_yaw_added_recovery:
@@ -377,6 +379,7 @@ def full_flow_sequential_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
             effective_yaw_i += added_yaw
 
@@ -406,6 +409,7 @@ def full_flow_sequential_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
 
         # NOTE: exponential
@@ -554,6 +558,7 @@ def cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
                 scale=2.0,
             )
             effective_yaw_i += added_yaw
@@ -584,6 +589,7 @@ def cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
                 scale=2.0,
             )
 
@@ -784,6 +790,7 @@ def full_flow_cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
                 scale=2.0,
             )
             effective_yaw_i += added_yaw
@@ -814,6 +821,7 @@ def full_flow_cc_solver(
                 turb_Cts[:, :, i:i+1],
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
                 scale=2.0,
             )
 
@@ -948,6 +956,7 @@ def turbopark_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
             effective_yaw_i += added_yaw
 
@@ -1009,6 +1018,7 @@ def turbopark_solver(
                 ct_i,
                 TSR_i,
                 axial_induction_i,
+                flow_field.wind_shear,
             )
 
         if model_manager.enable_yaw_added_recovery:

diff --git a/floris/simulation/wake_deflection/gauss.py b/floris/simulation/wake_deflection/gauss.py
@@ -257,6 +257,7 @@ def wake_added_yaw(
     ct_i,
     tip_speed_ratio,
     axial_induction_i,
+    wind_shear,
     scale=1.0,
 ):
     """
@@ -284,7 +285,7 @@ def wake_added_yaw(
     eps_gain = 0.2
     eps = eps_gain * D  # Use set value
 
-    vel_top = ((HH + D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
+    vel_top = ((HH + D / 2) / HH) ** wind_shear * np.ones((1, 1, 1, 1, 1))
     Gamma_top = gamma(
         D,
         vel_top,
@@ -293,7 +294,7 @@ def wake_added_yaw(
         scale,
     )
 
-    vel_bottom = ((HH - D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
+    vel_bottom = ((HH - D / 2) / HH) ** wind_shear * np.ones((1, 1, 1, 1, 1))
     Gamma_bottom = -1 * gamma(
         D,
         vel_bottom,
@@ -359,6 +360,7 @@ def calculate_transverse_velocity(
     ct_i,
     tsr_i,
     axial_induction_i,
+    wind_shear,
     scale=1.0,
 ):
     """
@@ -379,8 +381,7 @@ def calculate_transverse_velocity(
     eps_gain = 0.2
     eps = eps_gain * D  # Use set value
 
-    # TODO: wind sheer is hard-coded here but should be connected to the input
-    vel_top = ((HH + D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
+    vel_top = ((HH + D / 2) / HH) ** wind_shear * np.ones((1, 1, 1, 1, 1))
     Gamma_top = sind(yaw) * cosd(yaw) * gamma(
         D,
         vel_top,
@@ -389,7 +390,7 @@ def calculate_transverse_velocity(
         scale,
     )
 
-    vel_bottom = ((HH - D / 2) / HH) ** 0.12 * np.ones((1, 1, 1, 1, 1))
+    vel_bottom = ((HH - D / 2) / HH) ** wind_shear * np.ones((1, 1, 1, 1, 1))
     Gamma_bottom = -1 * sind(yaw) * cosd(yaw) * gamma(
         D,
         vel_bottom,