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

result = 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]