Wind Turbine Actuator Disc¶

Optimizing an Actuator Disc Model to Find Betz Limit for Wind Turbines¶

The Betz limit is the theoretical maximum amount of kinetic energy that a wind turbine can extract from the flow. This limit was derived analytically by Albert Betz in 1919, but it can also be found numerically using an optimizer and a simple actuator disc model for a wind-turbine.

The actuator disc model of a wind turbine.

The Actuator Disc Model¶

import scipy

class ActuatorDisc(om.ExplicitComponent):
    """Simple wind turbine model based on actuator disc theory"""

    def setup(self):

        # Inputs
        self.add_input('a', 0.5, desc="Induced Velocity Factor")
        self.add_input('Area', 10.0, units="m**2", desc="Rotor disc area")
        self.add_input('rho', 1.225, units="kg/m**3", desc="air density")
        self.add_input('Vu', 10.0, units="m/s", desc="Freestream air velocity, upstream of rotor")

        # Outputs
        self.add_output('Vr', 0.0, units="m/s",
                        desc="Air velocity at rotor exit plane")
        self.add_output('Vd', 0.0, units="m/s",
                        desc="Slipstream air velocity, downstream of rotor")
        self.add_output('Ct', 0.0, desc="Thrust Coefficient")
        self.add_output('thrust', 0.0, units="N",
                        desc="Thrust produced by the rotor")
        self.add_output('Cp', 0.0, desc="Power Coefficient")
        self.add_output('power', 0.0, units="W", desc="Power produced by the rotor")
        
        # Every output depends on `a`
        self.declare_partials(of='*', wrt='a', method='cs')

        # Other dependencies
        self.declare_partials(of='Vr', wrt=['Vu'], method='cs')
        self.declare_partials(of=['thrust', 'power'], wrt=['Area', 'rho', 'Vu'], method='cs')

    def compute(self, inputs, outputs):
        """ Considering the entire rotor as a single disc that extracts
        velocity uniformly from the incoming flow and converts it to
        power."""

        a = inputs['a']
        Vu = inputs['Vu']

        qA = .5 * inputs['rho'] * inputs['Area'] * Vu ** 2

        outputs['Vd'] = Vd = Vu * (1 - 2 * a)
        outputs['Vr'] = .5 * (Vu + Vd)

        outputs['Ct'] = Ct = 4 * a * (1 - a)
        outputs['thrust'] = Ct * qA

        outputs['Cp'] = Cp = Ct * (1 - a)
        outputs['power'] = Cp * qA * Vu

Build the problem and add the actuator disc model¶

prob = om.Problem()
prob.model.add_subsystem('a_disk', ActuatorDisc(),
                         promotes_inputs=['a', 'Area', 'rho', 'Vu']);

Define the driver and optimization problem¶

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'

prob.model.add_design_var('a', lower=0., upper=1.)
prob.model.add_design_var('Area', lower=0., upper=1.)

# negative one so we maximize the objective
prob.model.add_objective('a_disk.Cp', scaler=-1)

Call setup and provide the initial guess¶

prob.setup()

prob.set_val('a', .5)
prob.set_val('Area', 10.0, units='m**2')
prob.set_val('rho', 1.225, units='kg/m**3')
prob.set_val('Vu', 10.0, units='m/s')

Run the optimization¶

fail = prob.run_driver()

Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.592592590665925
            Iterations: 5
            Function evaluations: 6
            Gradient evaluations: 5
Optimization Complete
-----------------------------------

/usr/share/miniconda/envs/test/lib/python3.8/site-packages/openmdao/core/total_jac.py:1574: DerivativesWarning:Design variables [('Area', inds=[0])] have no impact on the constraints or objective.

View the inputs and outputs¶

prob.model.list_inputs(values=True, units=True)

4 Input(s) in 'model'

varname  val           units  
-------  ------------  -------
a_disk
  a      [0.33335528]  None   
  Area   [1.]          m**2   
  rho    [1.225]       kg/m**3
  Vu     [10.]         m/s    

/usr/share/miniconda/envs/test/lib/python3.8/site-packages/openmdao/core/system.py:3425: OpenMDAOWarning:<model> <class Group>: 'value' is deprecated and will be removed in 4.0. Please index in using 'val'

[('a_disk.a',
  {'units': None, 'val': array([0.33335528]), 'value': array([0.33335528])}),
 ('a_disk.Area', {'units': 'm**2', 'val': array([1.]), 'value': array([1.])}),
 ('a_disk.rho',
  {'units': 'kg/m**3', 'val': array([1.225]), 'value': array([1.225])}),
 ('a_disk.Vu', {'units': 'm/s', 'val': array([10.]), 'value': array([10.])})]

prob.model.list_outputs(values=True, units=True)

6 Explicit Output(s) in 'model'

varname   val             units
--------  --------------  -----
a_disk
  Vr      [6.6664472]     m/s  
  Vd      [3.33289439]    m/s  
  Ct      [0.88891815]    None 
  thrust  [54.44623668]   N    
  Cp      [0.59259259]    None 
  power   [362.96296178]  W    


0 Implicit Output(s) in 'model'

/usr/share/miniconda/envs/test/lib/python3.8/site-packages/openmdao/core/system.py:3567: OpenMDAOWarning:<model> <class Group>: 'value' is deprecated and will be removed in 4.0. Please index in using 'val'

[('a_disk.Vr',
  {'val': array([6.6664472]), 'units': 'm/s', 'value': array([6.6664472])}),
 ('a_disk.Vd',
  {'val': array([3.33289439]), 'units': 'm/s', 'value': array([3.33289439])}),
 ('a_disk.Ct',
  {'val': array([0.88891815]), 'units': None, 'value': array([0.88891815])}),
 ('a_disk.thrust',
  {'val': array([54.44623668]), 'units': 'N', 'value': array([54.44623668])}),
 ('a_disk.Cp',
  {'val': array([0.59259259]), 'units': None, 'value': array([0.59259259])}),
 ('a_disk.power',
  {'val': array([362.96296178]),
   'units': 'W',
   'value': array([362.96296178])})]