In [1]:
import numpy as np
from gpkit.shortcuts import *
import gpkit.interactive
%matplotlib inline
In [2]:
N_lift = Var("N_{lift}", 6.0, "-", "Wing loading multiplier")
pi = Var("\\pi", np.pi, "-", "Half of the circle constant")
sigma_max = Var("\\sigma_{max}", 250e6, "Pa", "Allowable stress, 6061-T6")
sigma_maxshear = Var("\\sigma_{max,shear}", 167e6, "Pa", "Allowable shear stress")
g = Var("g", 9.8, "m/s^2", "Gravitational constant")
w = Var("w", 0.5, "-", "Wing-box width/chord")
r_h = Var("r_h", 0.75, "-", "Wing strut taper parameter")
f_wadd = Var("f_{wadd}", 2, "-", "Wing added weight fraction")
W_fixed = Var("W_{fixed}", 14.7e3, "N", "Fixed weight")
C_Lmax = Var("C_{L,max}", 1.5, "-", "Maximum C_L, flaps down")
rho = Var("\\rho", 0.91, "kg/m^3", "Air density, 3000m")
rho_sl = Var("\\rho_{sl}", 1.23, "kg/m^3", "Air density, sea level")
rho_alum = Var("\\rho_{alum}", 2700, "kg/m^3", "Density of aluminum")
mu = Var("\\mu", 1.69e-5, "kg/m/s", "Dynamic viscosity, 3000m")
e = Var("e", 0.95, "-", "Wing spanwise efficiency")
A_prop = Var("A_{prop}", 0.785, "m^2", "Propeller disk area")
eta_eng = Var("\\eta_{eng}", 0.35, "-", "Engine efficiency")
eta_v = Var("\\eta_v", 0.85, "-", "Propeller viscous efficiency")
h_fuel = Var("h_{fuel}", 42e6, "J/kg", "fuel heating value")
V_sprint_reqt = Var("V_{sprintreqt}", 150, "m/s", "sprint speed requirement")
W_pay = Var("W_{pay}", 500*9.81, "N")
R_min = Var("R_{min}", 1e6, "m", "Minimum airplane range")
V_stallmax = Var("V_{stall,max}", 40, "m/s", "Stall speed")
# sweep variables
R_min = Var("R_{min}", 5e6, "m", "Minimum airplane range")
V_stallmax = Var("V_{stall,max}", 40, "m/s", "Stall speed")
In [3]:
V = Vec(3, "V", "m/s", "Flight speed")
C_L = Vec(3, "C_L", "-", "Wing lift coefficent")
C_D = Vec(3, "C_D", "-", "Wing drag coefficent")
C_Dfuse = Vec(3, "C_{D_{fuse}}", "-", "Fuselage drag coefficent")
C_Dp = Vec(3, "C_{D_p}", "-", "drag model parameter")
C_Di = Vec(3, "C_{D_i}", "-", "drag model parameter")
T = Vec(3, "T", "N", "Thrust force")
Re = Vec(3, "Re", "-", "Reynold's number")
W = Vec(3, "W", "N", "Aircraft weight")
eta_i = Vec(3, "\\eta_i", "-", "Aircraft efficiency")
eta_prop = Vec(3, "\\eta_{prop}", "-")
eta_0 = Vec(3, "\\eta_0", "-")
W_fuel = Vec(2, "W_{fuel}", "N", "Fuel weight")
z_bre = Vec(2, "z_{bre}", "-")
S = Var("S", "m^2", "Wing area")
R = Var("R", "m", "Airplane range")
A = Var("A", "-", "Aspect Ratio")
I_cap = Var("I_{cap}", "m^4", "Spar cap area moment of inertia per unit chord")
M_rbar = Var("\\bar{M}_r", "-")
P_max = Var("P_{max}", "W")
V_stall = Var("V_{stall}", "m/s")
nu = Var("\\nu", "-")
p = Var("p", "-")
q = Var("q", "-")
tau = Var("\\tau", "-")
t_cap = Var("t_{cap}", "-")
t_web = Var("t_{web}", "-")
W_cap = Var("W_{cap}", "N")
W_zfw = Var("W_{zfw}", "N", "Zero fuel weight")
W_eng = Var("W_{eng}", "N")
W_mto = Var("W_{mto}", "N", "Maximum takeoff weight")
W_pay = Var("W_{pay}", "N")
W_tw = Var("W_{tw}", "N")
W_web = Var("W_{web}", "N")
W_wing = Var("W_{wing}", "N")
Check the vector constraints:
In [4]:
W == 0.5*rho*C_L*S*V**2
Out[4]:
In [5]:
steady_level_flight = (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)
landing_fc = (W_mto <= 0.5*rho_sl*V_stall**2*C_Lmax*S,
V_stall <= V_stallmax)
sprint_fc = (P_max >= T[2]*V[2]/eta_0[2],
V[2] >= V_sprint_reqt)
drag_model = (C_D >= (0.05/S)*gpkit.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)))
propulsive_efficiency = (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)
# 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)
range_constraints = (R >= R_min,
z_bre >= g*R*T[:2]/(h_fuel*eta_0[:2]*W[:2]),
W_fuel/W[:2] >= z_bre_sum)
punits = gpkit.units.parse_expression('N/W^0.8083')
weight_relations = (W_pay >= 500*g*gpkit.units.kg,
W_tw >= W_fixed + W_pay + W_eng,
W_zfw >= W_tw + W_wing,
W_eng >= 0.0372*P_max**0.8083 * punits,
W_wing/f_wadd >= W_cap + W_web,
W[0] >= W_zfw + W_fuel[1],
W[1] >= W_zfw,
W_mto >= W[0] + W_fuel[0],
W[2] == W[0])
wunits = gpkit.units.m**-4
munits = gpkit.units.parse_expression('Pa*m**6')
wing_structural_model = (2*q >= 1 + p,
p >= 2.2,
tau <= 0.25,
M_rbar >= W_tw*A*p/(24*gpkit.units.N),
.92**2/2.*w*tau**2*t_cap >= I_cap * wunits + .92*w*tau*t_cap**2,
8 >= N_lift*M_rbar*A*q**2*tau/S/I_cap/sigma_max * munits,
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
)
In [6]:
eqns = (weight_relations + range_constraints + propulsive_efficiency
+ drag_model + steady_level_flight + landing_fc + sprint_fc + wing_structural_model)
m = gpkit.Model(W_fuel.sum(), eqns)
In [7]:
m.interact()
The "local model" is the power-law tangent to the Pareto frontier, gleaned from sensitivities.
In [8]:
m.solution["localmodel"]
Out[8]:
In [9]:
from gpkit.interactive.plotting import sensitivity_plot
_ = sensitivity_plot(m)
In [10]:
from string import Template
fuelupdate_js = Template("""
var W_eng = $W_eng,
lam = $lam
fuel.shearinner.scalex = 1-$tcap*10
fuel.shearinner.scaley = 1-$tweb*100
fuel.airfoil.scaley = $tau/0.13
fuel.fuse.scalex = $W_fus/24000
fuel.wing.scalex = $b/2/14
fuel.wing.scaley = $cr*1.21
""")
def fuelupdate_py(sol):
varstrs = ("p", "S", "A", "t_{cap}", "t_{web}", "w",
"\\tau", "W_{eng}", "W_{mto}", "W_{wing}")
p, S, A, t_cap, t_web, w, tau, W_eng, W_mto, W_wing = sol.getvars(*varstrs)
lam = 0.5*(p-1)
return fuelupdate_js.substitute(lam = lam,
b = (S*A)**0.5,
cr = 2/(1+lam)*(S/A)**0.5,
tcap = t_cap/tau,
tweb = t_web/w,
tau = tau,
W_eng = W_eng,
W_fus = W_mto - W_wing - W_eng)
fuelconstraint_js = """
fuel.engine1.scale = Math.pow(W_eng/3000, 2/3)
fuel.engine2.scale = Math.pow(W_eng/3000, 2/3)
fuel.engine1.y = 6*lam
fuel.engine2.y = 6*lam
fuel.wingrect.scaley = 1-lam
fuel.wingrect.y = -6 + 5*lam
fuel.wingtaper.scaley = lam
fuel.wingtaper.y = 5*lam
"""
In [11]:
gpkit.interactive.showcadtoon("fuel",
"position:absolute; height:0; right:0; top:24em;")
In [12]:
gpkit.interactive.ractorpy(m, fuelupdate_py,
{"V_{stall,max}": (20, 50, 0.3),
"R_{min}": (1e6, 1e7, 1e5),
"V_{sprintreqt}": (100, 200, 1)},
fuelconstraint_js)
In [13]:
gpkit.interactive.ractorjs("fuel", m, fuelupdate_py,
{"V_{stall,max}": (20, 50, 3),
"R_{min}": (1e6, 1e7, 2e6),
"V_{sprintreqt}": (100, 200, 20)},
fuelconstraint_js)
This concludes the aircraft example. Try playing around with the sliders up above until you're bored; then check out one of the other examples. Thanks for reading!
In [14]:
from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
"""Load default custom.css file from ipython profile"""
base = utils.path.get_ipython_dir()
csspath = os.path.join(base,'profile_default/static/custom/custom.css')
styles = "<style>\n%s\n</style>" % (open(csspath,'r').read())
return HTML(styles)
css_styling()
Out[14]: