Add swept (projected) area option for turbines

- Concentrates the thrust into a smaller (or larger) area
- Swept diameter defaults to normal diameter

thetis/options.py

@@ -462,6 +462,9 @@ class TidalTurbineOptions(FrozenHasTraits):
     name = 'Tidal turbine options'
     diameter = PositiveFloat(
         18., help='Turbine diameter').tag(config=True)
+    swept_diameter = PositiveFloat(
+        None, allow_none=True, help='Projected diameter of the turbine '
+                                    '(defaults to diameter if not provided)').tag(config=True)
     C_support = NonNegativeFloat(
         0., help='Thrust coefficient for support structure').tag(config=True)
     A_support = NonNegativeFloat(

thetis/turbines.py

@@ -14,22 +14,31 @@
 class TidalTurbine:
     def __init__(self, options, upwind_correction=False):
         self.diameter = options.diameter
+        self.swept_diameter = options.swept_diameter or self.diameter
         self.C_support = options.C_support
         self.A_support = options.A_support
         self.upwind_correction = upwind_correction
 
+        # Check that the parameter choices make some sense - these won't break
+        # the simulation but may give unexpected results if the choice isn't
+        # understood.
+        if self.swept_diameter != self.diameter:
+            log(INFO, 'Warning - swept_diameter and plan_diameter are not equal. '
+                      'Upwind correction may not be accurate.')
+
     def _thrust_area(self, uv):
         C_T = self.thrust_coefficient(uv)
-        A_T = pi * self.diameter**2 / 4
-        fric = C_T * A_T
+        A_T = pi * self.swept_diameter**2 / 4
+        swept_area_factor = self.diameter**2/self.swept_diameter**2
+        fric = C_T * A_T * swept_area_factor
         if self.C_support:
             fric += self.C_support * self.A_support
         return fric
 
     def velocity_correction(self, uv, depth):
         fric = self._thrust_area(uv)
         if self.upwind_correction:
-            return 0.5*(1+sqrt(1-fric/(self.diameter*depth)))
+            return 0.5*(1+sqrt(1-fric/(self.swept_diameter*depth)))
         else:
             return 1
 
@@ -41,7 +50,7 @@ def friction_coefficient(self, uv, depth):
     def power(self, uv, depth):
         # ratio of discrete to upstream velocity (NOTE: should include support drag!)
         alpha = self.velocity_correction(uv, depth)
-        A_T = pi * self.diameter**2 / 4
+        A_T = pi * self.diameter**2 / 4  # power is based on true turbine diameter
         uv3 = dot(uv, uv)**1.5 / alpha**3  # upwind cubed velocity
         C_P = self.power_coefficient(uv3**(1/3))
         # this assumes the velocity through the turbine does not change due to the support (is this correct?)
@@ -156,7 +165,7 @@ def add_turbines(self, coordinates):
         """
         x = SpatialCoordinate(self.mesh)
 
-        radius = self.turbine.diameter * 0.5
+        radius = self.turbine.swept_diameter * 0.5
         for coord in coordinates:
             dx0 = (x[0] - coord[0])/radius
             dx1 = (x[1] - coord[1])/radius