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

import scipy
import openmdao.api as om


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.)

# 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
-----------------------------------

View the inputs and outputs

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

4 Input(s) in 'model'

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

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

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

6 Explicit Output(s) in 'model'

varname   val              units  prom_name    
--------  ---------------  -----  -------------
a_disk
  Vr      [6.6664472]      m/s    a_disk.Vr    
  Vd      [3.33289439]     m/s    a_disk.Vd    
  Ct      [0.88891815]     None   a_disk.Ct    
  thrust  [544.46236677]   N      a_disk.thrust
  Cp      [0.59259259]     None   a_disk.Cp    
  power   [3629.62961783]  W      a_disk.power 


0 Implicit Output(s) in 'model'

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

# Verify the correct outputs

# minimum value
print(prob.get_val('a_disk.Cp'))
print(prob.get_val('a'))

[0.59259259]
[0.33335528]