Minimize fuel needed for a plane that can sprint and land quickly.
First we'll to import GPkit and turn on $\LaTeX$ printing for GPkit variables and equations. Then we'll create a few helper functions, though there should really be a better way to do this.
In [1]:
import numpy as np
from gpkit import Model, Variable, units
In [2]:
class Landing(Model):
def setup(self):
S = Variable("S", "m^2", "Wing area")
W_mto = Variable("W_{mto}", "N", "Maximum takeoff weight")
rho_sl = Variable("\\rho_{sl}", 1.23, "kg/m^3", "Air density, sea level")
C_Lmax = Variable("C_{L,max}", 1.5, "-", "Maximum C_L, flaps down")
V_stallmax = Variable("V_{stall,max}", 40, "m/s", "Stall speed")
V_stall = Variable("V_{stall}", "m/s", "Stall speed")
constraints = [W_mto <= 0.5*rho_sl*V_stall**2*C_Lmax*S,
V_stall <= V_stallmax]
cost = 1/W_mto
return cost, constraints
def test(self):
self.substitutions.update({"S": 10})
self.solve(verbosity=0)
Landing().test()
In [3]:
class Cruising(Model):
def setup(self):
rho = Variable("\\rho", 0.91, "kg/m^3", "Air density, 3000m")
mu = Variable("\\mu", 1.69e-5, "kg/m/s", "Dynamic viscosity, 3000m")
A = Variable("A", "-", "Aspect Ratio")
S = Variable("S", "m^2", "Wing area")
V = Variable("V", "m/s", "Flight speed")
C_L = Variable("C_L", "-", "Wing lift coefficent")
C_D = Variable("C_D", "-", "Wing drag coefficent")
T = Variable("T", "N", "Thrust force")
Re = Variable("Re", "-", "Reynold's number")
W = Variable("W", "N", "Aircraft weight")
constraints = [W == 0.5*rho*C_L*S*V**2,
T >= 0.5*rho*C_D*S*V**2,
Re == (rho/mu)*V*(S/A)**0.5]
cost = T/W
return cost, constraints
def test(self):
self.substitutions.update({"S":10, "C_L":1, "C_D":1, "V":10, "A":1})
self.solve(verbosity=0)
Cruising().test()
In [4]:
class Drag(Model):
def setup(self):
pi = Variable("\\pi", np.pi, "-", "Half of the circle constant")
e = Variable("e", 0.95, "-", "Wing spanwise efficiency")
S = Variable("S", "m^2", "Wing area")
A = Variable("A", "-", "Aspect Ratio")
tau = Variable("\\tau", "-")
C_L = Variable("C_L", "-", "Wing lift coefficent")
C_D = Variable("C_D", "-", "Wing drag coefficent")
C_Dp = Variable("C_{D_p}", "-", "fit to xrotor drag")
Re = Variable("Re", "-", "Reynold's number")
constraints = (C_D >= (0.05/S)*units.m**2 + C_Dp + C_L**2/(pi*e*A),
1 >= (2.56*C_L**5.88/(Re**1.54*tau**3.32*C_Dp**2.62) +
3.8e-9*tau**6.23/(C_L**0.92*Re**1.38*C_Dp**9.57) +
2.2e-3*Re**0.14*tau**0.033/(C_L**0.01*C_Dp**0.73) +
1.19e4*C_L**9.78*tau**1.76/(Re*C_Dp**0.91) +
6.14e-6*C_L**6.53/(Re**0.99*tau**0.52*C_Dp**5.19)))
cost = C_D.prod()
return cost, constraints
def test(self):
self.substitutions.update({"C_D":10, "C_L":1, "A":1,
"\\tau": 1, "Re":1e4})
self.solve(verbosity=0)
Drag().test()
In [5]:
class Engine(Model):
def setup(self):
T = Variable("T", "N", "Thrust force")
rho = Variable("\\rho", 0.91, "kg/m^3", "Air density, 3000m")
V = Variable("V", "m/s", "Flight speed")
P = Variable("P", "W", "engine power")
W_eng = Variable("W_{eng}", "N")
A_prop = Variable("A_{prop}", 0.785, "m^2", "Propeller disk area")
eta_eng = Variable("\\eta_{eng}", 0.35, "-", "Engine efficiency")
eta_v = Variable("\\eta_v", 0.85, "-", "Propeller viscous efficiency")
eta_i = Variable("\\eta_i", "-", "Aircraft efficiency")
eta_prop = Variable("\\eta_{prop}", "-", "Propeller efficiency")
eta_0 = Variable("\\eta_0", "-")
constraints = (eta_0 <= eta_eng*eta_prop,
eta_prop <= eta_i*eta_v,
4*eta_i + T*eta_i**2/(0.5*rho*V**2*A_prop) <= 4,
P >= T*V/eta_0,
W_eng >= 0.0372*P**0.8083 * units('N/W^0.8083'))
cost = W_eng
return cost, constraints
def test(self):
self.substitutions.update({"T":1e3, "V":20})
self.solve(verbosity=0)
Engine().test()
In [6]:
class FuelBurn(Model):
def setup(self):
g = Variable("g", 9.8, "m/s^2", "Gravitational constant")
eta_0 = Variable("\\eta_0", 0.5, "-")
W = Variable("W", "N", "Aircraft weight")
T = Variable("T", "N", "Thrust force")
R_min = Variable("R_{min}", 1e6, "m", "Minimum airplane range")
h_fuel = Variable("h_{fuel}", 42e6, "J/kg", "fuel heating value")
R = Variable("R", "m", "Airplane range")
z_bre = Variable("z_{bre}", "-")
W_fuel = Variable("W_{fuel}", "N", "Fuel weight")
# 4th order taylor approximation for e^x
z_bre_sum = 0
for i in range(1,5):
z_bre_sum += z_bre**i/np.math.factorial(i)
constraints = (R >= R_min,
z_bre >= g*R*T/(h_fuel*eta_0*W),
W_fuel/W >= z_bre_sum)
cost = W_fuel
return cost, constraints
def test(self):
self.substitutions.update({"T":1e3, "V":20, "W":1e3})
self.solve(verbosity=0)
FuelBurn().test()
In [7]:
class WingBox(Model):
def setup(self):
g = Variable("g", 9.8, "m/s^2", "Gravitational constant")
tau = Variable("\\tau", "-")
A = Variable("A", "-", "Aspect Ratio")
S = Variable("S", "m^2", "Wing area")
W = Variable("W", "N", "Aircraft weight")
W_wing = Variable("W_{wing}", "N")
W_tw = Variable("W_{tw}", "N")
f_wadd = Variable("f_{wadd}", 2, "-", "Wing added weight fraction")
N_lift = Variable("N_{lift}", 6.0, "-", "Wing loading multiplier")
sigma_max = Variable("\\sigma_{max}", 250e6, "Pa", "Allowable stress, 6061-T6")
sigma_maxshear = Variable("\\sigma_{max,shear}", 167e6, "Pa", "Allowable shear stress")
w = Variable("w", 0.5, "-", "Wing-box width/chord")
r_h = Variable("r_h", 0.75, "-", "Wing strut taper parameter")
I_cap = Variable("I_{cap}", "m^4", "Area moment of inertia per unit chord")
M_rbar = Variable("\\bar{M}_r", "-")
rho_alum = Variable("\\rho_{alum}", 2700, "kg/m^3", "Density of aluminum")
nu = Variable("\\nu", "-")
p = Variable("p", "-")
q = Variable("q", "-")
t_cap = Variable("t_{cap}", "-")
t_web = Variable("t_{web}", "-")
W_cap = Variable("W_{cap}", "N")
W_web = Variable("W_{web}", "N")
constraints = (2*q >= 1 + p,
p >= 2.2,
tau <= 0.25,
M_rbar >= W_tw*A*p/(24*units.N),
.92**2/2.*w*tau**2*t_cap >= I_cap * units.m**-4 + .92*w*tau*t_cap**2,
8 >= N_lift*M_rbar*A*q**2*tau/S/I_cap/sigma_max * units('Pa*m**6'),
12 >= A*W_tw*N_lift*q**2/tau/S/t_web/sigma_maxshear,
nu**3.94 >= .86*p**-2.38 + .14*p**0.56,
W_cap >= 8*rho_alum*g*w*t_cap*S**1.5*nu/3/A**.5,
W_web >= 8*rho_alum*g*r_h*tau*t_web*S**1.5*nu/3/A**.5,
W_wing/f_wadd >= W_cap + W_web)
cost = W_wing
return cost, constraints
def test(self):
self.substitutions.update({"\\tau":0.25, "A":1, "S":10,
"W":1e3, "W_{tw}":5e2})
self.solve(verbosity=0)
WingBox().test()
In [8]:
class FlightSegment(Model):
def setup(self):
landing = Landing()
cruising = Cruising()
drag = Drag()
engine = Engine()
fuel = FuelBurn()
spar = WingBox()
W = cruising["W"]
W_tw = spar["W_{tw}"]
W_wing = spar["W_{wing}"]
W_mto = landing["W_{mto}"]
W_fuel = fuel["W_{fuel}"]
g = Variable("g", 9.8, "m/s^2", "Gravitational constant")
W_zfw = Variable("W_{zfw}", "N", "Zero fuel weight")
W_eng = Variable("W_{eng}", 2e3, "N")
m_pay = Variable("m_{pay}", 500, "kg")
W_fixed = Variable("W_{fixed}", 14.7e3, "N", "Fixed weight")
constraints = [W_tw >= W_fixed + g*m_pay + W_eng,
W_zfw >= W_tw + W_wing,
W >= W_zfw,
W_mto >= W + W_fuel]
model = Model(W_fuel, constraints)
for subm in [landing, cruising, drag, engine, fuel, spar]:
model = model & subm
return model
def test(self):
self.solve(verbosity=0)
print FlightSegment().solve(verbosity=0).table(["cost", "freevariables"])
In [9]:
class Mission(Model):
def setup(self):
there = FlightSegment()
back = FlightSegment()
sprint = FlightSegment()
self.flightmodelnames = [there.modelname, back.modelname, sprint.modelname]
V_sprint = Variable("V_{sprintreqt}", 150, "m/s", "sprint speed requirement")
model = Model(there["W_{fuel}"] + back["W_{fuel}"],
[there["W"] >= there["W_{zfw}"] + back["W_{fuel}"],
sprint["W"] == there["W"],
sprint["V"] >= V_sprint])
exvars = ["V", "C_L", "C_D", "C_{D_{fuse}}", "C_{D_p}", "C_{D_i}", "T",
"Re", "W", "\\eta_i", "\\eta_{prop}", "\\eta_0", "W_{fuel}", "z_{bre}"]
for fl in there, back, sprint:
fl.substitutions = {k: v for k, v in fl.substitutions.items()
if k.nomstr in exvars}
model = model.merge(fl, excluded=exvars)
return model
ap = Mission()
ap.solve(verbosity=0)
modelnames = [ap.modelname] + [ap.modelname+flname for flname in ap.flightmodelnames]
print ap.solution.table(["freevariables"], included_models=modelnames)
In [ ]:
from IPython.paths import get_ipython_dir
from IPython.core.display import HTML
import os
def css_styling():
"""Load default custom.css file from ipython profile"""
base = get_ipython_dir()
styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
return HTML(styles)
css_styling()