Aircraft wing design¶
$$ \newcommand{\half}{{\frac{1}{2}}} \newcommand{\AR}{{\text{AR}}} \newcommand{\Re}{{\text{Re}}} \newcommand{\FDA}{{(C_DS)_\text{fuselage}}} \newcommand{\argmin}{\operatornamewithlimits{argmin}} \newcommand{\minimize}{{\text{minimize}}} $$ ASE7030: Convex Optimization, Inha University.
Jong-Han Kim (jonghank@inha.ac.kr)
This section presents a simple wing design example, where the challenge is to size a wing with total area $S$, span $b$, and aspect ratio $\AR=b^2/S$. These parameters should be chosen to minimize the total drag,
\begin{align*} D=\half \rho V^2C_DS \end{align*}The drag coefficient is modeled as the sum of fuselage parasite drag, wing parasite drag, and induced drag,
\begin{align*} C_D = \frac{\FDA}{S} + kC_f \frac{S_{\text{wet}}}{S} + \frac{C_L^2}{\pi e \AR } \end{align*}where $\FDA$ is the fuselage drag area, $k$ is a form factor that accounts for pressure drag, $S_\text{wet}/S$ is the wetted area ratio, and $e$ is the Oswald efficiency factor.
For a fully turbulent boundary layer, the skin friction coefficient $C_f$ can be approximated as
\begin{align*} C_f = \frac{0.074}{\Re^{0.2}} \end{align*}where $\Re = \rho V c/\mu$ is the Reynolds number at the mean chord $c=\sqrt{S/\AR}$. The total aircraft weight $W$ is modeled as the sum of a fixed weight $W_0$ and the wing weight,
\begin{align*} W = W_0 + W_w \end{align*}The wing weight is modeled as
\begin{align*} W_w & = 45.42S + 8.71\times 10^{-5}\frac{N_\text{lift}b^3\sqrt{W_0 W}}{S\tau} \\ & = 45.42S + 8.71\times 10^{-5}\frac{N_\text{lift}\AR^{3/2}\sqrt{W_0 WS}}{\tau} \end{align*}where $N_\text{lift}$ is the ultimate load factor for structural sizing, and $\tau$ is the airfoil thickness-to-chord ratio.
The weight equations are coupled to the drag equations by the constraint that lift equals weight,
\begin{align*} W = \half \rho V^2 C_L S \end{align*}Finally, for safe landing, the aircraft should be capable of flying at a reduced speed $V_\text{min}$, subject to a stall constraint,
\begin{align*} W \le \half \rho V_\min^2 C_{L,\max} S \end{align*}We must choose values of $S$, $\AR$, and $V$ that minimize drag, subject to all the relations in the preceding text. The constant parameters are given below.
\begin{align*} \begin{array}{c|c|c} \text{Quantity} & \text{Value} & \text{Description} \\ \hline \hline \FDA & 0.0306 \text{m}^2 & \text{Fuselage drag area} \\ \hline\hline \rho & 1.23 \text{kg}/\text{m}^3 & \text{Air density} \\ \hline \mu & 1.78\times 10^{-5}\text{kg/m/s} & \text{Air viscosity} \\ \hline S_\text{wet}/S & 2.05 & \text{Wetted area ratio} \\ \hline k & 1.2 & \text{Form factor} \\ \hline e & 0.96 & \text{Oswald efficiency factor} \\ \hline W_0 & 4940 \text{N} & \text{Aircraft weight excluding wing} \\ \hline N_\text{lift} & 2.5 & \text{Ultimate load factor} \\ \hline \tau & 0.12 & \text{Airfoil thickness-to-chord ratio} \\ \hline V_\min & 22 \text{m/s} & \text{Desired landing speed} \\ \hline C_{L,\max} & 2.0 & \text{Maximum $C_L$ with flaps down} \end{array} \end{align*}In [1]:
FDA = 0.0306
rho = 1.23
mu = 1.78e-5
SwetS = 2.05
k = 1.2
e = 0.96
W0 = 4940
Nlift = 2.5
tau = 0.12
Vmin = 22
CLmax = 2.0
The problem can be formulated as follows.
\begin{align*} \underset{\AR,C_D,C_f,C_L,\Re,S,V,W,W_w}{\minimize} \quad & \half \rho V^2 C_D S \\ \text{subject to} \quad & C_f = \frac{0.074}{\Re^{0.2}} \\ & C_D = \frac{\FDA}{S} + kC_f \frac{S_{\text{wet}}}{S} + \frac{C_L^2}{\pi e \AR } \\ & W = \frac{1}{2}\rho V^2 C_L S \\ & W = W_0 + W_w \\ & W_w = 45.42S + 8.71\times 10^{-5}\frac{N_\text{lift}\AR^{3/2}\sqrt{W_o WS}}{\tau}\\ & W \le \half \rho V_\min^2 C_{L,\max} S \\ & \Re = \frac{\rho V}{\mu} \sqrt{\frac{S}{\AR}} \end{align*}Relaxing the posynomial equalities (2nd, 4th, 5th) leads to the following geometric programming.
\begin{align*} \underset{\AR,C_D,C_f,C_L,\Re,S,V,W,W_w}{\minimize} \quad & V^2 C_D S \\ \text{subject to} \quad & \frac{0.074}{C_f\Re^{0.2}} = 1\\ & \frac{\FDA}{C_DS} + \frac{kC_f}{C_D} \frac{S_{\text{wet}}}{S} + \frac{C_L^2}{\pi e C_D\AR } \le 1 \\ & \frac{2W}{\rho V^2 C_L S} = 1\\ & \frac{W_0}{W} + \frac{W_w}{W} \le 1\\ & 45.42\frac{S}{W_w} + 8.71\times 10^{-5}\frac{N_\text{lift}\AR^{3/2}\sqrt{W_o WS}}{{W_w}\tau} \le 1\\ & \frac{2W}{\rho V_\min^2 SC_{L,\max}} \le 1 \\ & \frac{\rho V}{\mu \Re} \sqrt{\frac{S}{\AR}} = 1 \end{align*}The change of variables by
\begin{align*} a &= \log \AR \\ s &= \log S \\ c_d & = \log C_D \\ c_l &= \log C_L \\ c_f &= \log C_f \\ r &= \log \Re \\ w &= \log W \\ w_w &= \log W_w \\ v &= \log V \end{align*}and taking the logarithms of the cost and the constraint, we have the following equivalent convex problem.
\begin{align*} \underset{a,c_d,c_f,c_l,r,s,v,w,w_w}{\minimize} \quad & s + c_d + 2v \\ \text{subject to} \quad & c_f +0.2r - \log 0.074 = 0 \\ & w - 2v - c_l -s + \log \frac{2}{\rho} = 0 \\ & v - r + 0.5s - 0.5a + \log\frac{\rho}{\mu} = 0 \\ & w - s + \log \frac{2}{\rho V_\min^2C_{L,\max} } \le 0 \\ & \log \left\{ \exp\left( -c_d -s + \log\FDA \right) + \exp\left( c_f - c_d + \log\left(k\frac{S_\text{wet}}{S}\right) \right) + \exp\left( 2c_l - c_d - a - \log\left(\pi e\right) \right) \right\} \le 0 \\ & \log \{ \exp\left( \log W_0 - w\right) + \exp\left( w_w - w\right) \} \le 0 \\ & \log \left\{ \exp\left(s-w_w + \log 45.42\right) + \exp\left(1.5a+ 0.5w + 0.5s-w_w+ \log \frac{8.71\times 10^{-5}N_\text{lift}\sqrt{W_0}}{\tau}\right)\right\} \le 0 \\ \end{align*}Note that all the equality constraint functions are affine and all the inequality constraint functions are convex.
In [2]:
import numpy as np
# minimize c^T x
# subject to Ax == b
# Cx <= d
# log_sum_exp(Ex+f) <= 0
# log_sum_exp(Gx+h) <= 0
# log_sum_exp(Px+q) <= 0
import cvxpy as cp
# x: [ a, c_d, c_f, c_l, r, s, v, w, w_w]
x = cp.Variable(9)
# min c^T x
# a, c_d, c_f, c_l, r, s, v, w, w_w#
c = np.array([ 0, 1, 0, 0, 0, 1, 2, 0, 0]) # s +c_d +2v
obj = cp.Minimize(c.T@x)
# Ax == b
# a, c_d, c_f, c_l, r, s, v, w, w_w#
A = np.array([[ 0, 0, 1, 0, .2, 0, 0, 0, 0], # c_f +.2r
[ 0, 0, 0, -1, 0, -1, -2, 1, 0], # w -2v -c_l -s
[-.5, 0, 0, 0, -1, .5, 1, 0, 0]]) # v -r +.5s -.5a
b = -np.array([-np.log(0.074), np.log(2/rho), np.log(rho/mu)])
# Cx <= d
# a, c_d, c_f, c_l, r, s, v, w, w_w#
C = np.array([ 0, 0, 0, 0, 0, -1, 0, 1, 0]) # w -s
d = -np.log(2/rho/Vmin/Vmin/CLmax)
# log_sum_exp(Ex+f) <= 0
# a, c_d, c_f, c_l, r, s, v, w, w_w#
E = np.array([[ 0, -1, 0, 0, 0, -1, 0, 0, 0], # -c_d -s
[ 0, -1, 1, 0, 0, 0, 0, 0, 0], # c_f -c_d
[ -1, -1, 0, 2, 0, 0, 0, 0, 0]]) # 2c_l -c_d -a
f = np.array([np.log(FDA), np.log(k*SwetS), -np.log(np.pi*e)])
# log_sum_exp(Gx+h) <= 0
# a, c_d, c_f, c_l, r, s, v, w, w_w#
G = np.array([[ 0, 0, 0, 0, 0, 0, 0, -1, 0], # -w
[ 0, 0, 0, 0, 0, 0, 0, -1, 1]]) # w_w -w
h = np.array([np.log(W0), 0])
# log_sum_exp(Px+q) <= 0
# a, c_d, c_f, c_l, r, s, v, w, w_w#
P = np.array([[ 0, 0, 0, 0, 0, 1, 0, 0, -1], # s -w_w
[1.5, 0, 0, 0, 0, .5, 0, .5, -1]]) # 1.5a + .5w + .5s -w_w
q = np.array([np.log(45.42), np.log(8.71e-5*Nlift*np.sqrt(W0)/tau)])
constr = [A@x == b, C@x <= d,
cp.log_sum_exp(E@x+f) <= 0,
cp.log_sum_exp(G@x+h) <= 0,
cp.log_sum_exp(P@x+q) <= 0]
cp.Problem(obj, constr).solve(verbose=True)
AR = np.exp(x.value[0])
C_D = np.exp(x.value[1])
C_f = np.exp(x.value[2])
C_L = np.exp(x.value[3])
Re = np.exp(x.value[4])
S = np.exp(x.value[5])
V = np.exp(x.value[6])
W = np.exp(x.value[7])
W_w = np.exp(x.value[8])
print(f'AR = {AR}')
print(f'S = {S}')
print(f'C_D = {C_D}')
print(f'C_L = {C_L}')
print(f'C_f = {C_f}')
print(f'Re = {Re}')
print(f'W = {W}')
print(f'W_w = {W_w}')
print(f'V = {V}')
===============================================================================
CVXPY
v1.3.1
===============================================================================
(CVXPY) May 04 09:46:15 AM: Your problem has 9 variables, 5 constraints, and 0 parameters.
(CVXPY) May 04 09:46:15 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) May 04 09:46:15 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) May 04 09:46:15 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) May 04 09:46:15 AM: Compiling problem (target solver=ECOS).
(CVXPY) May 04 09:46:15 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS
(CVXPY) May 04 09:46:15 AM: Applying reduction Dcp2Cone
(CVXPY) May 04 09:46:16 AM: Applying reduction CvxAttr2Constr
(CVXPY) May 04 09:46:16 AM: Applying reduction ConeMatrixStuffing
(CVXPY) May 04 09:46:16 AM: Applying reduction ECOS
(CVXPY) May 04 09:46:16 AM: Finished problem compilation (took 2.119e-01 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) May 04 09:46:16 AM: Invoking solver ECOS to obtain a solution.
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) May 04 09:46:16 AM: Problem status: optimal
(CVXPY) May 04 09:46:16 AM: Optimal value: 6.027e+00
(CVXPY) May 04 09:46:16 AM: Compilation took 2.119e-01 seconds
(CVXPY) May 04 09:46:16 AM: Solver (including time spent in interface) took 1.050e-03 seconds
AR = 12.6972479450141
S = 12.075070874016546
C_D = 0.0230980869270328
C_L = 0.651221745014052
C_f = 0.0038574654699476585
Re = 2598055.1883249437
W = 7188.531040088981
W_w = 2248.5311457240487
V = 38.55433592085682
Contents partially taken from Hoburg and Abbeel, "Geometric Programming for Aircraft Design Optimization," AIAA Journal, 2014.