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
gpkit.init_ipynb_printing()
Next declare constants (including take-off and cruising speed, which will be 'swept' over a range of values):
constants = {
'pi': (np.pi, "half of the circle constant"),
'CDA0': (0.03062702, "m^2", "fuselage drag area"),
'rho': (1.23, "kg m^3", "density of air"),
'mu': (1.78e-5, "kg/m*s", "viscosity of air"),
'S_wetratio': (2.05, "wetted area ratio"),
'k': (1.2, "form factor"),
'e': (0.96, "Oswald efficiency factor"),
'W_0': (4940, "N", "aircraft weight excluding wing"),
'N_ult': (2.5, "ultimate load factor"),
'tau': (0.12, "airfoil thickness to chord ratio"),
'C_Lmax': (2.0, "max CL with flaps down"),
'V': ('sweep', np.linspace(45, 55, 11), "m/s", "cruising speed"),
'V_min': ('sweep', np.linspace(20, 25, 11), "m/s", "takeoff speed"),
}
gpkit.monify_up(globals(), constants)
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):
free_variables = {
'A': "aspect ratio",
'S': ["m^2", "wing area"],
'C_D': "wing drag coefficient",
'C_L': "wing lift coefficent",
'C_f': "skin friction coefficient",
'Re': "Reynold's number",
'W': ["N", "aircraft weight"],
'W_w': ["N", "wing weight"],
'D': ["N", "cruising drag"],
}
gpkit.monify_up(globals(), free_variables)
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 model for overall drag
W_w_strc = 8.71e-5*(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*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 declare our objective (minimizing cruising drag) and write the rest of the system's equations:
gp = gpkit.GP( # Minimize
D, # cruising drag
[ # subject to
D == 0.5*rho*S*C_D*V**2, # the definition of drag,
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,
] + equations, substitutions=constants) # and the equations and constants we defined earlier.
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' Sweeping 2 variables over 121 passes Solving took 2.13 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:
gp.plot_frontiers(['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!
gp.plot_frontiers(['W_w', 'S', 'A'], 1, 3, (15,5))
This example was written by Edward Burnell (eburn@mit.edu), who welcomes any comments or suggestions.