Optimizing an Actuator Disk Model to Find Betz Limit for Wind Turbines

Preliminary note:

This case is taken from OpenMDAO. OpenMDAO is an open-source computing platform for system analysis and multidisciplinary optimization developed by the NASA. Its philosophy shares some of the goals of CoSApp. Thus, this example is also an opportunity to compare both libraries.

OpenMDAO is licensed under Apache License.

Case description

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.

Creating the elements

from cosapp.base import System

class ActuatorDisc(System):
    """Simple wind turbine model based on actuator disc theory
    """
    def setup(self):
        # Inputs
        self.add_inward('a', 0.5, desc="Induced Velocity Factor")
        self.add_inward('area', 10.0, unit="m**2", desc="Rotor disc area")
        self.add_inward('rho', 1.225, unit="kg/m**3", desc="air density")
        self.add_inward('Vu', 10.0, unit="m/s",
            desc="Freestream air velocity, upstream of rotor")
        # Outputs
        self.add_outward('Vr', 0.0, unit="m/s",
            desc="Air velocity at rotor exit plane")
        self.add_outward('Vd', 0.0, unit="m/s",
            desc="Slipstream air velocity, downstream of rotor")
        self.add_outward('Ct', 0.0, desc="Thrust Coefficient")
        self.add_outward('thrust', 0.0, unit="N",
            desc="Thrust produced by the rotor")
        self.add_outward('Cp', 0.0, desc="Power Coefficient")
        self.add_outward('power', 0.0, unit="W",
            desc="Power produced by the rotor")

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

        qA = 0.5 * self.rho * self.area * Vu**2

        self.Vd = Vd = Vu * (1 - 2 * a)
        self.Vr = 0.5 * (Vu + Vd)

        self.Ct = Ct = 4 * a * (1 - a)
        self.thrust = Ct * qA

        self.Cp = Cp = Ct * (1 - a)
        self.power = Cp * qA * Vu

Setting the system

# Build model
disc = ActuatorDisc('disc')

disc.a = 1.0
disc.area = 0.1
disc.rho = 1.225
disc.Vu = 10.0

Solving the optimization problem

from cosapp.drivers import Optimizer

# Setup optimization problem
optim = disc.add_driver(Optimizer('optim', method='SLSQP'))

min_area = 1
max_area = 10

optim.add_unknown('a', lower_bound=0, upper_bound=1)
optim.add_unknown('area', lower_bound=min_area, upper_bound=max_area)
optim.set_maximum('power')

# Initialize unknowns
disc.a = 0.1
disc.area = 1.0

disc.run_drivers()

Validation

Comparison of numerical solution vs. exact optimal values:

def rel_error(actual, expected):
    return abs(actual / expected - 1)

optimum = {
    'a': 1 / 3,
    'Cp': 16 / 27,
    'area': 10,
}

solution = dict(disc.inwards.items())
solution.update(disc.outwards.items())

for varname, exact in optimum.items():
    actual = disc[varname]
    print(f"disc.{varname} = {actual}\t(error = {rel_error(actual, exact)})")

disc.a = 0.3333333286070798     (error = 1.4178760499028442e-08)
disc.Cp = 0.5925925925925923    (error = 3.3306690738754696e-16)
disc.area = 10.0        (error = 0.0)

Plot solution

import numpy as np
import itertools, copy

axes = {
    'a': np.linspace(0, 1, 21),
    'area': np.linspace(min_area, max_area, 19),
}
x = axes['a']
y = axes['area']
z = []

# Create new system to avoid side effects
s = ActuatorDisc('s')
for varname, value in disc.inwards.items():
    s[varname] = copy.copy(value)

for s.a, s.area in itertools.product(x, y):
    s.run_once()
    z.append(s.power)

z = np.reshape(z, (x.size, -1))

# Plot results
import plotly.graph_objects as go

template = '<br>'.join([
    'a: %{x}',
    'area: %{y}',
    'power: %{z:.1f}',
])

traces = [
    go.Surface(
        x = axes['a'],
        y = axes['area'],
        z = z.T,
        name = 'power',
        hovertemplate = template + '<extra></extra>',
        colorscale = 'balance',
    ),
    go.Scatter3d(
        x = [solution['a']],
        y = [solution['area']],
        z = [solution['power']],
        name = 'solution',
        marker=dict(size=6, opacity=1, color='red'),
        hovertemplate = template,
    ),
]

fig = go.Figure(data=traces)

fig.update_layout(
    scene = {
        "xaxis": {"title": "a"},
        "yaxis": {"title": "area"},
        "zaxis": {"title": "power"},
    },
    height = 600,
)

fig.show()