Skip to content

Commit

Permalink
Minor update to GPfit
Browse files Browse the repository at this point in the history
  • Loading branch information
@courtin
courtin committed Apr 18, 2018
1 parent b49b7b1 commit 9b5e8b3
Show file tree
Hide file tree
Showing 4 changed files with 22 additions and 77 deletions.
29 changes: 7 additions & 22 deletions gpkitmodels/GP/aircraft/motor/motor_test.py
View file
@@ -1,6 +1,6 @@
from gpkit import Model, parse_variables, SignomialsEnabled, SignomialEquality, units
from motor import Propulsor, Motor, MotorPerf
from gpkitmodels.GP.aircraft.prop.propeller import Propeller, ActuatorProp, SimpleQProp, MultiElementProp
from gpkitmodels.GP.aircraft.prop.propeller import Propeller, ActuatorProp, MultiElementProp
from gpkitmodels.GP.aircraft.wing.wing_test import FlightState

class Propulsor_Test(Model):
Expand All @@ -12,13 +12,12 @@ def setup(self):
p = Propulsor()
pp = p.flight_model(p,fs)
pp.substitutions[pp.prop.T] = 100
#pp.substitutions[pp.prop.AR_b] = 15
self.cost = 1./pp.motor.etam + p.W/(1000*units('lbf')) + 1./pp.prop.eta

return fs,p,pp

class Actuator_Propulsor_Test(Model):
"""Propulsor Test Model
"""Propulsor Test Model w/ Actuator Disk Propeller
"""

def setup(self):
Expand All @@ -27,14 +26,13 @@ def setup(self):
p = Propulsor()
pp = p.flight_model(p,fs)
pp.substitutions[pp.prop.T] = 100
#pp.substitutions[pp.prop.AR_b] = 15
self.cost = pp.motor.Pelec/(1000*units('W')) + p.W/(1000*units('lbf'))

return fs,p,pp


class MultiElement_Propulsor_Test(Model):
"""Propulsor Test Model
"""Propulsor Test Model w/ Blade Element Propeller
"""

def setup(self):
Expand All @@ -43,33 +41,22 @@ def setup(self):
p = Propulsor()
pp = p.flight_model(p,fs)
pp.substitutions[pp.prop.T] = 100
#pp.substitutions[pp.prop.AR_b] = 15
self.cost = pp.motor.Pelec/(1000*units('W')) + p.W/(1000*units('lbf'))

return fs,p,pp

def actuator_propulsor_test():

test = Actuator_Propulsor_Test()
#sol = test.debug()
sol = test.solve()
#sol = test.solve()
#print sol.table()

test.solve()

def ME_propulsor_test():
print "Blade Element Model Test"
test = MultiElement_Propulsor_Test()
#sol = test.debug()
sol = test.localsolve()
#sol = test.solve()
print sol.table()
def propulsor_test():

def propulsor_test():
test = Propulsor_Test()
#sol = test.debug()
#sol = test.localsolve(iteration_limit = 400)
sol = test.solve()
#print sol.table()

class Motor_P_Test(Model):
def setup(self):
Expand Down Expand Up @@ -111,8 +98,7 @@ def setup(self):
def motor_test():
test = Motor_P_Test()
sol = test.solve()
#sol = test.debug()
#print sol.table()


def test():
motor_test()
Expand All @@ -122,4 +108,3 @@ def test():

if __name__ == "__main__":
test()
#speed_280_test()
4 changes: 2 additions & 2 deletions gpkitmodels/GP/aircraft/prop/arccos_fit.py
View file
Expand Up @@ -3,13 +3,13 @@
from numpy import arccos,arange
from gpfit.fit import fit
from numpy.random import random_sample
i = arange(0.0001,1,.001)
i = arange(0.0001,3,.001)
j = arccos(exp(-i))
#P = Vdd**2 + 30*Vdd*exp(-(Vth-0.06*Vdd)/0.039)
#u = vstack((Vdd,Vth))
x = log(i)
y = log(j)
K = 1
K = 3

cstrt, rmsErr = fit(x,y,K,"SMA")
print rmsErr
10 changes: 1 addition & 9 deletions gpkitmodels/GP/aircraft/prop/prop_test.py
View file
Expand Up @@ -15,31 +15,23 @@ def simpleprop_test():
[fs, p, pp])
m.substitutions.update({"rho": 1.225, "V": 50, "T": 100, "omega":1000})
sol = m.solve()
#print sol.table()


def ME_eta_test():

fs = FlightState()
Propeller.flight_model = MultiElementProp
p = Propeller()

#p.substitutions[p.R] = 2
pp = p.flight_model(p,fs)
pp.substitutions[pp.T] = 100

#pp.substitutions[pp.omega] = 500

#pp.substitutions[pp.omega] = 1000
pp.cost = 1./pp.eta + pp.Q/(100.*units("N*m")) + p.T_m/(1000*units('N'))
#sol = pp.debug()
sol = pp.localsolve(iteration_limit = 400)
print sol.table()


def test():
"tests"
#simpleprop_test()
simpleprop_test()
ME_eta_test()
if __name__ == "__main__":
test()
Expand Down
56 changes: 12 additions & 44 deletions gpkitmodels/GP/aircraft/prop/propeller.py
View file
Expand Up @@ -70,24 +70,19 @@ class Blade_Element_Performance(Model):
va [m/s] Axial induced velocity
vt [m/s] Tangential induced velocity
G [m^2/s] Circulation
cl [-] local lift coefficient
cl [-] local lift coefficient
cd [-] local drag coefficient
cdp [-] local drag coefficient
B 2 [-] number of blades
r [m] local radius
lam_w [-] advance ratio
eps [-] blade efficiency
AR_b [-] blade aspect ratio
AR_b [-] blade aspect ratio
AR_b_max 50 [-] max blade aspect ratio
Ut [m/s] tangential freestreem
Re [-] blade reynolds number
f [-] intermediate tip loss variable
F [-] Prandtl tip loss factor
cl_max .6 [-] max airfoil cl
cl_max .6 [-] max airfoil cl
dr [m] length of blade element
FF [-] form factor
t_c .1 [-] airfoil t/c
C_f [-] skin friction coefficient
M [-] Mach number
a 295 [m/s] Speed of sound at altitude
Expand All @@ -107,7 +102,7 @@ def setup(self,static, state):
constraints = [TCS([Wa>=V + va]),
TCS([Wt + vt<=omega*r]),
TCS([G == (1./2.)*Wr*c*cl]),
F == (2./pi)*(1.02232*f**.0458729)**(1./.1), #This is the GP fit of arccos(exp(-f))
F == (2./pi)*(1.01116*f**.0379556)**(1./.1), #This is the GP fit of arccos(exp(-f))
M == Wr/a,
lam_w == (r/R)*(Wa/Wt),
va == vt*(Wt/Wa),
Expand All @@ -118,10 +113,6 @@ def setup(self,static, state):
Re == Wr*c*rho/mu,
eta_i == (V/(omega*r))*(Wt/Wa),
TCS([f+(r/R)*B/(2*lam_w) <= (B/2.)*(1./lam_w)]),
#TCS([cd >= cl**2/(pi*AR_b*.85) + cdp]),
#TCS([FF >= 1. + .3*t_c + (100*t_c)**.4]),
#C_f == .027/(Re**(1./7.)),
#cdp >= FF*C_f*(2.+.4*t_c),

XfoilFit(fd, cd, [cl,Re], name="polar"),

Expand All @@ -146,16 +137,16 @@ class MultiElementProp(Model):
Variables
---------
Mtip .5 [-] tip mach number
Mtip .5 [-] Max tip mach number
omega_max 10000 [rpm] maximum rotation rate
eta [-] overall efficiency
omega [rpm] rotation rate
T [lbf] total thrust
Q [N*m] total torque
"""
def setup(self,static, state, N = 5):
exec parse_variables(MultiElementProp.__doc__)
T = self.T = Variable('T', 'lbf', 'Overall thrust')
Q = self.Q = Variable('Q', 'N*m', 'Overall torque')
eta = self.eta = Variable('eta', '-', 'Overall efficiency')
omega = self.omega = Variable('omega', 'rpm', 'rotation rate')
omega_max = Variable('omega_max',10000, 'rpm', 'rotation rate')
#Mtip = Variable('Mtip', .5, '-', 'Tip Mach')

with Vectorize(N):
blade = Blade_Element_Performance(static, state)

Expand All @@ -166,14 +157,12 @@ def setup(self,static, state, N = 5):


with SignomialsEnabled():
# constraints = [blade.r <= static.R]
for n in range(1,N):
constraints += [TCS([blade["r"][n] >= blade["r"][n-1] + static.R/N]),
blade["eta_i"][n] == blade["eta_i"][n-1],
#blade["c"][n] == static["c"][n],
]
constraints += [TCS([T <= sum(blade["dT"])]),
TCS([Q >= sum(blade["dQ"][n] for n in range(N))]),
TCS([Q >= sum(blade["dQ"])]),
eta == state.V*T/(omega*Q),
blade["M"][-1] <= Mtip,
static.T_m>=T,
Expand Down Expand Up @@ -205,24 +194,3 @@ def setup(self, N = 5):
exec parse_variables(Propeller.__doc__)
self.N = N
return [W >= K*T_m*R**2]



class Multi_Element_Propeller(Model):
""" Propeller Model
Variables
---------
R [ft] prop radius
W [lbf] prop weight
K 4e-4 [1/ft^2] prop weight scaling factor
T_m [lbf] prop max static thrust
"""

flight_model = MultiElementProp
def setup(self):
exec parse_variables(Multi_Element_Propeller.__doc__)

constraints = [W >= K*T_m*R**2]

return constraints

0 comments on commit 9b5e8b3

Please sign in to comment.