actuator_disc_sparseJ.py

from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float

import numpy as np
from itertools import combinations

class ActuatorDisc(Component):
    """Simple wind turbine model based on actuator disc theory"""

    # inputs
    a = Float(.5, iotype="in", desc="Induced Velocity Factor")
    Area = Float(10, iotype="in", desc="Rotor disc area", units="m**2", low=0)
    rho = Float(1.225, iotype="in", desc="air density", units="kg/m**3")
    Vu = Float(10, iotype="in", desc="Freestream air velocity, upstream of rotor", units="m/s")

    # outputs
    Vr = Float(iotype="out", desc="Air velocity at rotor exit plane", units="m/s")
    Vd = Float(iotype="out", desc="Slipstream air velocity, dowstream of rotor", units="m/s")
    Ct = Float(iotype="out", desc="Thrust Coefficient")
    thrust = Float(iotype="out", desc="Thrust produced by the rotor", units="N")
    Cp = Float(iotype="out", desc="Power Coefficient")
    power = Float(iotype="out", desc="Power produced by the rotor", units="W")

    def execute(self):
        # we use 'a' and 'V0' a lot, so make method local variables

        a = self.a
        Vu = self.Vu

        qA = .5*self.rho*self.Area*Vu**2
        """
        rho = .5*self.rho*self.Area*Vu**2
        area = .5*self.rho*Vu**2
        Vu = self.rho*self.Area*Vu
        a = 0
        """

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

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

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

    def provideJ(self):
        self.J = {}

        # pre-compute commonly needed quantities
        a_times_area = self.a*self.Area
        rho_times_vu = self.rho*self.Vu
        one_minus_a = 1 - self.a 
        a_area_rho_vu = a_times_area*rho_times_vu


        self.J["Vr"] = {}
        self.J["Vr"]["a"] = - self.Vu
        self.J["Vr"]["Vu"] = 1 - self.a

        self.J["Vd"] = {}
        self.J["Vd"]["a"] = -2*self.Vu
        self.J["Vd"]["Vu"] = 1 - 2*self.a

        self.J["Ct"] = {}
        self.J["Ct"]["a"] = 4 - 8*self.a

        self.J["thrust"] = {}
        self.J["thrust"]["a"] = -2.0*self.Area*self.Vu**2*self.a*self.rho + 2.0*self.Area*self.Vu**2*self.rho*one_minus_a
        self.J["thrust"]["Area"] = 2.0*self.Vu**2*self.a*self.rho*one_minus_a
        self.J["thrust"]["rho"] = 2.0*a_times_area*self.Vu**2*(one_minus_a)
        self.J["thrust"]["Vu"] = 4.0*a_area_rho_vu*(one_minus_a)

        self.J["Cp"] = {}
        self.J["Cp"]["a"] = 4*self.a*(2*self.a - 2) + 4*(one_minus_a)**2

        self.J["power"] = {}
        self.J["power"]["a"] = 2.0*self.Area*self.Vu**3*self.a*self.rho*(2*self.a - 2) + 2.0*self.Area*self.Vu**3*self.rho*one_minus_a**2
        self.J["power"]["Area"] = 2.0*self.Vu**3*self.a*self.rho*one_minus_a**2
        self.J["power"]["rho"] = 2.0*a_times_area*self.Vu**3*(one_minus_a)**2
        self.J["power"]["Vu"] = 6.0*self.Area*self.Vu**2*self.a*self.rho*one_minus_a**2

    def list_deriv_vars(self):
        input_keys = ('a', 'Area', 'rho', 'Vu')
        output_keys = ('Vr', 'Vd','Ct','thrust','Cp','power',)
        return input_keys, output_keys


    def apply_deriv(self, arg, result):
        for name in result:
            result[name] += sum([self.J[name][var]*arg[var] \
                            for var in self.J[name] if var in arg])

    def apply_derivT(self, arg, result):
        for name in result:
            result[name] += sum([self.J[var][name]*arg[var] for var in self.J \
                                 if (var in arg and name in self.J[var])])

if __name__ == "__main__":

    comp = ActuatorDisc()
    comp.run()

    comp.check_gradient(mode="adjoint")