codykarcher
commented
8 years ago Another option returning a posynomial drag polar instead of a constant. Note, this one takes about twice as many iterations to converge:
import numpy as np
from gpkit import Variable, Model
from gpkit.constraints.set import ConstraintSet
import subprocess
import matplotlib.pyplot as plt
# Constants
k = Variable("k", 1.2, "-", "form factor")
e = Variable("e", 0.95, "-", "Oswald efficiency factor")
mu = Variable("\\mu", 1.78e-5, "kg/m/s", "viscosity of air")
pi = Variable("\\pi", np.pi, "-", "half of the circle constant")
rho = Variable("\\rho", 1.23, "kg/m^3", "density of air")
tau = Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
N_ult = Variable("N_{ult}", 3.8, "-", "ultimate load factor")
V_min = Variable("V_{min}", 22, "m/s", "takeoff speed")
C_Lmax = Variable("C_{L,max}", 1.5, "-", "max CL with flaps down")
S_wetratio = Variable("(\\frac{S}{S_{wet}})", 2.05, "-", "wetted area ratio")
W_W_coeff1 = Variable("W_{W_{coeff1}}", 8.71e-5, "1/m",
"Wing Weight Coefficent 1")
W_W_coeff2 = Variable("W_{W_{coeff2}}", 45.24, "Pa",
"Wing Weight Coefficent 2")
CDA0 = Variable("(CDA0)", 0.031, "m^2", "fuselage drag area")
W_0 = Variable("W_0", 4940.0, "N", "aircraft weight excluding wing")
# Free Variables
D = Variable("D", "N", "total drag force")
A = Variable("A", "-", "aspect ratio")
S = Variable("S", "m^2", "total wing area")
V = Variable("V", "m/s", "cruising speed")
W = Variable("W", "N", "total aircraft weight")
Re = Variable("Re", "-", "Reynold's number")
C_D = Variable("C_D", "-", "Drag coefficient of wing")
C_L = Variable("C_L", "-", "Lift coefficent of wing")
C_f = Variable("C_f", "-", "skin friction coefficient")
W_w = Variable("W_w", "N", "wing weight")
C_D_fuse = Variable("C_D_fuse", "-", "Fuselage drag coefficient")
C_D_p = Variable("C_D_p", "-", "Profile drag coefficient of wing")
C_D_ind = Variable("C_D_ind", "-", "Induced drag coefficient of wing")
constraints = []
# Drag model
constraints += [C_D_p >= 0.05,
C_D_fuse == CDA0/S,
C_D_ind == C_L**2/(pi*A*e),
C_D >= C_D_fuse + C_D_p + C_D_ind]
# Wing weight model
W_w_strc = W_W_coeff1*(N_ult*A**1.5*(W_0*W*S)**0.5)/tau
W_w_surf = W_W_coeff2 * S
constraints += [W_w >= W_w_surf + W_w_strc]
# and the rest of the models
constraints += [D >= 0.5*rho*S*C_D*V**2,
Re <= (rho/mu)*V*(S/A)**0.5,
W <= 0.5*rho*S*C_L*V**2,
W <= 0.5*rho*S*C_Lmax*V_min**2,
W >= W_0 + W_w]
class XFOIL():
def __init__(self):
self.pathname = "/Users/codykarcher/Xfoil/bin/./xfoil"
def cd_model(self, x0, max_iter=100):
# File XFOIL will save data to
polarfile = 'polars.txt'
# Remove any existing file. Prevents XFOIL Error
ls_res = subprocess.check_output(["ls -a"], shell = True)
ls = ls_res.split()
if polarfile in ls_res:
subprocess.call("rm " + polarfile, shell = True)
# Open XFOIL and run case
proc = subprocess.Popen([self.pathname],
stdout=subprocess.PIPE, stdin=subprocess.PIPE)
proc.stdin.write('naca2412 \n' +
'oper \n' +
"iter %d\n" %(max_iter)+
'visc \n' +
"%.0f \n" %(x0["Re"]) +
'pacc \n' +
'\n' +
'\n' +
"cl 0\n" +
"cl %.3f \n" %(0.5*x0["C_L"]) +
"cl %.3f \n" %(x0["C_L"]) +
'pwrt \n' +
polarfile +
'\n' +
'q')
stdout_val = proc.communicate()[0]
proc.stdin.close()
# Read data from XFOIL output
data = np.genfromtxt(polarfile, skip_header=12)
data2 = data[data[:, 0].argsort()]
alpha = data2[:, 0]
cl = data2[:, 1]
cd = data2[:, 2]
cdp = data2[:, 3]
cm = data2[:, 4]
xtr_top = data2[:, 5]
xtr_bot = data2[:, 6]
# Calculate coefficients for Drag polar
CD_0 = cd[0]
K = cd[-1]/cl[-1]**2
# Return constraint
return ConstraintSet([C_D_p >= CD_0 + K*C_L])
class SimpleCDModel(ConstraintSet):
def as_posyslt1(self):
raise ValueError("SimpleCDModel is not allowed to solve as a GP.")
# Returns constraint list with appended XFOIL constraint as a gp constraint set
def as_gpconstr(self, x0):
if x0:
xfoil=XFOIL()
self[0] = xfoil.cd_model(x0)
return super(SimpleCDModel, self).as_gpconstr(x0)
constraints = SimpleCDModel(constraints)
m = Model(D, constraints)
sol = m.localsolve('cvxopt', verbosity=1)Cost
----
490.4 [N]
Free Variables
--------------
A : 6.467 aspect ratio
C_D : 0.02994 Drag coefficient of wing
C_D_fuse : 0.002095 Fuselage drag coefficient
C_D_ind : 0.008425 Induced drag coefficient of wing
C_D_p : 0.01941 Profile drag coefficient of wing
C_L : 0.4032 Lift coefficent of wing
D : 490.4 [N] total drag force
Re : 1.423e+06 Reynold's number
S : 14.8 [m**2] total wing area
V : 42.43 [m/s] cruising speed
W : 6606 [N] total aircraft weight
W_w : 1666 [N] wing weight
Constants
---------
(CDA0) : 0.031 [m**2] fuselage drag area
C_{L,max} : 1.5 max CL with flaps down
N_{ult} : 3.8 ultimate load factor
V_{min} : 22 [m/s] takeoff speed
W_0 : 4940 [N] aircraft weight excluding wing
W_{W_{coeff1}} : 8.71e-05 [1/m] Wing Weight Coefficent 1
W_{W_{coeff2}} : 45.24 [Pa] Wing Weight Coefficent 2
\mu : 1.78e-05 [kg/m/s] viscosity of air
\pi : 3.142 half of the circle constant
\rho : 1.23 [kg/m**3] density of air
\tau : 0.12 airfoil thickness to chord ratio
e : 0.95 Oswald efficiency factor
Sensitivities
-------------
W_0 : 1.024 aircraft weight excluding wing
W_{W_{coeff1}} : 0.1876 Wing Weight Coefficent 1
N_{ult} : 0.1876 ultimate load factor
W_{W_{coeff2}} : 0.126 Wing Weight Coefficent 2
(CDA0) : 0.06999 fuselage drag area
C_{L,max} : -0.1498 max CL with flaps down
\rho : -0.1498 density of air
\tau : -0.1876 airfoil thickness to chord ratio
e : -0.2814 Oswald efficiency factor
\pi : -0.2814 half of the circle constant
V_{min} : -0.2997 takeoff speed I'm also rather surprised that the two results are so different (~40%).
bqpd
whoburg
I'm planning to add this as an example to the docs:
any issues or suggested improvements?