Coupling RotorSE with FLORIS resulting in messy code - how to improve?

[pgebraad]

I have a working code of a coupling of RotorSE/CCBlade with my FLORIS wake mode. I have tried to use some FUSED-wind standards. It is illustrated in attached figure. The code is below.

RotorSE/CCBlade computes the axial induction and C_P of each rotor, and feeds those factors into the FLORIS wake model, which calculates effective wind speeds and powers at each turbine. Depending on thpse effective wind speeds, the C_P and axial induction may change -- this feedback of wind speeds is not implemented yet however. Now, it works for below-rated only (constant TSR).

The code is very 'awkward' and messy I think. Mainly because of three reasons:

• The assembly is built 'inline', and not as a general class. This is because RotorSE is a single-turbine model, while FLORIS uses the wt_layout structure with several turbines. Depending on the number of turbines in the wind farm, a different number of CCBlade models should be added to the assembly. I am not sure how to do this in a more generalized way with OpenMDAO.
• It still contains separate definitions of the turbine for the FLORIS wake model and for RotorSE model. Should this always be the case (different parts of the model need different information), or should there be some unified turbine definition that is connected to each of the models? Is there a Fused-Wind 'philosophy' for this?
• RotorSE does not follow the same naming standards as Fused-WIND, I will work to resolve this.

If you could give me some feedback on how I could create a better implementation of this idea, that would be great!

The code consists of three files: ExampleCall_florisRotorSE.py - Builds and runs an assembly of CCBladed models (one for each turbine), defines all the turbine properties, and runs the FLORIS wake model. floris.py - contains the FLORIS wake model florisRotorSE.py - contains some building blocks needed for the coupling

============== ExampleCall_florisRotorSE.py =================================

Builds and runs an assembly of CCBladed models and FLORIS wake model

from openmdao.main.api import Assembly from floris import FLORIS from florisRotorSE import AeroelasticHAWTVT_CCBlade, CCBladeAeroelasticHAWTVT, CCBladeCoefficients, AeroelasticHAWTVT_floris import os import numpy as np

Define flow properties

freestream_wind_speed = 8.0 # m/s air_density = 1.1716 # kg/m^3 wind_direction = 30 # deg viscosity = 1.81206e-5

define turbine positions

turbineX = np.array([1164.7, 947.2, 1682.4, 1464.9, 1982.6, 2200.1]) turbineY = np.array([1024.7, 1335.3, 1387.2, 1697.8, 2060.3, 1749.7])

create wind plant consisting of AeroelasticHAWTVT_floris turbines

florisWindPlant = FLORIS() for turbI in range(0,turbineX.size): wt = AeroelasticHAWTVT_floris()

wt.name = 'turbine%s' % turbI
wt.position = np.array([turbineX[turbI], turbineY[turbI]])

# AeroelasticHAWTVT inputs
wt.turbine_name = 'NREL 5-MW baseline turbine'

wt.rotor_diameter = 126.4
wt.rotor_area = 0.25*np.pi*wt.rotor_diameter*wt.rotor_diameter
wt.yaw = 0.0

florisWindPlant.wt_layout.add_wt(wt)

create wind plant consisting of AeroelasticHAWTVT_floris turbines

FLORISplusCCBlade = Assembly() FLORISplusCCBlade.add('floris',florisWindPlant)

add CCBlade model for each turbine, and define workflow (first all CCBlade runs, then FLORIS)

workflow = [] for turbineI in range(0,turbineX.size): ccbladeName = 'CCblade%s' % turbineI FLORISplusCCBlade.add(ccbladeName,CCBladeCoefficients()) workflow.append(ccbladeName) workflow.append('floris') print workflow

FLORISplusCCBlade.driver.workflow.add(workflow)

for turbineI in range(0,turbineX.size):

# connect CCblade-predicted axial induction and CP to FLORIS
florisTurbine = "floris.wt_layout.turbine%s" % turbineI
FLORISplusCCBlade.connect("CCblade%s.CP" % turbineI, "%s.CP" % florisTurbine)
FLORISplusCCBlade.connect("CCblade%s.axial_induction" % turbI, "%s.axial_induction" % florisTurbine)

# define rotor properties for CCBlade
ccbladeName = 'CCblade%s' % turbineI
ccblade = getattr(FLORISplusCCBlade,ccbladeName)
wt = ccblade.turbine

wt.name = 'turbine%s' % turbI
wt.position = np.array([turbineX[turbI], turbineY[turbI]])
wt.turbine_name = 'NREL 5-MW baseline turbine'
wt.rotor_diameter = 126.4
wt.rotor_area = 0.25*np.pi*wt.rotor_diameter*wt.rotor_diameter # TODO: remove double definitions
wt.tip_radius = 63.2 # TODO: remove double definitions
wt.hub_radius = 1.5
wt.nblades = 3  # number of blades
wt.hub_height = 90.0
wt.tilt_angle = 5.0
wt.cone_angle = 2.5
# CCblade inputs
wt.nSector = 8
wt.radial_locations = np.array([2.8667, 5.6000, 8.3333, 11.7500, 15.8500, 19.9500, 24.0500, 28.1500,
                                32.2500, 36.3500, 40.4500, 44.5500, 48.6500, 52.7500, 56.1667, 58.9000, 61.6333])
wt.chord = np.array([3.542, 3.854, 4.167, 4.557, 4.652, 4.458, 4.249, 4.007, 3.748,
                     3.502, 3.256, 3.010, 2.764, 2.518, 2.313, 2.086, 1.419])
wt.theta = np.array([13.308, 13.308, 13.308, 13.308, 11.480, 10.162, 9.011, 7.795,
                     6.544, 5.361, 4.188, 3.125, 2.319, 1.526, 0.863, 0.370, 0.106])
# define airfoil
basePath = os.path.join(os.getcwd(), '5MW_AFFiles')
airfoil_types = ['Cylinder1.dat','Cylinder2.dat','DU40_A17.dat','DU35_A17.dat','DU30_A17.dat','DU25_A17.dat','DU21_A17.dat','NACA64_A17.dat']
af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7] # place at appropriate radial stations
wt.airfoil_files = [basePath]*len(wt.radial_locations)
for i in range(len(wt.radial_locations)):
    wt.airfoil_files[i] = os.path.join(basePath,airfoil_types[af_idx[i]])
# flow
    wt.wind_speed_hub = np.array([freestream_wind_speed])
    wt.air_density = air_density
    wt.viscosity = viscosity
# control DOFs
wt.yaw = 0.0
wt.pitch = np.array([0.0])
tsr = 7.55
wt.rotor_speed = wt.wind_speed_hub*tsr/wt.tip_radius * 30.0/np.pi

ccblade.turbine = wt

setattr(FLORISplusCCBlade, ccbladeName, ccblade)

Define flow properties in FLORIS

FLORISplusCCBlade.floris.wind_speed = freestream_wind_speed # m/s FLORISplusCCBlade.floris.air_density = air_density # kg/m^3 FLORISplusCCBlade.floris.wind_direction = wind_direction # deg

FLORISplusCCBlade.run()

print the CP and axial induction of each turbine

for turbineI in range(0,turbineX.size): florisTurbine = "turbine%s" % turbineI print getattr(FLORISplusCCBlade.floris.wt_layout,florisTurbine).axial_induction print getattr(FLORISplusCCBlade.floris.wt_layout,florisTurbine).CP

============= floris.py ====================== from fusedwind.plant_flow.comp import GenericWindFarm, GenericFlowModel from openmdao.lib.datatypes.api import Array, Float, VarTree from openmdao.main.api import VariableTree

import numpy as np

class FLORISParameters(VariableTree): """Container of FLORIS wake parameters"""

pP = Float(1.88, iotype='in')
ke = Float(0.065, iotype='in')
keCorrDA = Float(0, iotype='in')
kd = Float(0.15, iotype='in')
me = Array(np.array([-0.5,0.22,1.0]), iotype='in')
MU = Array(np.array([0.5,1.0,5.5]), iotype='in')
ad = Float(-4.5, iotype='in')
bd = Float(-0.01, iotype='in')
aU = Float(5.0, iotype='in', units='deg')
bU = Float(1.66, iotype='in')

class FLORIS(GenericWindFarm, GenericFlowModel): """Evaluates the FLORIS model and gives the FLORIS-predicted powers of the turbines at locations turbineX, turbineY, and, optionally, the FLORIS-predicted velocities at locations (velX,velY)"""

parameters = VarTree(FLORISParameters(), iotype='in')

def execute(self):

    # rename inputs and outputs
    Vinf = self.wind_speed                               # from standard FUSED-WIND GenericWindFarm component
    windDirection = self.wind_direction*np.pi/180        # from standard FUSED-WIND GenericWindFarm component
    rho = self.air_density

    rotorDiameter = self.wt_layout.wt_array(attr='rotor_diameter') # from standard FUSED-WIND GenericWindFarm component
    rotorArea = self.wt_layout.wt_array(attr='rotor_area')         # from standard FUSED-WIND GenericWindFarm component
    yaw = self.wt_layout.wt_array(attr='yaw')
    Cp = self.wt_layout.wt_array(attr='CP')
    axialInd = self.wt_layout.wt_array(attr='axial_induction')

    positions = self.wt_layout.wt_array(attr='position')
    turbineX = positions[:, 0]
    turbineY = positions[:, 1]

    if self.ws_positions.any():
        velX = self.ws_positions[:, 0]
        velY = self.ws_positions[:, 1]
    else:
        velX = np.zeros([0,0])
        velY = np.zeros([0,0])

    pP = self.parameters.pP
    ke = self.parameters.ke
    keCorrDA = self.parameters.keCorrDA
    kd = self.parameters.kd
    me = self.parameters.me
    MU = self.parameters.MU
    ad = self.parameters.ad
    bd = self.parameters.bd
    aU = self.parameters.aU
    bU = self.parameters.bU

    # find size of arrays
    nTurbines = turbineX.size
    nLocations = velX.size

    # convert to downwind-crosswind coordinates
    rotationMatrix = np.array([(np.cos(-windDirection), -np.sin(-windDirection)),
                               (np.sin(-windDirection), np.cos(-windDirection))])
    turbineLocations = np.dot(rotationMatrix,np.array([turbineX,turbineY]))
    turbineX = turbineLocations[0]
    turbineY = turbineLocations[1]

    if velX.size>0:
        locations = np.dot(rotationMatrix,np.array([velX,velY]))
        velX = locations[0]
        velY = locations[1]

    # calculate y-location of wake centers
    wakeCentersY = np.zeros((nLocations,nTurbines))
    wakeCentersYT = np.zeros((nTurbines,nTurbines))
    for turb in range(0,nTurbines):
        wakeAngleInit = 0.5*np.power(np.cos(yaw[turb]),2)*np.sin(yaw[turb])*4*axialInd[turb]*(1-axialInd[turb])
        for loc in range(0,nLocations):  # at velX-locations
            deltax = np.maximum(velX[loc]-turbineX[turb],0)
            factor = (2*kd*deltax/rotorDiameter[turb])+1
            wakeCentersY[loc,turb] = turbineY[turb]
            wakeCentersY[loc,turb] = wakeCentersY[loc,turb] + ad+bd*deltax  # rotation-induced deflection
            wakeCentersY[loc,turb] = wakeCentersY[loc,turb] + \
                                     (wakeAngleInit*(15*np.power(factor,4)+np.power(wakeAngleInit,2))/((30*kd*np.power(factor,5))/rotorDiameter[turb]))-(wakeAngleInit*rotorDiameter[turb]*(15+np.power(wakeAngleInit,4))/(30*kd)) # yaw-induced deflection
        for turbI in range(0,nTurbines):  # at turbineX-locations
            deltax = np.maximum(turbineX[turbI]-turbineX[turb],0)
            factor = (2*kd*deltax/rotorDiameter[turb])+1
            wakeCentersYT[turbI,turb] = turbineY[turb]
            wakeCentersYT[turbI,turb] = wakeCentersYT[turbI,turb] + ad+bd*deltax  # rotation-induced deflection
            wakeCentersYT[turbI,turb] = wakeCentersYT[turbI,turb] + \
                                        (wakeAngleInit*(15*np.power(factor,4)+np.power(wakeAngleInit,4))/((30*kd*np.power(factor,5))/rotorDiameter[turb]))-(wakeAngleInit*rotorDiameter[turb]*(15+np.power(wakeAngleInit,4))/(30*kd)) # yaw-induced deflection

    # calculate wake zone diameters at velX-locations
    wakeDiameters = np.zeros((nLocations,nTurbines,3))
    wakeDiametersT = np.zeros((nTurbines,nTurbines,3))
    for turb in range(0,nTurbines):
        for loc in range(0,nLocations):  # at velX-locations
            deltax = velX[loc]-turbineX[turb]
            for zone in range(0,3):
                wakeDiameters[loc,turb,zone] = rotorDiameter[turb]+2*ke*me[zone]*np.maximum(deltax,0)
        for turbI in range(0,nTurbines):  # at turbineX-locations
            deltax = turbineX[turbI]-turbineX[turb]
            for zone in range(0,3):
                wakeDiametersT[turbI,turb,zone] = np.maximum(rotorDiameter[turb]+2*ke*me[zone]*deltax,0)

    # calculate overlap areas at rotors
    # wakeOverlapT(TURBI,TURB,ZONEI) = overlap area of zone ZONEI of wake
    # of turbine TURB with rotor of turbine TURBI
    wakeOverlapT = calcOverlapAreas(turbineX,turbineY,rotorDiameter,wakeDiametersT,wakeCentersYT)

    # make overlap relative to rotor area (maximum value should be 1)
    wakeOverlapTRel = wakeOverlapT
    for turb in range(0,nTurbines):
        wakeOverlapTRel[turb] = wakeOverlapTRel[turb]/rotorArea[turb]

    # array effects with full or partial wake overlap:
    # use overlap area of zone 1+2 of upstream turbines to correct ke
    # Note: array effects only taken into account in calculating
    # velocity deficits, in order not to over-complicate code
    # (avoid loops in calculating overlaps)

    keUncorrected = ke
    ke = np.zeros(nTurbines)
    for turb in range(0,nTurbines):
        s = np.sum(wakeOverlapTRel[turb,:,0]+wakeOverlapTRel[turb,:,1])
        ke[turb] = keUncorrected*(1+s*keCorrDA)

    # calculate velocities in full flow field (optional)
    self.ws_array = np.tile(Vinf,nLocations)
    for turb in range(0,nTurbines):
        mU = MU/np.cos(aU*np.pi/180+bU*yaw[turb])
        for loc in range(0,nLocations):
            deltax = velX[loc] - turbineX[turb]
            radiusLoc = abs(velY[loc]-wakeCentersY[loc,turb])
            axialIndAndNearRotor = 2*axialInd[turb]

            if deltax > 0 and radiusLoc < wakeDiameters[loc,turb,0]/2.0:    # check if in zone 1
                reductionFactor = axialIndAndNearRotor*\
                                  np.power((rotorDiameter[turb]/(rotorDiameter[turb]+2*ke[turb]*(mU[0])*np.maximum(0,deltax))),2)
            elif deltax > 0 and radiusLoc < wakeDiameters[loc,turb,1]/2.0:    # check if in zone 2
                reductionFactor = axialIndAndNearRotor*\
                                  np.power((rotorDiameter[turb]/(rotorDiameter[turb]+2*ke[turb]*(mU[1])*np.maximum(0,deltax))),2)
            elif deltax > 0 and radiusLoc < wakeDiameters[loc,turb,2]/2.0:    # check if in zone 3
                reductionFactor = axialIndAndNearRotor*\
                                  np.power((rotorDiameter[turb]/(rotorDiameter[turb]+2*ke[turb]*(mU[2])*np.maximum(0,deltax))),2)
            elif deltax <= 0 and radiusLoc < rotorDiameter[turb]/2.0:     # check if axial induction zone in front of rotor
                reductionFactor = axialIndAndNearRotor*(0.5+np.arctan(2.0*np.minimum(0,deltax)/(rotorDiameter[turb]))/np.pi)
            else:
                reductionFactor = 0
            self.ws_array[loc] *= (1-reductionFactor)

    # find effective wind speeds at downstream turbines, then predict power downstream turbine
    self.velocitiesTurbines = np.tile(Vinf,nTurbines)

    for turbI in range(0,nTurbines):

        # find overlap-area weighted effect of each wake zone
        wakeEffCoeff = 0
        for turb in range(0,nTurbines):

            wakeEffCoeffPerZone = 0
            deltax = turbineX[turbI] - turbineX[turb]

            if deltax > 0:
                mU = MU / np.cos(aU*np.pi/180 + bU*yaw[turb])
                for zone in range(0,3):
                    wakeEffCoeffPerZone = wakeEffCoeffPerZone + np.power((rotorDiameter[turb])/(rotorDiameter[turb]+2*ke[turb]*mU[zone]*deltax),2.0) * wakeOverlapTRel[turbI,turb,zone]

                wakeEffCoeff = wakeEffCoeff + np.power(axialInd[turb]*wakeEffCoeffPerZone,2.0)

        wakeEffCoeff = (1 - 2 * np.sqrt(wakeEffCoeff))

        # multiply the inflow speed with the wake coefficients to find effective wind speed at turbine
        self.velocitiesTurbines[turbI] *= wakeEffCoeff

    # find turbine powers
    self.wt_power = np.power(self.velocitiesTurbines,3.0) * (0.5*rho*rotorArea*Cp*np.power(np.cos(yaw),pP))
    self.wt_power /= 1000  # in kW
    self.power = np.sum(self.wt_power)

def calcOverlapAreas(turbineX,turbineY,rotorDiameter,wakeDiameters,wakeCenters): """calculate overlap of rotors and wake zones (wake zone location defined by wake center and wake diameter) turbineX,turbineY is x,y-location of center of rotor

wakeOverlap(TURBI,TURB,ZONEI) = overlap area of zone ZONEI of wake of turbine TURB with rotor of downstream turbine
TURBI"""

nTurbines = turbineY.size

wakeOverlap = np.zeros((nTurbines,nTurbines,3))

for turb in range(0,nTurbines):
    for turbI in range(0,nTurbines):
        if turbineX[turbI] > turbineX[turb]:
            OVdYd = wakeCenters[turbI,turb]-turbineY[turbI]
            OVr = rotorDiameter[turbI]/2
            for zone in range(0,3):
                OVR = wakeDiameters[turbI,turb,zone]/2
                OVdYd = abs(OVdYd)
                if OVdYd != 0:
                    OVL = (-np.power(OVr,2.0)+np.power(OVR,2.0)+np.power(OVdYd,2.0))/(2.0*OVdYd)
                else:
                    OVL = 0

                OVz = np.power(OVR,2.0)-np.power(OVL,2.0)

                if OVz > 0:
                    OVz = np.sqrt(OVz)
                else:
                    OVz = 0

                if OVdYd < (OVr+OVR):
                    if OVL < OVR and (OVdYd-OVL) < OVr:
                        wakeOverlap[turbI,turb,zone] = np.power(OVR,2.0)*np.arccos(OVL/OVR) + np.power(OVr,2.0)*np.arccos((OVdYd-OVL)/OVr) - OVdYd*OVz
                    elif OVR > OVr:
                        wakeOverlap[turbI,turb,zone] = np.pi*np.power(OVr,2.0)
                    else:
                        wakeOverlap[turbI,turb,zone] = np.pi*np.power(OVR,2.0)
                else:
                    wakeOverlap[turbI,turb,zone] = 0

for turb in range(0,nTurbines):
    for turbI in range(0,nTurbines):
        wakeOverlap[turbI,turb,2] = wakeOverlap[turbI,turb,2]-wakeOverlap[turbI,turb,1]
        wakeOverlap[turbI,turb,1] = wakeOverlap[turbI,turb,1]-wakeOverlap[turbI,turb,0]

return wakeOverlap

================= florisRotorSE.py ==================== from openmdao.main.api import Assembly, Component from openmdao.lib.datatypes.api import Array, Float, Bool, Int, List, Str, VarTree

from rotorse.rotoraerodefaults import CCBlade from rotorse.rotoraero import Coefficients from floris import FLORIS from fusedwind.turbine.turbine_vt import AeroelasticHAWTVT from fusedwind.plant_flow.comp import GenericWindFarm from floris import FLORIS

import os import numpy as np

class AeroelasticHAWTVT_floris(AeroelasticHAWTVT):

inherits from AeroelasticHAWTVT:

# turbine_name
# rotor_diameter
# and other parameters that are not used:
#   nblades
#   hub_height
#   tilt_angle
#   cone_angle
#   hub_radius

CP = Array(np.array([0.5926]), iotype='in',desc='power coefficient')
axial_induction = Array(np.array([0.333]), iotype='in',desc='axial induction')
position = Array(iotype='in', units='m',desc='position')
rotor_area = Float(iotype='in', units='m*m', desc='rotor swept area')
name = Str(iotype='in', desc='turbine tag in wind plant')
yaw = Float(iotype='in', desc='yaw error', units='deg')

class AeroelasticHAWTVT_CCBlade(AeroelasticHAWTVT): """ AeroelasticHAWTVT with extra CCblade inputs and control inputs"""

# inherits from AeroelasticHAWTVT:
# turbine_name
# rotor_diameter
# nblades = 3 --> B
# hub_height = 90.0 --> hubHt
# tilt_angle --> tilt
# cone_angle --> precone
# hub_radius --> Rhub

# flow
tip_radius = Float(iotype='in', units='m', desc='tip radius') # Rtip
wind_speed_hub = Array(iotype='in', units='m/s', desc='hub height wind speed') # Uhub
air_density = Float(1.225, iotype='in', units='kg/m**3', desc='density of air') # rho
dynamic_viscosity = Float(1.81206e-5, iotype='in', units='kg/(m*s)', desc='dynamic viscosity of air') # mu
shear_exponent = Float(0.2, iotype='in', desc='shear exponent') # shearExp

# CCblade inputs
radial_locations = Array(iotype='in', units='m', desc='radial locations where blade is defined (should be increasing and not go all the way to hub or tip)') #r
chord = Array(iotype='in', units='m', desc='chord length at each section')
theta = Array(iotype='in', units='deg', desc='twist angle at each section (positive decreases angle of attack)')
precurve = Array(iotype='in', units='m', desc='precurve at each section')
precurveTip = Float(0.0, iotype='in', units='m', desc='precurve at tip')
airfoil_files = List(Str, iotype='in', desc='names of airfoil file')
nSector = Int(4, iotype='in', desc='number of sectors to divide rotor face into in computing thrust and power')
tiploss = Bool(True, iotype='in', desc='include Prandtl tip loss model')
hubloss = Bool(True, iotype='in', desc='include Prandtl hub loss model')
wakerotation = Bool(True, iotype='in', desc='include effect of wake rotation (i.e., tangential induction factor is nonzero)')
usecd = Bool(True, iotype='in', desc='use drag coefficient in computing induction factors')

# control DOFs
yaw = Float(iotype='in', desc='yaw error', units='deg')
pitch = Array(iotype='in', desc='blade pitch angle', units='deg')
rotor_speed = Array(iotype='in', desc='rotor speed', units='rpm') #Omega

class CCBladeAeroelasticHAWTVT(Assembly): """ CCBlade that takes AeroelasticHAWTVT_CCBlade as input """

turbine = VarTree(AeroelasticHAWTVT_CCBlade(), iotype='in')
power = Array(iotype='out', desc='rotor power', units='W')
thrust = Array(iotype='out', desc='aerodynamic thrust on rotor', units='N')
torque = Array(iotype='out', desc='aerodynamic torque on rotor', units='N*m')

def configure(self):

    self.add("CCBlade", CCBlade())

    self.driver.workflow.add(['CCBlade'])

    self.CCBlade.run_case = 'power'

    self.connect("turbine.nblades", "CCBlade.B")
    self.connect("turbine.hub_height", "CCBlade.hubHt")
    self.connect("turbine.tilt_angle", "CCBlade.tilt")
    self.connect("turbine.cone_angle", "CCBlade.precone")
    self.connect("turbine.hub_radius", "CCBlade.Rhub")
    self.connect("turbine.tip_radius", "CCBlade.Rtip")
    self.connect("turbine.wind_speed_hub", "CCBlade.Uhub")
    self.connect("turbine.air_density", "CCBlade.rho")
    self.connect("turbine.dynamic_viscosity", "CCBlade.mu")
    self.connect("turbine.shear_exponent", "CCBlade.shearExp")
    self.connect("turbine.radial_locations", "CCBlade.r")
    self.connect("turbine.chord", "CCBlade.chord")
    self.connect("turbine.theta", "CCBlade.theta")
    self.connect("turbine.precurve", "CCBlade.precurve")
    self.connect("turbine.precurveTip", "CCBlade.precurveTip")
    self.connect("turbine.airfoil_files", "CCBlade.airfoil_files")
    self.connect("turbine.nSector", "CCBlade.nSector")
    self.connect("turbine.tiploss", "CCBlade.tiploss")
    self.connect("turbine.hubloss", "CCBlade.hubloss")
    self.connect("turbine.wakerotation", "CCBlade.wakerotation")
    self.connect("turbine.usecd", "CCBlade.usecd")
    self.connect("turbine.yaw", "CCBlade.yaw")
    self.connect("turbine.pitch", "CCBlade.pitch")
    self.connect("turbine.rotor_speed", "CCBlade.Omega")

    self.connect("CCBlade.P", "power")
    self.connect("CCBlade.T", "thrust")
    self.connect("CCBlade.Q", "torque")

class CCBladeCoefficients(Assembly): """ Couples components with CCBlade so that it outputs CP and CT coefficients and axial-induction factor takes AeroelasticHAWTVT_CCBlade as input """

turbine = VarTree(AeroelasticHAWTVT_CCBlade(), iotype='in')

CP = Array(iotype='out', desc='power coefficient')
axial_induction = Array(iotype='out', desc='axial induction')

def configure(self):
    """ Creates a new Assembly containing a Paraboloid component"""

    self.add("CCBlade", CCBladeAeroelasticHAWTVT())
    self.add("Coefficients", Coefficients())
    self.add("CTtoAxialInd", CTtoAxialInd())

    self.driver.workflow.add(['CCBlade', 'Coefficients', 'CTtoAxialInd'])

    self.connect("turbine", "CCBlade.turbine")
    self.connect("CCBlade.power", "Coefficients.P")
    self.connect("CCBlade.thrust", "Coefficients.T")
    self.connect("CCBlade.torque", "Coefficients.Q")
    self.connect("turbine.wind_speed_hub", "Coefficients.V")
    self.connect("turbine.tip_radius", "Coefficients.R")
    self.connect("turbine.air_density", "Coefficients.rho")
    self.connect("Coefficients.CT", "CTtoAxialInd.CT")
    self.connect("Coefficients.CP", "CP")
    self.connect("CTtoAxialInd.axial_induction", "axial_induction")

    self.create_passthrough('Coefficients.CT')
    self.create_passthrough('Coefficients.CQ')

class CTtoAxialInd(Component): """Convert thrust coefficient to axial induction factor"""

CT = Array(iotype='in')
axial_induction = Array(iotype='out')

def execute(self):
    self.axial_induction = 0.5*(1-np.sqrt(1-self.CT))

[fzahle]

Could you put this code in a gist instead, since its now a mix of markdown formatted text and code

[pgebraad]

Frederik and Pierre suggested to use an approach with a case iterator and a aggregator. I will try that - with that. I am closing the issue.