Commit 2c452ad4 authored by Frederik Zahle's avatar Frederik Zahle
Browse files

cleanup

parent b6a73ba2
......@@ -8,23 +8,23 @@ import numpy as np
from scipy.optimize import minimize
from openmdao.main.api import Component
from openmdao.lib.datatypes.api import Float, Array, Bool
from openmdao.lib.datatypes.api import Float, Array, Bool, Int
class SEAMRotor(Component):
class SEAMBladeStructure(Component):
Nsections = Float(iotype='in', desc='number of sections')
Nsections = Int(iotype='in', desc='number of sections')
Neq = Float(1.e7, iotype='in', desc='')
WohlerExpFlap = Float(iotype='in', desc='')
PMtarget = Float(iotype='in', desc='')
Diameter = Float(iotype='in', units='m', desc='') #[m]
rotor_diameter = Float(iotype='in', units='m', desc='') #[m]
# RootChord = Float(iotype='in', units='m', desc='') #[m]
# MaxChord = Float(iotype='in', units='m', desc='') #[m]
MaxChordrR = Float(iotype='in', units='m', desc='') #[m]
OverallMaxFlap = Float(iotype='in', desc='')
OverallMaxEdge = Float(iotype='in', desc='')
overallMaxFlap = Float(iotype='in', desc='')
overallMaxEdge = Float(iotype='in', desc='')
TIF_FLext = Float(iotype='in', desc='') # Tech Impr Factor _ flap extreme
TIF_EDext = Float(iotype='in', desc='')
......@@ -48,33 +48,39 @@ class SEAMRotor(Component):
BladeWeight = Float(iotype = 'out', units = 'kg', desc = 'BladeMass' )
def execute(self):
volumen = 0. # [m^3]
r = np.linspace(0, self.Diameter/2., self.Nsections)
volumen = 0.
r = np.linspace(0, self.rotor_diameter/2., self.Nsections)
C = np.zeros(r.shape)
thick = np.zeros(r.shape)
RootChord = (self.Diameter/2.)/25.
RootChord = (self.rotor_diameter/2.)/25.
if RootChord < 1:
RootChord = 1
MaxChord = 0.001*(self.Diameter/2.)**2-0.0354*(self.Diameter/2.)+1.6635 #Empirical
MaxChord = 0.001*(self.rotor_diameter/2.)**2-0.0354*(self.rotor_diameter/2.)+1.6635 #Empirical
if MaxChord < 1.5:
RootChord = 1.5
for i in range(0,int(self.Nsections)):
for i in range(self.Nsections):
# Calculating the chord and thickness
if r[i] > (self.MaxChordrR*self.Diameter/2.):
C[i] = MaxChord-(MaxChord)/((self.Diameter/2.)-self.MaxChordrR*self.Diameter/2.)*(r[i]-(self.MaxChordrR*self.Diameter/2.))
thick[i] = 0.4*MaxChord-(0.4*MaxChord-0.05*MaxChord)/(self.Diameter/2.-self.MaxChordrR*self.Diameter/2.)*(r[i]-self.MaxChordrR*self.Diameter/2.)
if r[i] > (self.MaxChordrR*self.rotor_diameter/2.):
C[i] = MaxChord - (MaxChord)/((self.rotor_diameter/2.) - \
self.MaxChordrR*self.rotor_diameter/2.)*(r[i] - \
(self.MaxChordrR*self.rotor_diameter/2.))
thick[i] = 0.4*MaxChord - \
(0.4*MaxChord-0.05*MaxChord)/(self.rotor_diameter/2. - \
self.MaxChordrR*self.rotor_diameter/2.)*(r[i] - \
self.MaxChordrR*self.rotor_diameter/2.)
if thick[i]<0.001:
thick[i] = 0.001
if C[i]<0.001:
C[i] = 0.001
else:
C[i] = MaxChord-(MaxChord)/((self.Diameter/2.)-self.MaxChordrR*(self.Diameter/2.))*(self.MaxChordrR*(self.Diameter/2.)-r[i])
C[i] = MaxChord-(MaxChord)/((self.rotor_diameter/2.) - \
self.MaxChordrR*(self.rotor_diameter/2.)) * \
(self.MaxChordrR*(self.rotor_diameter/2.)-r[i])
thick[i] = RootChord
if thick[i]<0.001:
......@@ -82,80 +88,66 @@ class SEAMRotor(Component):
if C[i]<0.001:
C[i] = 0.001
# Printing radius, thickness and chord
print 'radius\n', r
print
print 'thickness\n', thick
print
print 'chord\n', C
print
self.thick = thick
self.chord = C
# Calculating load flapwise extreme
Mext_flap = (self.OverallMaxFlap-1.75*self.OverallMaxFlap*r/(self.Diameter/2.)+0.75*self.OverallMaxFlap*r/(self.Diameter/2.)*r/(self.Diameter/2.))*self.TIF_FLext
Mext_flap = (self.overallMaxFlap-1.75*self.overallMaxFlap*r/(self.rotor_diameter/2.) +\
0.75*self.overallMaxFlap*r/(self.rotor_diameter/2.) * r / \
(self.rotor_diameter/2.))*self.TIF_FLext
# Calculating load edgewise extreme
Mext_edge = (self.OverallMaxEdge-1.75*self.OverallMaxEdge*r/(self.Diameter/2.)+0.75*self.OverallMaxEdge*r/(self.Diameter/2.)*r/(self.Diameter/2.))*self.TIF_EDext
Mext_edge = (self.overallMaxEdge-1.75*self.overallMaxEdge*r/(self.rotor_diameter/2.) + \
0.75*self.overallMaxEdge*r/(self.rotor_diameter/2.) * r / \
(self.rotor_diameter/2.))*self.TIF_EDext
# Calculating load flapwise fatigue
Mfat_flap = (self.FlapLEQ-1.75*self.FlapLEQ*r/(self.Diameter/2.)+0.75*self.FlapLEQ*r/(self.Diameter/2.)*r/(self.Diameter/2.))*self.TIF_FLfat
Mfat_flap = (self.FlapLEQ-1.75*self.FlapLEQ*r/(self.rotor_diameter/2.) + \
0.75*self.FlapLEQ*r/(self.rotor_diameter/2.)*r / \
(self.rotor_diameter/2.))*self.TIF_FLfat
# Calculating load edgewise fatigue
Mfat_edge = (self.EdgeLEQ-1.75*self.EdgeLEQ*r/(self.Diameter/2.)+0.75*self.EdgeLEQ*r/(self.Diameter/2.)*r/(self.Diameter/2.))*self.TIF_EDext
# Printing the 4 load cases
print 'Mext_flap\n', Mext_flap
print
print 'Mext_edge\n', Mext_edge
print
print 'Mfat_flap\n', Mfat_flap
print
print 'Mfat_edge\n', Mfat_edge
print
Mfat_edge = (self.EdgeLEQ-1.75*self.EdgeLEQ*r/(self.rotor_diameter/2.) + \
0.75*self.EdgeLEQ*r/(self.rotor_diameter/2.)*r/(self.rotor_diameter/2.))*self.TIF_EDext
self.Mext_flap = Mext_flap
self.Mext_edge = Mext_edge
self.Mfat_flap = Mfat_flap
self.Mfat_edge = Mfat_edge
# Calculating thickness flapwise extreme
text_flap = np.zeros(self.Nsections)
for i in range(int(self.Nsections)):
res = minimize(self.solve_text_flap, 0.03, args = (C[i], thick[i], Mext_flap[i]), bounds = None, method = 'COBYLA', tol = 1.e-12, options={'maxiter':10000})
for i in range(self.Nsections):
norm = self.solve_text_flap(0.01, C[i], thick[i], Mext_flap[i], 1.)
res = minimize(self.solve_text_flap, 0.01, args = (C[i], thick[i], Mext_flap[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
text_flap[i] = res['x']
#print 'done', text_flap[i]
# if not res['success']: print 'WARNING', res
if not res['success']: print 'WARNING enter text_flap', i,res
# Calculating thickness edgewise extreme
text_edge = np.zeros(self.Nsections)
for i in range(int(self.Nsections)):
res = minimize(self.solve_text_edge, 0.03, args = (C[i], thick[i], Mext_edge[i]), bounds = None, method = 'COBYLA', tol = 1.e-8)
for i in range(self.Nsections):
norm = self.solve_text_edge(0.01, C[i], thick[i], Mext_flap[i], 1.)
res = minimize(self.solve_text_edge, 0.01, args = (C[i], thick[i], Mext_edge[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
text_edge[i] = res['x']
# if not res['success']: print 'WARNING', res
if not res['success']: print 'WARNING solve_text_edge', i,res
# Calculating thickness flapwise fatigue
tfat_flap = np.zeros(self.Nsections)
for i in range(int(self.Nsections)):
res = minimize(self.solve_tfat_flap, 0.03, args = (C[i], thick[i], Mfat_flap[i]), bounds = None, method = 'COBYLA', tol = 1.e-8, options={'maxiter':10000})
for i in range(self.Nsections):
norm = self.solve_tfat_flap(0.01, C[i], thick[i], Mext_flap[i], 1.)
res = minimize(self.solve_tfat_flap, 0.01, args = (C[i], thick[i], Mfat_flap[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
tfat_flap[i] = res['x']
#print 'done', tfat_flap[i]
# if not res['success']: print 'WARNING', res
if not res['success']: print 'WARNING solve_tfat_flap', i, res
# Calculating thickness edgewise fatigue
tfat_edge = np.zeros(self.Nsections)
for i in range(int(self.Nsections)):
res = minimize(self.solve_tfat_edge, 0.03, args = (C[i], thick[i], Mfat_edge[i]), bounds = None, method = 'COBYLA', tol = 1.e-8)
for i in range(self.Nsections):
norm = self.solve_tfat_edge(0.01, C[i], thick[i], Mext_flap[i], 1.)
res = minimize(self.solve_tfat_edge, 0.01, args = (C[i], thick[i], Mfat_edge[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
tfat_edge[i] = res['x']
# if not res['success']: print 'WARNING', res
# Printing the 4 thicknesses
print 'text_flap\n', text_flap
print
print 'text_edge\n', text_edge
print
print 'tfat_flap\n', tfat_flap
print
print 'tfat_edge\n', tfat_edge
print
self.x = np.array(range(1,22))
if not res['success']: print 'WARNING solve_tfat_edge', i,res
self.text_flap = text_flap
self.text_edge = text_edge
self.tfat_flap = tfat_flap
......@@ -165,11 +157,8 @@ class SEAMRotor(Component):
# Calculating the final thickness for flap and edge direction
tfinal_flap = np.maximum(text_flap, tfat_flap)
tfinal_edge = np.maximum(text_edge, tfat_edge)
print 'tfinal_flap', tfinal_flap
print
print 'tfinal_edge', tfinal_edge
print
self.tfinal_flap = tfat_flap
self.tfinal_edge = tfat_edge
# Calculating different costs and mass
for i in range(0,int(self.Nsections)):
if i>0:
......@@ -177,28 +166,28 @@ class SEAMRotor(Component):
+2*self.sc_frac_edge*thick[i]*tfinal_edge[i])+(2*self.sc_frac_flap*C[i-1]*tfinal_flap[i-1]
+2*self.sc_frac_edge*thick[i-1]*tfinal_edge[i-1]))/2.
self.BladeWeight = self.BladeDens*volumen
self.BladeWeight = self.BladeDens*volumen
print 'volumen', volumen
print 'BladeWeigth', self.BladeWeight
rotor = self.BladeWeight*self.BladeCostPerMass/1e6 # Meuro
print 'rotor', rotor
# Scaling laws for hub, spinner, pitch
HubMass = 0.954*self.BladeWeight+5680.3
print 'HubMass', HubMass
HubCost = HubMass*self.HubCostPerMass/1e6
print 'HubCost', HubCost
#PitchCost
SpinnerMass = 18.5*self.Diameter-520.5
print 'SpinnerMass', SpinnerMass
SpinnerCost = SpinnerMass*self.SpinnerCostPerMass/1e6
print 'SpinnerCost', SpinnerCost
#
# rotor = self.BladeWeight*self.BladeCostPerMass/1e6 # Meuro
# print 'rotor', rotor
# # Scaling laws for hub, spinner, pitch
# HubMass = 0.954*self.BladeWeight+5680.3
# print 'HubMass', HubMass
#
# HubCost = HubMass*self.HubCostPerMass/1e6
# print 'HubCost', HubCost
#
# #PitchCost
#
# SpinnerMass = 18.5*self.rotor_diameter-520.5
# print 'SpinnerMass', SpinnerMass
#
# SpinnerCost = SpinnerMass*self.SpinnerCostPerMass/1e6
# print 'SpinnerCost', SpinnerCost
#Rotor = 3*(self.BladeWeight*self.BladeCostPerMass/1e6)+HubCost+PitchCost+SpinnerCost
......@@ -207,32 +196,32 @@ class SEAMRotor(Component):
# Defining functions in order to solve for the thickness
# Solving for t in flap direction, extremes
def solve_text_flap(self, t, C, thick, Mext_flap):
def solve_text_flap(self, t, C, thick, Mext_flap, norm):
Ine = (2./3.)*self.sc_frac_flap*C*t**3-self.sc_frac_flap*C*thick*t**2+self.sc_frac_flap*C*thick**2/2.*t
W = Ine/(thick/2.)
sext = self.SF_blade*1e3*Mext_flap/W/1e6
return abs(sext - self.Slim_ext_blade)
# print 'text_flap', abs(sext - self.Slim_ext_blade)
return abs(sext - self.Slim_ext_blade) / norm
#Solving for t in edge direction, extremes
def solve_text_edge(self, t, C, thick, Mext_edge):
def solve_text_edge(self, t, C, thick, Mext_edge, norm):
Ine = (2/3.)*self.sc_frac_edge*thick*t**3-self.sc_frac_edge*thick*C*t**2+self.sc_frac_edge*thick*C**2/2.*t
W = Ine/(C/2.);
sext = self.SF_blade*1.e3*Mext_edge/W/1.e6
return abs(sext - self.Slim_ext_blade)
return abs(sext - self.Slim_ext_blade) / norm
# Solving for t in flap direction, fatigue
def solve_tfat_flap(self, t, C, thick, Mfat_flap):
def solve_tfat_flap(self, t, C, thick, Mfat_flap, norm):
Ine = (2/3.)*self.sc_frac_flap*C*t**3-self.sc_frac_flap*C*thick*t**2+self.sc_frac_flap*C*thick**2/2.*t
W = Ine/(thick/2.)
sfat = np.maximum(1.e-6, self.SF_blade*1.e3*Mfat_flap/ W / 1.e6)
PM = self.Neq/(pow(10, (self.Slim_fat_blade - self.WohlerExpFlap*np.log10(sfat))))
return abs(PM - self.PMtarget)
return abs(PM - self.PMtarget) / norm
# Solving for t in edge direction, fatigue
def solve_tfat_edge(self, t, C, thick, Mfat_edge):
def solve_tfat_edge(self, t, C, thick, Mfat_edge, norm):
Ine = (2/3.)*self.sc_frac_edge*thick*t**3-self.sc_frac_edge*thick*C*t**2+self.sc_frac_edge*thick*C**2/2.*t
W = Ine/(C/2.)
sfat = np.maximum(1.e-6, self.SF_blade*1.e3*Mfat_edge/W/1.e6)
PM = self.Neq/(pow(10, (self.Slim_fat_blade - self.WohlerExpFlap*np.log10(sfat))))
return abs(PM - self.PMtarget)
return abs(PM - self.PMtarget) / norm
import numpy as np
import unittest
from SEAMRotor.SEAMRotor import SEAMRotor
from SEAMRotor.SEAMRotor import SEAMBladeStructure
class SEAMRotorTestCase(unittest.TestCase):
......@@ -20,13 +20,13 @@ class SEAMRotorTestCase(unittest.TestCase):
if __name__ == "__main__":
# unittest.main()
top = SEAMRotor()
top = SEAMBladeStructure()
top.Nsections = 21.
top.Neq = 1e7
top.WohlerExpFlap = 10.0
top.PMtarget = 1.0
top.Diameter = 177. #[m]
top.rotor_diameter = 177. #[m]
# top.RootChord = 2.0 #[m]
# top.MaxChord = 2.5 #[m]
top.MaxChordrR = 0.2 #[m]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment