Wind Turbine Actuator Disc
Contents
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.59259259]
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
------- ------------ -------
a_disk
a [0.33335528] None
Area [10.] m**2
rho [1.225] kg/m**3
Vu [10.] m/s
[('a_disk.a', {'units': None, 'val': array([0.33335528])}),
('a_disk.Area', {'units': 'm**2', 'val': array([10.])}),
('a_disk.rho', {'units': 'kg/m**3', 'val': array([1.225])}),
('a_disk.Vu', {'units': 'm/s', 'val': array([10.])})]
prob.model.list_outputs(val=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 [544.46236677] N
Cp [0.59259259] None
power [3629.62961783] W
0 Implicit Output(s) in 'model'
[('a_disk.Vr', {'units': 'm/s', 'val': array([6.6664472])}),
('a_disk.Vd', {'units': 'm/s', 'val': array([3.33289439])}),
('a_disk.Ct', {'units': None, 'val': array([0.88891815])}),
('a_disk.thrust', {'units': 'N', 'val': array([544.46236677])}),
('a_disk.Cp', {'units': None, 'val': array([0.59259259])}),
('a_disk.power', {'units': 'W', 'val': array([3629.62961783])})]
# Verify the correct outputs
# minimum value
print(prob.get_val('a_disk.Cp'))
print(prob.get_val('a'))
[0.59259259]
[0.33335528]