Benchmark2.py

# Benchmark2.py
#
# NREL 5MW Wind Turbine Analysis
# Using parameters specified by the NREL report for the 5MW tower, perform 
# aerodynamic, structural and cost analysis on the turbine to verify that the 
# numbers are reasonable.
#
# Author: Lewis Li (lewisli@stanford.edu)
# Original Date: November 1st 2015

import os
from math import pi
import numpy as np
from pyopt_driver.pyopt_driver import pyOptDriver
from openmdao.lib.casehandlers.api import DumpCaseRecorder

from commonse.utilities import check_gradient_unit_test
from rotorse.rotoraero import Coefficients, SetupRunVarSpeed, \
		RegulatedPowerCurve, AEP
from rotorse.rotoraerodefaults import RotorAeroVSVPWithCCBlade,GeometrySpline, \
		CCBladeGeometry, CCBlade, CSMDrivetrain, WeibullCDF, \
		WeibullWithMeanCDF, RayleighCDF
from rotorse.precomp import Profile, Orthotropic2DMaterial, CompositeSection

from drivese.hub import HubSE
from drivese.drive import Drive4pt
from towerse.tower import TowerSE


from rotorse.rotor import RotorSE 

from WindDistribution import CalculateAEPConstantWind
from WindDistribution import CalculateAEPWeibull
from WindDistribution import ComputeScaleFunction
from WindDistribution import EstimateCapacity
from WindDistribution import ComputeLCOE

from openmdao.main.api import VariableTree, Component, Assembly, set_as_top
from openmdao.main.datatypes.api import Int, Float, Array, VarTree, Bool, Slot

from commonse.WindWaveDrag import FluidLoads, AeroHydroLoads, TowerWindDrag, TowerWaveDrag
from commonse.environment import WindBase, WaveBase  # , SoilBase
from commonse import Tube
from fusedwind.turbine.tower import TowerFromCSProps
from fusedwind.interface import implement_base
from commonse.UtilizationSupplement import fatigue, hoopStressEurocode, shellBucklingEurocode, \
    bucklingGL, vonMisesStressUtilization

import matplotlib.pyplot as plt

import frame3dd

################################################################################
### 1. Aerodynamic and structural performance using RotorSE

BladeLength = 34
HubHeight = 65

ReferenceBladeLength = 68;
rotor = RotorSE()
# -------------------

# === blade grid ===
# (Array): initial aerodynamic grid on unit radius
rotor.initial_aero_grid = np.array([0.02222276, 0.06666667, 0.11111057, \
	0.16666667, 0.23333333, 0.3, 0.36666667, 0.43333333, 0.5, 0.56666667, \
	0.63333333, 0.7, 0.76666667, 0.83333333, 0.88888943, 0.93333333, \
    0.97777724]) 

 # (Array): initial structural grid on unit radius
rotor.initial_str_grid = np.array([0.0, 0.00492790457512, 0.00652942887106, 
	0.00813095316699, 0.00983257273154, 0.0114340970275, 0.0130356213234, 
	0.02222276, 0.024446481932, 0.026048006228, 0.06666667, 0.089508406455,
    0.11111057, 0.146462614229, 0.16666667, 0.195309105255, 0.23333333, 
    0.276686558545, 0.3, 0.333640766319,0.36666667, 0.400404310407, 0.43333333, 
    0.5, 0.520818918408, 0.56666667, 0.602196371696, 0.63333333,
    0.667358391486, 0.683573824984, 0.7, 0.73242031601, 0.76666667, 0.83333333, 
    0.88888943, 0.93333333, 0.97777724, 1.0]) 

# (Int): first idx in r_aero_unit of non-cylindrical section, 
# constant twist inboard of here
rotor.idx_cylinder_aero = 3  

# (Int): first idx in r_str_unit of non-cylindrical section
rotor.idx_cylinder_str = 14  

# (Float): hub location as fraction of radius
rotor.hubFraction = 0.025  
# ------------------

# === blade geometry ===
# (Array): new aerodynamic grid on unit radius
rotor.r_aero = np.array([0.02222276, 0.06666667, 0.11111057, 0.2, 0.23333333, 
	0.3, 0.36666667, 0.43333333, 0.5, 0.56666667, 0.63333333, 0.64, 0.7, 
	0.83333333, 0.88888943, 0.93333333, 0.97777724])  

# (Float): location of max chord on unit radius
rotor.r_max_chord = 0.23577

# (Array, m): chord at control points. defined at hub, then at linearly spaced
# locations from r_max_chord to tip
ReferenceChord = [3.2612, 4.5709, 3.3178, 1.4621]
rotor.chord_sub = [x * np.true_divide(BladeLength,ReferenceBladeLength) for x in ReferenceChord]

# (Array, deg): twist at control points.  defined at linearly spaced locations 
# from r[idx_cylinder] to tip
rotor.theta_sub = [13.2783, 7.46036, 2.89317, -0.0878099]  

# (Array, m): precurve at control points.  defined at same locations at chord, 
# starting at 2nd control point (root must be zero precurve)
rotor.precurve_sub = [0.0, 0.0, 0.0] 

# (Array, m): adjustment to precurve to account for curvature from loading
rotor.delta_precurve_sub = [0.0, 0.0, 0.0]  

# (Array, m): spar cap thickness parameters
rotor.sparT = [0.05, 0.047754, 0.045376, 0.031085, 0.0061398] 

# (Array, m): trailing-edge thickness parameters
rotor.teT = [0.1, 0.09569, 0.06569, 0.02569, 0.00569]  

# (Float, m): blade length (if not precurved or swept) 
# otherwise length of blade before curvature
rotor.bladeLength = BladeLength  

# (Float, m): adjustment to blade length to account for curvature from loading
rotor.delta_bladeLength = 0.0  
rotor.precone = 2.5  # (Float, deg): precone angle
rotor.tilt = 5.0  # (Float, deg): shaft tilt
rotor.yaw = 0.0  # (Float, deg): yaw error
rotor.nBlades = 3  # (Int): number of blades
# ------------------

# === airfoil files ===
basepath = os.path.join(os.path.dirname(os.path.realpath(__file__)), '5MW_AFFiles')

# load all airfoils
airfoil_types = [0]*8
airfoil_types[0] = os.path.join(basepath, 'Cylinder1.dat')
airfoil_types[1] = os.path.join(basepath, 'Cylinder2.dat')
airfoil_types[2] = os.path.join(basepath, 'DU40_A17.dat')
airfoil_types[3] = os.path.join(basepath, 'DU35_A17.dat')
airfoil_types[4] = os.path.join(basepath, 'DU30_A17.dat')
airfoil_types[5] = os.path.join(basepath, 'DU25_A17.dat')
airfoil_types[6] = os.path.join(basepath, 'DU21_A17.dat')
airfoil_types[7] = os.path.join(basepath, 'NACA64_A17.dat')

# place at appropriate radial stations
af_idx = [0, 0, 1, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7, 7, 7]

n = len(af_idx)
af = [0]*n
for i in range(n):
    af[i] = airfoil_types[af_idx[i]]
rotor.airfoil_files = af  # (List): names of airfoil file
# ----------------------

# === atmosphere ===
rotor.rho = 1.225  # (Float, kg/m**3): density of air
rotor.mu = 1.81206e-5  # (Float, kg/m/s): dynamic viscosity of air
rotor.shearExp = 0.25  # (Float): shear exponent
rotor.hubHt = HubHeight  # (Float, m): hub height
rotor.turbine_class = 'I'  # (Enum): IEC turbine class
rotor.turbulence_class = 'B'  # (Enum): IEC turbulence class class
rotor.cdf_reference_height_wind_speed = 30.0 
rotor.g = 9.81  # (Float, m/s**2): acceleration of gravity
# ----------------------

# === control ===
rotor.control.Vin = 3.0  # (Float, m/s): cut-in wind speed
rotor.control.Vout = 26.0  # (Float, m/s): cut-out wind speed
rotor.control.ratedPower = 1.5e6  # (Float, W): rated power
rotor.control.minOmega = 0.0  # (Float, rpm): minimum allowed rotor rotation speed
rotor.control.maxOmega = 24.0  # (Float, rpm): maximum allowed rotor rotation speed
rotor.control.tsr = 7  # (Float): tip-speed ratio in Region 2 (should be optimized externally)
rotor.control.pitch = 0.0  # (Float, deg): pitch angle in region 2 (and region 3 for fixed pitch machines)
rotor.pitch_extreme = 0.0  # (Float, deg): worst-case pitch at survival wind condition
rotor.azimuth_extreme = 0.0  # (Float, deg): worst-case azimuth at survival wind condition
rotor.VfactorPC = 0.7  # (Float): fraction of rated speed at which the deflection is assumed to representative throughout the power curve calculation
# ----------------------

# === aero and structural analysis options ===
rotor.nSector = 4  # (Int): number of sectors to divide rotor face into in computing thrust and power
rotor.npts_coarse_power_curve = 20  # (Int): number of points to evaluate aero analysis at
rotor.npts_spline_power_curve = 200  # (Int): number of points to use in fitting spline to power curve
rotor.AEP_loss_factor = 1.0  # (Float): availability and other losses (soiling, array, etc.)
rotor.drivetrainType = 'geared'  # (Enum)
rotor.nF = 5  # (Int): number of natural frequencies to compute
rotor.dynamic_amplication_tip_deflection = 1.35  # (Float): a dynamic amplification factor to adjust the static deflection calculation
# ----------------------

# === materials and composite layup  ===
basepath = os.path.join(os.path.dirname(os.path.realpath(__file__)), '5MW_PrecompFiles')

materials = Orthotropic2DMaterial.listFromPreCompFile(os.path.join(basepath, 'materials.inp'))

ncomp = len(rotor.initial_str_grid)
upper = [0]*ncomp
lower = [0]*ncomp
webs = [0]*ncomp
profile = [0]*ncomp

# (Array): array of leading-edge positions from a reference blade axis 
# (usually blade pitch axis). locations are normalized by the local chord 
# length. e.g. leLoc[i] = 0.2 means leading edge is 0.2*chord[i] from reference 
# axis.  positive in -x direction for airfoil-aligned coordinate system
rotor.leLoc = np.array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.498, 0.497, 
	0.465, 0.447, 0.43, 0.411, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 
	0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4])  

# (Array): index of sector for spar (PreComp definition of sector)
rotor.sector_idx_strain_spar = [2]*ncomp  

# (Array): index of sector for trailing-edge (PreComp definition of sector)
rotor.sector_idx_strain_te = [3]*ncomp  

web1 = np.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.4114, 0.4102, 
	0.4094, 0.3876, 0.3755, 0.3639, 0.345, 0.3342, 0.3313, 0.3274, 0.323, 
	0.3206, 0.3172, 0.3138, 0.3104, 0.307, 0.3003, 0.2982, 0.2935, 0.2899, 
	0.2867, 0.2833, 0.2817, 0.2799, 0.2767, 0.2731, 0.2664, 0.2607, 0.2562, 
	0.1886, -1.0])

web2 = np.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.5886, 0.5868, 
	0.5854, 0.5508, 0.5315, 0.5131, 0.4831, 0.4658, 0.4687, 0.4726, 0.477, 
	0.4794, 0.4828, 0.4862, 0.4896, 0.493, 0.4997, 0.5018, 0.5065, 0.5101, 
	0.5133, 0.5167, 0.5183, 0.5201, 0.5233, 0.5269, 0.5336, 0.5393, 0.5438, 
	0.6114, -1.0])
web3 = np.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 
	-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 
	1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 
	-1.0, -1.0])

# (Array, m): chord distribution for reference section, thickness of structural 
# layup scaled with reference thickness (fixed t/c for this case)
rotor.chord_str_ref = np.array([3.2612, 3.3100915356, 3.32587052924, 
	3.34159388653, 3.35823798667, 3.37384375335, 3.38939112914, 3.4774055542, 
	3.49839685, 3.51343645709, 3.87017220335, 4.04645623801, 4.19408216643,
    4.47641008477, 4.55844487985, 4.57383098262, 4.57285771934, 4.51914315648, 
    4.47677655262, 4.40075650022, 4.31069949379, 4.20483735936, 4.08985563932, 
    3.82931757126, 3.74220276467, 3.54415796922, 3.38732428502, 3.24931446473, 
    3.23421422609, 3.22701537997, 3.21972125648, 3.08979310611, 2.95152261813, 
    2.330753331, 2.05553464181, 1.82577817774, 1.5860853279, 1.4621])* \
	np.true_divide(BladeLength,ReferenceBladeLength)


for i in range(ncomp):
    webLoc = []
    if web1[i] != -1:
        webLoc.append(web1[i])
    if web2[i] != -1:
        webLoc.append(web2[i])
    if web3[i] != -1:
        webLoc.append(web3[i])

    upper[i], lower[i], webs[i] = CompositeSection.initFromPreCompLayupFile\
    (os.path.join(basepath, 'layup_' + str(i+1) + '.inp'), webLoc, materials)
    profile[i] = Profile.initFromPreCompFile(os.path.join(basepath, 'shape_' \
    	+ str(i+1) + '.inp'))

# (List): list of all Orthotropic2DMaterial objects used in defining the geometry
rotor.materials = materials 

# (List): list of CompositeSection objections defining the properties for upper surface
rotor.upperCS = upper  

# (List): list of CompositeSection objections defining the properties for lower surface
rotor.lowerCS = lower  

# (List): list of CompositeSection objections defining the properties for shear webs
rotor.websCS = webs  

# (List): airfoil shape at each radial position
rotor.profile = profile  
# --------------------------------------


# === fatigue ===

# (Array): nondimensional radial locations of damage equivalent moments
rotor.rstar_damage = np.array([0.000, 0.022, 0.067, 0.111, 0.167, 0.233, 0.300,
 0.367, 0.433, 0.500, 0.567, 0.633, 0.700, 0.767, 0.833, 0.889, 0.933, 0.978]) 


# (Array, N*m): damage equivalent moments about blade c.s. x-direction
rotor.Mxb_damage = 1e3*np.array([2.3743E+003, 2.0834E+003, 1.8108E+003, 
	1.5705E+003, 1.3104E+003, 1.0488E+003, 8.2367E+002, 6.3407E+002, 
	4.7727E+002, 3.4804E+002, 2.4458E+002, 1.6339E+002, 1.0252E+002,
	 5.7842E+001, 2.7349E+001, 1.1262E+001, 3.8549E+000, 4.4738E-001])  

# (Array, N*m): damage equivalent moments about blade c.s. y-direction
rotor.Myb_damage = 1e3*np.array([2.7732E+003, 2.8155E+003, 2.6004E+003, 
	2.3933E+003, 2.1371E+003, 1.8459E+003, 1.5582E+003, 1.2896E+003, 
	1.0427E+003, 8.2015E+002, 6.2449E+002, 4.5229E+002, 3.0658E+002, 
	1.8746E+002, 9.6475E+001, 4.2677E+001, 1.5409E+001, 1.8426E+000])  

rotor.strain_ult_spar = 1.0e-2  # (Float): ultimate strain in spar cap

# (Float): uptimate strain in trailing-edge panels, note that I am putting a 
# factor of two for the damage part only.
rotor.strain_ult_te = 2500*1e-6 * 2  
rotor.eta_damage = 1.35*1.3*1.0  # (Float): safety factor for fatigue
rotor.m_damage = 10.0  # (Float): slope of S-N curve for fatigue analysis

# (Float): number of cycles used in fatigue analysis  
rotor.N_damage = 365*24*3600*20.0  
# ----------------

# from myutilities import plt

# === run and outputs ===
rotor.run()

# Evaluate AEP Using Lewis' Functions
# Weibull Wind Parameters
WindReferenceHeight = 50
WindReferenceMeanVelocity = 7.5
WeibullShapeFactor = 2.0
ShearFactor = 0.25

PowerCurve = rotor.P/1e6
PowerCurveVelocity = rotor.V

HubHeight = rotor.hubHt

AEP,WeibullScale = CalculateAEPWeibull(PowerCurve,PowerCurveVelocity, HubHeight, \
  	BladeLength,WeibullShapeFactor, WindReferenceHeight, \
  	WindReferenceMeanVelocity, ShearFactor)

NamePlateCapacity = EstimateCapacity(PowerCurve,PowerCurveVelocity, rotor.ratedConditions.V)

#AEP = CalculateAEPConstantWind(PowerCurve, PowerCurveVelocity, 7.5)

print 'AEP = %d MWH' %(AEP)
print 'NamePlateCapacity = %fMW' %(NamePlateCapacity)
print 'diameter =', rotor.diameter
print 'ratedConditions.V =', rotor.ratedConditions.V
print 'ratedConditions.Omega =', rotor.ratedConditions.Omega
print 'ratedConditions.pitch =', rotor.ratedConditions.pitch
print 'mass_one_blade =', rotor.mass_one_blade
print 'mass_all_blades =', rotor.mass_all_blades
print 'I_all_blades =', rotor.I_all_blades
print 'freq =', rotor.freq
print 'tip_deflection =', rotor.tip_deflection
print 'root_bending_moment =', rotor.root_bending_moment
print '##########################################'

################################################################################
### 2. Hub Sizing 
# Specify hub parameters based off rotor

# Load default hub model
hubS = HubSE()
hubS.rotor_diameter = rotor.Rtip*2 # m
hubS.blade_number  = rotor.nBlades
hubS.blade_root_diameter = rotor.chord_sub[0]*1.25
hubS.L_rb = rotor.hubFraction*rotor.diameter
hubS.MB1_location = np.array([-0.5, 0.0, 0.0])
hubS.machine_rating = rotor.control.ratedPower
hubS.blade_mass = rotor.mass_one_blade
hubS.rotor_bending_moment = rotor.root_bending_moment

hubS.run()
print "Estimate of Hub Component Sizes for the Baseline 5 MW Reference Turbine"
print "Hub Components"
print '  Hub: {0:8.1f} kg'.format(hubS.hub.mass)  # 31644.47
print '  Pitch system: {0:8.1f} kg'.format(hubS.pitchSystem.mass) # 17003.98
print '  Nose cone: {0:8.1f} kg'.format(hubS.spinner.mass) # 1810.50

RotorTotalWeight = rotor.mass_all_blades + hubS.spinner.mass + \
hubS.hub.mass + hubS.pitchSystem.mass

print 'Rotor Total Weight = %d kg' %RotorTotalWeight

print '##########################################'

################################################################################
### 3. Drive train + Nacelle Mass estimation
# NREL 5 MW Rotor Variables
nace = Drive4pt()
nace.rotor_diameter = rotor.Rtip *2 # m
nace.rotor_speed = rotor.ratedConditions.Omega # #rpm m/s
nace.machine_rating = hubS.machine_rating/1000
nace.DrivetrainEfficiency = 0.95

 # 6.35e6 #4365248.74 # Nm
nace.rotor_torque =  rotor.ratedConditions.Q
nace.rotor_thrust = rotor.ratedConditions.T # N
nace.rotor_mass = 0.0 #accounted for in F_z # kg

nace.rotor_bending_moment_x = rotor.Mxyz_0[0]
nace.rotor_bending_moment_y = rotor.Mxyz_0[1]
nace.rotor_bending_moment_z = rotor.Mxyz_0[2]

nace.rotor_force_x = rotor.Fxyz_0[0] # N
nace.rotor_force_y = rotor.Fxyz_0[1]
nace.rotor_force_z = rotor.Fxyz_0[2] # N

# NREL 5 MW Drivetrain variables
nace.drivetrain_design = 'geared' # geared 3-stage Gearbox with induction generator machine
nace.gear_ratio = 96.76 # 97:1 as listed in the 5 MW reference document
nace.gear_configuration = 'eep' # epicyclic-epicyclic-parallel

nace.crane = True # onboard crane present
nace.shaft_angle = 5.0 #deg
nace.shaft_ratio = 0.10
nace.Np = [3,3,1]
nace.ratio_type = 'optimal'
nace.shaft_type = 'normal'
nace.uptower_transformer=False
nace.shrink_disc_mass = 333.3*nace.machine_rating/1000.0 # estimated
nace.mb1Type = 'CARB'
nace.mb2Type = 'SRB'
nace.flange_length = 0.5 #m
nace.overhang = 5.0
nace.gearbox_cm = 0.1
nace.hss_length = 1.5
nace.check_fatigue = 0 #0 if no fatigue check, 1 if parameterized fatigue check, 2 if known loads inputs
nace.blade_number=rotor.nBlades
nace.cut_in=rotor.control.Vin #cut-in m/s
nace.cut_out=rotor.control.Vout #cut-out m/s
nace.Vrated=rotor.ratedConditions.V #rated windspeed m/s
nace.weibull_k = WeibullShapeFactor # windepeed distribution shape parameter

# Might need to change this...
nace.weibull_A = WeibullScale  # windspeed distribution scale parameter

nace.T_life=20. #design life in years
nace.IEC_Class_Letter = 'B'
nace.L_rb = hubS.L_rb # length from hub center to main bearing, leave zero if unknown

# NREL 5 MW Tower Variables
nace.tower_top_diameter = 3.78 # m

nace.run()

print "Estimate of Nacelle Component Sizes for the NREL 5 MW Reference Turbine"
print 'Low speed shaft: {0:8.1f} kg'.format(nace.lowSpeedShaft.mass)
print 'Main bearings: {0:8.1f} kg'.format(nace.mainBearing.mass + nace.secondBearing.mass)
print 'Gearbox: {0:8.1f} kg'.format(nace.gearbox.mass)
print 'High speed shaft & brakes: {0:8.1f} kg'.format(nace.highSpeedSide.mass)
print 'Generator: {0:8.1f} kg'.format(nace.generator.mass)
print 'Variable speed electronics: {0:8.1f} kg'.format(nace.above_yaw_massAdder.vs_electronics_mass)
print 'Overall mainframe:{0:8.1f} kg'.format(nace.above_yaw_massAdder.mainframe_mass)
print '     Bedplate: {0:8.1f} kg'.format(nace.bedplate.mass)
print 'Electrical connections: {0:8.1f} kg'.format(nace.above_yaw_massAdder.electrical_mass)
print 'HVAC system: {0:8.1f} kg'.format(nace.above_yaw_massAdder.hvac_mass )
print 'Nacelle cover: {0:8.1f} kg'.format(nace.above_yaw_massAdder.cover_mass)
print 'Yaw system: {0:8.1f} kg'.format(nace.yawSystem.mass)
print 'Overall nacelle: {0:8.1f} kg'.format(nace.nacelle_mass, nace.nacelle_cm[0], nace.nacelle_cm[1], nace.nacelle_cm[2], nace.nacelle_I[0], nace.nacelle_I[1], nace.nacelle_I[2]  )


################################################################################
### 4. Tower Mass

print '##########################################'

# --- tower setup ------
from commonse.environment import PowerWind

tower = set_as_top(TowerSE())

# ---- tower ------
tower.replace('wind1', PowerWind())
tower.replace('wind2', PowerWind())
# onshore (no waves)

# --- geometry ----
ReferenceTowerHeight = 95
tower.z_param = [0.0, HubHeight*0.5, HubHeight]
TowerRatio = np.true_divide(HubHeight,ReferenceTowerHeight)

tower.d_param = [6.0*TowerRatio, 4.935*TowerRatio, 3.87*TowerRatio]
tower.t_param = [0.027*1.3*TowerRatio, 0.023*1.3*TowerRatio, 0.019*1.3*TowerRatio]
n = 10

tower.z_full = np.linspace(0.0, HubHeight, n)
tower.L_reinforced = 15.0*np.ones(n)  # [m] buckling length
tower.theta_stress = 0.0*np.ones(n)
tower.yaw = 0.0

# --- material props ---
tower.E = 210e9*np.ones(n)
tower.G = 80.8e9*np.ones(n)
tower.rho = 8500.0*np.ones(n)
tower.sigma_y = 450.0e6*np.ones(n)

# --- spring reaction data.  Use float('inf') for rigid constraints. ---
tower.kidx = [0]  # applied at base
tower.kx = [float('inf')]
tower.ky = [float('inf')]
tower.kz = [float('inf')]
tower.ktx = [float('inf')]
tower.kty = [float('inf')]
tower.ktz = [float('inf')]

# --- extra mass ----
tower.midx = [n-1]  # RNA mass at top
tower.m = [0.8]
tower.mIxx = [1.14930678e+08]
tower.mIyy = [2.20354030e+07]
tower.mIzz = [1.87597425e+07]
tower.mIxy = [0.00000000e+00]
tower.mIxz = [5.03710467e+05]
tower.mIyz = [0.00000000e+00]
tower.mrhox = [-1.13197635]
tower.mrhoy = [0.]
tower.mrhoz = [0.50875268]
tower.addGravityLoadForExtraMass = False
# -----------

# --- wind ---
tower.wind_zref = 90.0
tower.wind_z0 = 0.0
tower.wind1.shearExp = 0.14
tower.wind2.shearExp = 0.14
# ---------------

# if addGravityLoadForExtraMass=True be sure not to double count by adding those force here also
# # --- loading case 1: max Thrust ---
tower.wind_Uref1 = 11.73732
tower.plidx1 = [n-1]  # at tower top
tower.Fx1 = [0.19620519]
tower.Fy1 = [0.]
tower.Fz1 = [-2914124.84400512]
tower.Mxx1 = [3963732.76208099]
tower.Myy1 = [-2275104.79420872]
tower.Mzz1 = [-346781.68192839]
# # ---------------

# # --- loading case 2: max wind speed ---
tower.wind_Uref2 = 70.0
tower.plidx1 = [n-1]  # at tower top
tower.Fx1 = [930198.60063279]
tower.Fy1 = [0.]
tower.Fz1 = [-2883106.12368949]
tower.Mxx1 = [-1683669.22411597]
tower.Myy1 = [-2522475.34625363]
tower.Mzz1 = [147301.97023764]
# # ---------------

# # --- run ---
tower.run()

print 'mass (kg) =', tower.mass
print 'f1 (Hz) =', tower.f1
print 'f2 (Hz) =', tower.f2
print 'top_deflection1 (m) =', tower.top_deflection1
print 'top_deflection2 (m) =', tower.top_deflection2

################################################################################
## 5. Turbine captial costs analysis

from turbine_costsse.turbine_costsse.turbine_costsse import Turbine_CostsSE

turbine = Turbine_CostsSE()

# NREL 5 MW turbine component masses based on Sunderland model approach
# Rotor
turbine.blade_mass = rotor.mass_one_blade  # inline with the windpact estimates
turbine.hub_mass = hubS.hub.mass
turbine.pitch_system_mass = hubS.pitchSystem.mass
turbine.spinner_mass = hubS.spinner.mass

# Drivetrain and Nacelle
turbine.low_speed_shaft_mass = nace.lowSpeedShaft.mass
turbine.main_bearing_mass=nace.mainBearing.mass 
turbine.second_bearing_mass = nace.secondBearing.mass
turbine.gearbox_mass = nace.gearbox.mass
turbine.high_speed_side_mass = nace.highSpeedSide.mass
turbine.generator_mass = nace.generator.mass
turbine.bedplate_mass = nace.bedplate.mass
turbine.yaw_system_mass = nace.yawSystem.mass

# Tower
turbine.tower_mass = tower.mass*0.5

# Additional non-mass cost model input variables
turbine.machine_rating = hubS.machine_rating/1000
turbine.advanced = False
turbine.blade_number = rotor.nBlades
turbine.drivetrain_design = 'geared'
turbine.crane = False
turbine.offshore = False

# Target year for analysis results
turbine.year = 2010
turbine.month =  12

turbine.run()

print '##########################################'
print "Overall rotor cost with 3 advanced blades is ${0:.2f} USD".format(turbine.rotorCC.cost)
print "Blade cost is ${0:.2f} USD".format(turbine.rotorCC.bladeCC.cost)
print "Hub cost is ${0:.2f} USD".format(turbine.rotorCC.hubCC.cost)
print "Pitch system cost is ${0:.2f} USD".format(turbine.rotorCC.pitchSysCC.cost)
print "Spinner cost is ${0:.2f} USD".format(turbine.rotorCC.spinnerCC.cost)
print
print "Overall nacelle cost is ${0:.2f} USD".format(turbine.nacelleCC.cost)
print "LSS cost is ${0:.2f} USD".format(turbine.nacelleCC.lssCC.cost)
print "Main bearings cost is ${0:.2f} USD".format(turbine.nacelleCC.bearingsCC.cost)
print "Gearbox cost is ${0:.2f} USD".format(turbine.nacelleCC.gearboxCC.cost)
print "Hight speed side cost is ${0:.2f} USD".format(turbine.nacelleCC.hssCC.cost)
print "Generator cost is ${0:.2f} USD".format(turbine.nacelleCC.generatorCC.cost)
print "Bedplate cost is ${0:.2f} USD".format(turbine.nacelleCC.bedplateCC.cost)
print "Yaw system cost is ${0:.2f} USD".format(turbine.nacelleCC.yawSysCC.cost)
print
print "Tower cost is ${0:.2f} USD".format(turbine.towerCC.cost)
print
print "The overall turbine cost is ${0:.2f} USD".format(turbine.turbine_cost)
print


################################################################################
## 5. Operating Expenses

# A simple test of nrel_csm_bos model
from plant_costsse.nrel_csm_bos.nrel_csm_bos import bos_csm_assembly

bos = bos_csm_assembly()

# Set input parameters
bos = bos_csm_assembly()
bos.machine_rating = hubS.machine_rating/1000
bos.rotor_diameter = rotor.diameter
bos.turbine_cost = turbine.turbine_cost
bos.hub_height = HubHeight
bos.turbine_number = 1
bos.sea_depth = 0
bos.year = 2009
bos.month = 12
bos.multiplier = 1.0
bos.run()

print "BOS cost per turbine: ${0:.2f} USD".format(bos.bos_costs / bos.turbine_number)

from plant_costsse.nrel_csm_opex.nrel_csm_opex import opex_csm_assembly
om = opex_csm_assembly()

om.machine_rating = rotor.control.ratedPower/1000  # Need to manipulate input or underlying component will not execute
om.net_aep = AEP*10e4
om.sea_depth = 0
om.year = 2009
om.month = 12
om.turbine_number = 100

om.run()

print "Average annual operational expenditures for an on wind plant with 100 turbines"
print "OPEX offshore ${:.2f}: USD".format(om.avg_annual_opex)
print "Preventative OPEX by turbine: ${:.2f} USD".format(om.opex_breakdown.preventative_opex / om.turbine_number)
print "Corrective OPEX by turbine: ${:.2f} USD".format(om.opex_breakdown.corrective_opex / om.turbine_number)
print "Land Lease OPEX by turbine: ${:.2f} USD".format(om.opex_breakdown.lease_opex / om.turbine_number)

CapitalCost = turbine.turbine_cost + bos.bos_costs / bos.turbine_number
OperatingCost = om.opex_breakdown.preventative_opex / om.turbine_number + \
om.opex_breakdown.lease_opex / om.turbine_number + \
om.opex_breakdown.corrective_opex / om.turbine_number

Years = 25
DiscountRate = 0.08

LCOE = ComputeLCOE(AEP, CapitalCost, OperatingCost, DiscountRate, Years)

print "Levelized Cost of Energy over %d years \
is $%f/kWH" %(Years,LCOE/1000)