Optimizing an Actuator Disk 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 disk model for a wind-turbine.

from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp, ExplicitComponent

class ActuatorDisc(ExplicitComponent):

def setup(self):

# Inputs
self.add_input('Area', 10.0, units="m**2", desc="Rotor disc area")
self.add_input('Vu', 10.0, units="m/s", desc="Freestream air velocity, upstream of rotor")

# Outputs
desc="Air velocity at rotor exit plane")
desc="Slipstream air velocity, downstream of rotor")
desc="Thrust produced by the rotor")
self.add_output('power', 0.0, units="W", desc="Power produced by the rotor")

self.declare_partials('Vr', ['a', 'Vu'])
self.declare_partials('Vd', 'a')
self.declare_partials('Ct', 'a')
self.declare_partials('thrust', ['a', 'Area', 'rho', 'Vu'])
self.declare_partials('Cp', 'a')
self.declare_partials('power', ['a', 'Area', 'rho', 'Vu'])

def compute(self, inputs, outputs):

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

def compute_partials(self, inputs, J):

a = inputs['a']
Vu = inputs['Vu']
Area = inputs['Area']
rho = inputs['rho']

# pre-compute commonly needed quantities
a_times_area = a * Area
one_minus_a = 1.0 - a
a_area_rho_vu = a_times_area * rho * Vu

J['Vr', 'a'] = -Vu
J['Vr', 'Vu'] = one_minus_a

J['Vd', 'a'] = -2.0 * Vu

J['Ct', 'a'] = 4.0 - 8.0 * a

J['thrust', 'a'] = .5 * rho * Vu**2 * Area * J['Ct', 'a']
J['thrust', 'Area'] = 2.0 * Vu**2 * a * rho * one_minus_a
J['thrust', 'rho'] = 2.0 * a_times_area * Vu ** 2 * (one_minus_a)
J['thrust', 'Vu'] = 4.0 * a_area_rho_vu * (one_minus_a)

J['Cp', 'a'] = 4.0 * a * (2.0 * a - 2.0) + 4.0 * (one_minus_a)**2

J['power', 'a'] = 2.0 * Area * Vu**3 * a * rho * (
2.0 * a - 2.0) + 2.0 * Area * Vu**3 * rho * one_minus_a ** 2
J['power', 'Area'] = 2.0 * Vu**3 * a * rho * one_minus_a ** 2
J['power', 'rho'] = 2.0 * a_times_area * Vu ** 3 * (one_minus_a)**2
J['power', 'Vu'] = 6.0 * Area * Vu**2 * a * rho * one_minus_a**2

# build the model
prob = Problem()

promotes_inputs=['a', 'Area', 'rho', 'Vu'])

# setup the optimization
prob.driver = ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'

# negative one so we maximize the objective

prob.setup()
prob.run_driver()
# minimum value
print(prob['a_disk.Cp'])
print(prob['a'])

Optimization terminated successfully.    (Exit mode 0)
Current function value: -0.592592592584
Iterations: 7
Function evaluations: 9
[ 0.33333478]