Curled wake model wind farm solver implementation

examples/examples_curled_wake_model/01_test_cwm_solve.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+from floris import FlorisModel
+
+
+fmodel = FlorisModel("../inputs/cwm.yaml")
+
+# Changing the wind farm layout uses FLORIS' set method to a two-turbine layout
+fmodel.set(layout_x=[100, 500.0, 1000.0], layout_y=[0.0, 0.0, 0.0])
+
+
+single_condition = False
+
+if single_condition:
+    fmodel.set(
+        wind_directions=np.array([270.0]),
+        wind_speeds=np.array([8.0]),
+        turbulence_intensities=np.array([0.06])
+    )
+else:
+    fmodel.set(
+        wind_directions=np.array([270.0, 180.0]),
+        wind_speeds=np.array([9.0, 8.0]),
+        turbulence_intensities=np.array([0.06, 0.06])
+    )
+
+fmodel.run()
+
+# There are functions to get either the power of each turbine, or the farm power
+turbine_powers = fmodel.get_turbine_powers() / 1000.0
+farm_power = fmodel.get_farm_power() / 1000.0
+
+print("Turbine powers")
+print(turbine_powers)
+print("Shape: ", turbine_powers.shape)
+
+print("Farm power")
+print(farm_power)
+print("Shape: ", farm_power.shape)

examples/examples_curled_wake_model/02_random_layout.py

@@ -0,0 +1,93 @@
+import numpy as np
+
+from floris import FlorisModel
+
+
+def main():
+    # The FlorisModel class is the entry point for most usage.
+    # Initialize using an input yaml file
+    fmodel = FlorisModel("../inputs/cwm.yaml")
+    np.random.seed(0)
+
+    # Changing the wind farm layout uses FLORIS' set method to a two-turbine layout
+    #fmodel.set(layout_x=[100, 500.0], layout_y=[300.0, 100.0])
+
+#    layout_x = [100, 100, 500, 600, 800, 1200]
+#    layout_y = [0, 100, 500, 0, 300, 400]
+
+    # Example
+    D = 126.  # Turbine diameter [m]
+
+    # Test a random layout generator
+    N= 20  # Number of turbines
+    xd = 5 * D  # Minimum distance in x-direction
+    yd = 3 * D  # Minimum distance in y-direction
+    xlim = (0, np.sqrt(N*2.5) * xd)  # X limits
+    ylim = (0, np.sqrt(N*2.5) * yd)  # Y limits
+    layout_x, layout_y = random_layout(N=N, xd=5*D, yd=3*D, xlim=xlim, ylim=ylim)
+    print(layout_x, layout_y)
+
+    yaw_angles = np.random.default_rng(0).uniform(-20, 20, size=(1, N))
+    #yaw_angles = -20 * np.ones((1, N))
+    print("Yaw angles:", yaw_angles)
+
+    fmodel.set(layout_x=layout_x, layout_y=layout_y,
+               yaw_angles=yaw_angles,
+               )
+
+
+    # Single wind condition
+    fmodel.set(
+        wind_directions=np.array([270.0]),
+        wind_speeds=np.array([8.0]),
+        turbulence_intensities=np.array([0.06])
+    )
+
+    # After the set method, the run method is called to perform the simulation
+    fmodel.run()
+
+    # There are functions to get either the power of each turbine, or the farm power
+    turbine_powers = fmodel.get_turbine_powers() / 1000.0
+    farm_power = fmodel.get_farm_power() / 1000.0
+
+    print("Turbine powers")
+    print(turbine_powers)
+    print("Shape: ", turbine_powers.shape)
+
+    print("Farm power")
+    print(farm_power)
+    print("Shape: ", farm_power.shape)
+
+
+
+def random_layout(N, xd, yd, xlim=(0, 100), ylim=(0, 100), max_attempts=100000):
+    layout_x = []
+    layout_y = []
+    attempts = 0
+
+    while len(layout_x) < N and attempts < max_attempts:
+        x_new = np.random.uniform(*xlim)
+        y_new = np.random.uniform(*ylim)
+
+        # Check distances with existing points
+        if layout_x:
+            dx = np.abs(np.array(layout_x) - x_new)
+            dy = np.abs(np.array(layout_y) - y_new)
+
+            # Valid if no other point is too close in BOTH directions
+            if np.any((dx < xd) & (dy < yd)):
+                attempts += 1
+                continue
+
+        layout_x.append(x_new)
+        layout_y.append(y_new)
+        attempts = 0  # reset since we succeeded
+
+    if len(layout_x) < N:
+        raise RuntimeError(f"Could not place {N} points after {max_attempts} attempts.")
+
+    return np.array(layout_x), np.array(layout_y)
+
+
+if __name__ == "__main__":
+    main()

examples/inputs/cwm.yaml

@@ -0,0 +1,66 @@
+
+name: CWM
+description: Three turbines using Curled Wake Model
+floris_version: v4
+
+logging:
+  console:
+    enable: true
+    level: WARNING
+  file:
+    enable: false
+    level: WARNING
+
+# TODO: Likely will need to overwrite this.
+solver:
+  type: turbine_grid
+  turbine_grid_points: 3
+
+farm:
+  layout_x:
+  - 0.0
+  - 630.0
+  - 1260.0
+  layout_y:
+  - 0.0
+  - 0.0
+  - 0.0
+  turbine_type:
+  - nrel_5MW
+
+flow_field:
+  air_density: 1.225
+  reference_wind_height: -1 # -1 is code for use the hub height
+  turbulence_intensities:
+  - 0.06
+  wind_directions:
+  - 270.0
+  wind_shear: 0.12
+  wind_speeds:
+  - 8.0
+  wind_veer: 10.0
+
+wake:
+  model_strings:
+    combination_model: none
+    deflection_model: none
+    turbulence_model: none
+    velocity_model: cwm
+
+  enable_secondary_steering: false
+  enable_yaw_added_recovery: false
+  enable_transverse_velocities: false
+  enable_active_wake_mixing: false
+
+  wake_deflection_parameters:
+    None: {}
+
+  wake_velocity_parameters:
+    cwm:
+      C: 4.0
+      dx_D: 0.08
+      dy_D: 0.1
+      dz_D: 0.1
+
+  wake_turbulence_parameters:
+    None: {}

floris/core/__init__.py

@@ -47,6 +47,7 @@
 from .wake import WakeModelManager
 from .solver import (
     cc_solver,
+    curled_wake_solver,
     empirical_gauss_solver,
     full_flow_cc_solver,
     full_flow_empirical_gauss_solver,

floris/core/core.py

@@ -12,6 +12,7 @@
 from floris.core import (
     BaseClass,
     cc_solver,
+    curled_wake_solver,
     empirical_gauss_solver,
     Farm,
     FlowField,
@@ -198,6 +199,13 @@ def steady_state_atmospheric_condition(self):
                 self.grid,
                 self.wake
             )
+        elif vel_model=="cwm":
+            curled_wake_solver(
+                self.farm,
+                self.flow_field,
+                self.grid,
+                self.wake
+            )
         else:
             sequential_solver(
                 self.farm,

floris/core/solver.py

@@ -24,9 +24,15 @@
     yaw_added_turbulence_mixing,
 )
 from floris.core.wake_velocity.empirical_gauss import awc_added_wake_mixing
-from floris.type_dec import NDArrayFloat
+from floris.type_dec import (
+    floris_float_type,
+    NDArrayFloat,
+)
 from floris.utilities import cosd
 
+# For the CWM
+from .solver_cwm import curled_wake_solver
+
 
 def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid):
     """