Based on "Geometric Programming for Aircraft Design Optimization" by Hoburg and Abbeel, and using GPkit for model representation and solver interfacing.
An airplane with that consumes as little fuel as possible when cruising (while remaining capable of taking off from the ground, with wings that won't break, etc.)
Start by importing Numpy and GPkit, and turning on LaTex printing for GPkit variables and equations:
import numpy as np
import gpkit
import gpkit.interactive
Next declare constants (including take-off and cruising speed, which will be 'swept' over a range of values):
k = gpkit.Variable("k", 1.2, "-", "form factor")
e = gpkit.Variable("e", 0.95, "-", "Oswald efficiency factor")
V = gpkit.Variable("V", ("sweep", np.linspace(45, 55, 5)), "m/s", "cruising speed")
pi = gpkit.Variable("\\pi", np.pi, "-", "half of the circle constant")
mu = gpkit.Variable("\\mu", 1.78e-5, "kg/m/s", "viscosity of air")
rho = gpkit.Variable("\\rho", 1.23, "kg/m^3", "density of air")
tau = gpkit.Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
W_0 = gpkit.Variable("W_0", 4940.0, "N", "aircraft weight excluding wing")
CDA0 = gpkit.Variable("(CDA0)", 0.031, "m^2", "fuselage drag area")
N_ult = gpkit.Variable("N_{ult}", 3.8, "-", "ultimate load factor")
V_min = gpkit.Variable("V_{min}", ("sweep", np.linspace(20, 25, 5)), "m/s", "takeoff speed")
C_Lmax = gpkit.Variable("C_{L,max}", 1.5, "-", "max CL with flaps down")
S_wetratio = gpkit.Variable("(\\frac{S}{S_{wet}})", 2.05, "-", "wetted area ratio")
Then declare the rest of the variables (there's no need to separate variables and constants like this, I just find it a little more readable):
D = gpkit.Variable("D", "N", "total drag force")
A = gpkit.Variable("A", "-", "aspect ratio")
S = gpkit.Variable("S", "m^2", "total wing area")
W = gpkit.Variable("W", "N", "total aircraft weight")
Re = gpkit.Variable("Re", "-", "Reynold's number")
W_w = gpkit.Variable("W_w", "N", "wing weight")
C_D = gpkit.Variable("C_D", "-", "Drag coefficient of wing")
C_L = gpkit.Variable("C_L", "-", "Lift coefficent of wing")
C_f = gpkit.Variable("C_f", "-", "skin friction coefficient")
LaTex makes it easier to check one's equations (note the use of == to form an equality constraint):
D == 0.5*rho*S*C_D*V**2
The drag and wing-weight models consist of several sub-models, so we'll write them here for convenience:
equations = []
C_D_fuse = CDA0/S # Fuselage viscous drag
C_D_wpar = k*C_f*S_wetratio # wing parasitic drag
C_D_ind = C_L**2/(pi*A*e) # and induced drag
equations += [C_D >= C_D_fuse + C_D_wpar + C_D_ind] # together make a drag model
W_w_strc = (8.71e-5*1/gpkit.units.m)*(N_ult*A**1.5*(W_0*W*S)**0.5)/tau # A model for the wing's structural members
W_w_surf = (45.24*gpkit.units.N/gpkit.units.m**2)*S # and the weight of the wing's 'skin'
equations += [W_w >= W_w_surf + W_w_strc] # form a model of total wing weight
Then we'll write the rest of the system's equations:
equations += [D >= 0.5*rho*S*C_D*V**2, # the definition of C_D
Re <= (rho/mu)*V*(S/A)**0.5, # the definition of Reynold's number
C_f >= 0.074/Re**0.2, # a skin friction model
W <= 0.5*rho*S*C_L*V**2, # a cruising lift constraint
W <= 0.5*rho*S*C_Lmax*V_min**2, # a takeoff lift constraint
W >= W_0 + W_w] # the definition of aircraft weight
gp = gpkit.GP(D, # Minimize cruising drag
equations) # subject to our model
We can check that we've got the right equations by examining the LaTex representation of the system:
gp
data = gp.solve()
Using solver 'mosek' Solving for 11 variables over 25 passes. Solving took 1.04 seconds.
GPkit uses matplotlib to visualize the optimal ('Pareto') frontier of a system.
We'll start up the iPython Notebook plotting environment so we can use them:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
Plotting any particular design variables is then quite straightforward:
gpkit.interactive.plotting.plot_frontiers(gp, [D, C_D, C_L], 1, 3, (15,5))
We can see that drag does not increase with cruising speed as drastically as might have been expected.
This is because (under the assumptions of our equations) increasing speed allows one to fly with a smaller and lighter wing.
But don't believe me, take a look at the graphs!
gpkit.interactive.plotting.plot_frontiers(gp, [W_w, S, A], 1, 3, (15,5))
This example was written by Edward Burnell (eburn@mit.edu), who welcomes any comments or suggestions.