Basic Wing Aerostructural Optimization Using an Optimization Framework

In this problem we will perform a (very basic, noncoupled) aerostructural optimization of a large UAV wing as pictured below.

In particular, our goal is to be introduced to using an optimization framework, in this case Isight. Within the framework we will connect a CAD package (CATIA), a finite element solver (ANSYS), an analysis script in Matlab (for the aerodynamics), and an external optimizer in Matlab. A student in ME 578 (CAD/Cam Applications) will provide the experience with the Isight/CAD/ANSYS setup, and you will provide the experience with the problem statement and optimization setup.

There are six design variables in this problem: the chord $c$ at the root and the tip of the wing, the airfoil thickness $t$ at the root and tip, and the skin thickness $t_s$ at the root and tip. It is assumed that the chord and thicknesses vary linearly between the root and tip.

$$x = [c_r, c_t, t_r, t_t, {t_s}_r, {t_s}_t]$$

A Matlab script is provided (at the bottom of this page) that, given the design variables, computes the lift-to-drag ratio of the wing $L/D$, the lift coefficient distribution at a maneuver condition $c_l$, the lift distribution at a maneuver condition (simplified to a pressure $p$ that varies spanwise at defined locations $yp$), the take-off weight of the aircraft $W_t$, and the zero-fuel weight of the aircraft excluding the structural weight of the wing $W_0$.

You will also need to pass the design variables to CAD to setup the geometry. While aerodynamics depends on the full blade shape, for the structural analysis you only need to model the hollow tapered box beam (gray part of figure above). The outer skin around the box beam carries little to no load. The main spar is approximated by the tapered box beam shown below. You can compute the weight in ANSYS, but it's probably easier to do it in CAD where you have the actual geometry. Remember to multiply by the acceleration of gravity, and a factor of 2 (because you are only modeling half of the wing in CAD).

The geometry from CAD and the applied load from Matlab must then be passed to the finite element solver. ANSYS will compute the wing structural weight $W_w$, and the stress $\sigma$ at a few locations along the span (or at a minimum just at the root). An overview of the connection graph is shown below.

The objective of the optimization is to maximize the range (distance the UAV can fly). Assuming a constant engine efficiency, the range is proportional to:

$$ J = \frac{L}{D} \ln{\frac{W_t}{W_w + W_0}}$$

The constraints are

\begin{align*} c_l & \le {c_l}_{max} \quad\text{(stall)}\\ \eta \sigma & \le \sigma_y \quad\text{(bending stress)}\\ t_t/c_t & \ge 6\% \quad\text{(manufacturing)}\\ {t_s}_r/t_r & \le 25\% \quad\text{(volume)}\\ {t_s}_t/t_t & \le 25\% \quad\text{(volume)}\\ \end{align*}

You will likely use a relatively coarse discretization in ANSYS, and as a result the function output will be a little noisy. The optimization will should still work reasonbly well but convergence will be slow. You'll welcome to run over night or manually stop it when you feel the design is "good enough".

Parameters that you will need in this analysis are listed below.

Description Variable Value
maximum lift coefficient ${c_l}_{max}$ 1.3
safety factor $\eta$ 1.5
yield stress $\sigma_y$ $1\times10^6$ psf
spar cap material density $\rho$ 25 slugs/ft$^3$
chord fraction of structural box $c_{box}$ 0.2
wing span b 20 ft
thickness of shear webs $t_w$ 0.01 ft

The other piece you will need is a reasonable starting point. Generally, we use our engineering knowledge to pick something rational. For an aircraft wing, important nondimensional parameters are its aspect ratio is defined as: $ AR = b^2/S$, and its tip-speed ratio as $\lambda = c_t/c_r$. Typical values for an efficient aircraft would be around $AR=9.0$ and $\lambda=2.0$. From those definitions, the formula for the area for this simple tapered wing ($S = 0.5(c_r+c_t)b$) and knowing the span we can compute a reasonable root and tip chord:

$$ c_r = \frac{2 b}{AR (1 + \lambda)} $$$$ c_t = \lambda c_r $$

Next, a typical airfoil section is somewhere around 13% thick at the root and 9% at the tip:

$$ t_r = 0.13 c_r $$$$ t_t = 0.09 c_t $$

Finally, the skin thickness is probably around 10% of the section thickness:

$$ {t_s}_r = 0.1 t_r $$$$ {t_s}_t = 0.1 t_t $$

These are not optimal values by any means, just reasonable starting points.

The simplified Matlab analysis code you will need is listed below (you don't need to worry about the details internal to the function, just need to pass the inputs and use the outputs):

In [ ]:
function [L_D, cl, yp, p, Wt, W0] = aero(cr, ct, tr, tt)
cr, ct: airfoil chord at root and tip
tr, tt: airfoil thickness at root and tip

L_D: lift to drag ratio (L/D)
cl: lift coefficient distribution at multiple points along wing (yp)
yp: locations along span where pressure is defined
p: spanwise pressure distribution at each yp location
Wt: total weight of aircraft
W0: zero fuel weight of wing (minus wing structural weight)

% --- parameters ---
b = 20.0;
W_b2 = 12.0;  % a typical value of W/b^2 for transport aircraft
CL = 0.5;  % cruise lift coefficient
nmvr = 2.0; % manuever load factor
N = 200;  % number of points in discretization
cf = 0.008;  % skin friction coefficient
K = 0.38;  % factor of visous portion of span efficiency
e_reduce = 0.96;  % span efficiency reduction to things like prop wash
w0 = 1.0;  % weight per unit area of nonstructural portions of wing
% ------------------

% setup discretization
y = linspace(0, b/2, N);
c = linspace(cr, ct, N);
t = linspace(tr, tt, N);
ycoarse = linspace(0, b/2, N/10);

% lift/weight
Wt = W_b2*b^2;
L = Wt;

% planform area
S = 0.5*(cr + ct)*b;

% dynamic pressure
q = L/(CL*S);

% assume cruise elliptic lift distribution
l0 = 4*L/(b*pi);
l = l0*sqrt(1 - (2*y/b).^2);

% % equivalent manuever load force at centroid
% F = nmvr*L;
% yF = 4/(3*pi)*b/2;

% lift coefficient distribution
cl = l./(q*c);

% maneuver lift distribution
clmvr = cl + (nmvr - 1)*CL;  % approximation. assumes elliptic planform
lmvr = clmvr*q.*c;

% constrain maneuver cl but let's make it coarser
cl = interp1(y, clmvr, ycoarse);

% turn into a pressure load for ANSYS
p = interp1(y, lmvr./c, ycoarse);
yp = ycoarse;

% viscous drag
tc = t./c;
k = 6.8*tc.^2 + 1.35*tc + 1;  % form factor
SwetS = 2*(1+0.2*tc);  % area ratio
CDp = 2*trapz(y, k*cf.*SwetS.*c)/S;  % chord weighted

% induced drag
AR = b^2/S;
e_inv = 1.0;  % b.c. assumed elliptic
e = 1/(1/e_inv + pi*AR*K*CDp);
e = e_reduce*e;
CDi = CL^2/(pi*AR*e);

% lift to drag ratio
L_D = CL/(CDp + CDi);

% non wing weight
W0 = 0.5*Wt + w0*S;  % assuming half of take off weight for fuel, and adding nonstructural weight proportional to area