Low-fidelity method for rapid aerostructural optimisation and design-space exploration of planar wings

Abstract

During early phases of wing design, analytic and low-fidelity methods are often used to identify promising design concepts. In many cases, solutions obtained using these methods provide intuition about the design space that is not easily obtained using higher-fidelity methods. This is especially true for aerostructural design. However, many analytic and low-fidelity aerostructural solutions are limited in application to wings with specific planforms and weight distributions. Here, a numerical method for minimising induced drag with structural constraints is presented that uses approximations that apply to unswept planar wings with arbitrary planforms and weight distributions. The method is applied to the National Aeronautics and Space Administration (NASA) Ikhana airframe to show how it can be used for rapid aerostructural optimisation and design-space exploration. The design space around the optimum solution is visualised, and the sensitivity of the optimum solution to changes in weight distribution, structural properties, wing loading and taper ratio is shown. The optimum lift distribution and wing-structure weight for the Ikhana airframe are shown to be in good agreement with analytic solutions. Whereas most modern high-fidelity solvers obtain solutions in a matter of hours, all of the solutions shown here can be obtained in a matter of seconds.

NOMENCLATURE

A

beam cross-sectional area

$A_{n}$

Fourier coefficients in the lifting-line solution for the section-lift distribution, Equation (1)

b

wingspan

$B_{n}$

Fourier coefficients in the lifting-line solution for the dimensionless section-lift distribution, Equation (1)

c

local wing section chord length

$C_{\delta}$

shape coefficient for the deflection-limited design, Equation (15)

$C_{\sigma}$

shape coefficient for the stress-limited design, Equation (5)

$D_{i}$

wing induced drag

$E$

modulus of elasticity of the beam material

h

height of the beam cross-section

I

beam section moment of inertia

K

scaling coefficient in the equation for the fuel distribution, Equation (21)

L

total wing lift

$\widetilde{L}$

local wing section lift

$\widetilde{M}_{b}$

local wing section bending moment

$n_{a}$

load factor, g

$n_{g}$

limiting load factor at the hard-landing design limit

$n_{m}$

limiting load factor at the manoeuvring-flight design limit

$R_{A}$

wing aspect ratio

$R_{T}$

wing taper ratio

S

wing planform area

$S_{b}$

proportionality coefficient between $\widetilde{W}_{s}(z)$ and $\widetilde{M}_{b}(z)$ having units of length squared, Equations (5) and (15)

$t_{\max}$

maximum thickness of the local aerofoil section

$V_{\infty}$

freestream airspeed

w

width of the beam cross-section

$w_{\max}$

maximum allowable width of the beam cross-section

W

aircraft gross weight

$W_{f}$

gross weight of fuel

$W_{n}$

aircraft net weight, defined as $W-W_{s}$

$W_{r}$

that portion of W_n carried at the wing root

$W_{s}$

total weight of the wing structure required to support the wing bending-moment distribution

$\widetilde{W}_{n}$

net weight of the wing per unit span, i.e. total wing weight per unit span less $\widetilde{W}_{s}$

$\widetilde{W}_{s}$

weight of the wing structure per unit span required to support the wing bending-moment distribution

z

spanwise co-ordinate relative to the midspan

$\gamma$

specific weight of the beam material

$\delta$

local wing deflection

$\delta_{\max}$

maximum wing deflection

$\theta$

change of variables for the spanwise co-ordinate, Equation (1)

$\rho$

air density

$\sigma_{\max}$

maximum longitudinal stress

1.0 INTRODUCTION

When designing a wing for minimum drag, low-fidelity tools are useful for rapid design-space exploration and for gaining important insight into how the design variables, parameters and constraints influence the optimum solution. Designers often rely on rules-of-thumb based on these insights during the conceptual and preliminary design phases. In many cases, low-fidelity solutions have been shown to be in good agreement with experimental data and computational fluid dynamics^(1–8), while providing significantly more mathematical and physical insight than higher-fidelity models. For example, the well-known elliptic lift distribution, which minimises induced drag on an unswept planar wing with fixed weight and wingspan, was first identified from analytic solutions based on lifting-line theory^(9,10) by Prandtl⁽⁹⁾ and later by Munk⁽¹¹⁾. The elliptic lift distribution remains a common benchmark in many mid- and high-fidelity computational studies^(12–19). However, the elliptic lift distribution does not minimise drag under all conditions^(20–29). In particular, when structural effects are considered, drag is typically minimised using a non-elliptic lift distribution that depends on the design constraints^{(12,13,16–18,30–49)} . Low-fidelity and analytic aerostructural methods are valuable for identifying these non-elliptic lift distributions and for understanding how structural considerations affect the minimum-drag solution.

There are many mid- and high-fidelity computational studies for minimising drag under structural constraints that include solutions with non-elliptic lift distributions^{(12,13,16–18,30–38)} . However, there are relatively few studies that approach this multidisciplinary problem from an analytic or low-fidelity point of view^(39–49). Prandtl seems to be the first do so, minimising induced drag with fixed lift and moment of inertia of gross lift⁽³⁹⁾. Jones later⁽⁴⁰⁾ sought to minimise induced drag under the constraints of fixed gross lift and root bending moment in cruise. Pate and German⁽⁴¹⁾ constrained the root bending moment at a given off-design lift coefficient. DeYoung⁽⁴²⁾ replaced Jones’ root-bending-moment constraint with a constraint on the bending moment at a prescribed spanwise location. Jones and Lasinski⁽⁴³⁾ constrained the integrated bending moment. Klein and Viswanathan^(44,45) considered both root and integrated bending moment⁽⁴⁴⁾ and included the effects of shear on the wing-structure weight⁽⁴⁵⁾. Löbert⁽⁴⁶⁾ introduced a constraint based on the ratio of the bending-moment distribution and the wing-section thickness. More recently, Phillips et al.^(47,48) and Taylor and Hunsaker⁽⁴⁹⁾ minimised induced drag under constraints of fixed gross weight^(47,49) , fixed net weight^(48,49) , fixed wing loading^(47–49) and fixed stall speed⁽⁴⁸⁾, including the effects of the planform shape on the wing-structure weight and the effects of the wing weight distribution on the bending moments.

Each of the studies in Refs. (39–49) includes assumptions that may not be representative of all aircraft. For example, Refs. (39,44,45) include assumptions about the proportionality between the wing-structure weight and the bending moments that correspond to rectangular wings. References (39–46) include the assumption that the bending moments are caused by the lift alone, which limits their application to wings with negligible structural or payload weight. The formulations given by Phillips et al.^(47,48) and Taylor and Hunsaker⁽⁴⁹⁾ are arguably more general than those given in Refs. (39–46). Still, to obtain analytic solutions, Phillips et al.^(47,48) and Taylor and Hunsaker⁽⁴⁹⁾ limited their results to specific wing planforms with a single ideal weight distribution.

The purpose of this paper is to present a low-fidelity numerical method that extends the work of Phillips et al.^(47,48) and Taylor and Hunsaker⁽⁴⁹⁾ to more practical aircraft configurations with arbitrary planforms and weight distributions. We will apply the method to a high-endurance unmanned aircraft configuration to demonstrate how it can be used for rapid conceptual design and for gaining intuition about the aerostructural design space. The present work builds on the approach taken by Prandtl⁽³⁹⁾ and Phillips et al.^(47,48) . Therefore, we will first briefly review the work of these authors.

2.0 ANALYTICAL FOUNDATION

Using Prandtl’s classical lifting-line theory^(9,10) , the dimensionless spanwise section-lift distribution on a finite wing with no dihedral or sweep immersed in a uniform flow can be written as⁽⁴⁷⁾

(1) \begin{equation}\frac{b\widetilde{L}(\theta)}{L}=\frac{4}{\pi}\left[sin{\theta}+\sum_{n=2}^{\infty}B_{n}\sin(n\theta)\right]\!;\quad B_{n}\equiv\frac{A_{n}}{A_{1}},\quad \theta\equiv\cos^{-1}(-2z/b)\end{equation}

where $B_{n}$ are normalised Fourier coefficients. Below stall, any lift distribution can be produced by a twisted wing of any planform if the correct twist distribution is used⁽⁵⁰⁾. Therefore, in this paper, the lift distribution and the planform are treated as independent parameters, related through the wing twist, which is assumed to be correctly designed to achieve the desired lift distribution. In steady-level flight, the drag induced by such a wing can be written as

(2) \begin{equation}D_{i}=\frac{2(W/b)^{2}}{\pi\rho V^{2}_{\infty}}\left(1+\sum_{n=2}^{\infty}nB_{n}^{2}\right)\end{equation}

where W is the wing weight, and b is the wingspan. Because this study focuses on minimising induced drag, we will neglect the effects of viscous drag.

Equation (2) reveals that induced drag depends on the weight, wingspan and lift distribution. For a fixed ratio of weight to wingspan, Equation (2) is minimised with a lift distribution having $B_{n}=0$ for all $n>1$ , which gives the well-known elliptic lift distribution. If weight and wingspan are allowed to vary, the induced drag can be reduced by increasing wingspan or decreasing wing weight. However, as wingspan increases, the weight of the wing structure required to support the bending moments also increases, which increases the total weight. Certain lift distributions that shift lift inboard can alleviate bending moments near the wingtips, allowing a higher wingspan with no increase in wing-structure weight. Therefore, to fully minimise Equation (2) for a given flight condition, the weight, wingspan and lift distribution must all be considered.

In 1933, Prandtl⁽³⁹⁾ identified a bell-shaped lift distribution having $B_{2}=0$ , $B_{3}=-1/3$ and $B_{n}=0$ for $n>3$ that minimises induced drag for rectangular wings under constraints of fixed gross weight and moment of inertia of gross weight. Prandtl assumed that the wing-structure weight distribution $\widetilde{W}_{s}(z)$ is related to the bending-moment distribution $\widetilde{M}_{b}(z)$ by a spanwise-invariant proportionality coefficient $S_{b}$ , i.e.

(3) \begin{equation}\widetilde{W}_{s}(z)=\frac{\widetilde{M}_{b}(z)}{S_{b}}\end{equation}

This assumption is best matched by a rectangular wing with a constant thickness-to-chord ratio⁽³⁹⁾. Prandtl also assumed that the bending-moment distribution is a function of the lift distribution alone. Under the constraints of these assumptions, Prandtl’s 1933 lift distribution allows an increase in wingspan of 22.5% and a reduction in induced drag of 11.1% when compared with that of the elliptic lift distribution with the same wing-structure weight. However, Prandtl acknowledged that his formulation of the problem may not be the most appropriate for practical wing designs⁽³⁹⁾.

Phillips et al.^(47,48) reformulated the problem with more practical assumptions and constraints. They pointed out that at each spanwise location, the wing bending moments are a function of the lift distribution, the net-weight distribution $\widetilde{W}_{n}(z)$ of all non-structural components carried in the wing, and the wing-structure weight distribution $\widetilde{W}_{s}(z)$ according to the relation⁽⁴⁷⁾

(4) \begin{equation}\widetilde{M}_{b}(z)=\int\limits^{b/2}_{z^{\prime}=z}[\widetilde{L}(z^{\prime})-n_{a}\widetilde{W}_{n}(z^{\prime})-n_{a}\widetilde{W}_{s}(z^{\prime})](z^{\prime}-z)dz^{\prime},\quad \text{for}\ z\geq0\end{equation}

where n _a is the load factor. The wing structure must be designed to support the bending moments during a high-load manoeuvre with a positive load limit n _m and during a hard landing with a negative load limit n _g . Assuming that all of the wing bending moments are supported by a single, vertically symmetric beam in pure bending with maximum allowable stress $\sigma_{\max}$ , the weight of the wing structure required to support the bending moments can be written⁽⁴⁷⁾

(5) \begin{equation}W_{s}=2\int\limits_{0}^{b/2}\frac{\left|\widetilde{M}_{b}(z)\right|}{S_{b}(z)}dz;\quad S_{b}(z)=\frac{C_{\sigma}[t_{\max}(z)/c(z)]c(z)\sigma_{\max}}{\gamma},\quad C_{\sigma}=\frac{2I(h/t_{\max})}{Ah^{2}}\end{equation}

where c(z) is the section chord-length distribution, $\gamma$ is the specific weight of the beam material, $t_{\max}/c$ is the maximum-thickness-to-chord ratio of the local aerofoil section, and $C_{\sigma}$ is a beam shape factor. A list of shape factors for common beam cross sections is given in Ref. (47). For deflection-limited designs, Equation (5) can be rewritten as⁽⁴⁷⁾

(6) \begin{equation}W_{s}=2\int\limits_{0}^{b/2}\frac{\left|\widetilde{M}_{b}(z)\right|}{S_{b}(z)}dz;\quad S_{b}(z)=\frac{C_{\sigma}E[t_{\max}(z)/c(z)]^{2}\delta_{\max}}{\gamma(W/S)^{2}}\frac{W^2}{b^{4}},\quad C_{\sigma}=\frac{8I(h/t_{\max})^{2}}{Ah^{2}}\end{equation}

where $C_{\delta}$ is the beam shape factor for the deflection-limited design and $\delta$ _max is the maximum allowable vertical wingtip deflection. Although vertical deflection limits are seldom explicitly enforced in practice, excessive vertical wingtip deflection can result in serious adverse effects, including wingtip strike at landing and dynamic instabilities during flight. Therefore, we will include both stress and vertical deflection limits in this paper. Nevertheless, the deflection limits in this paper are for structural sizing only. The static aeroelastic effects of structural bending and torsion are not explicitly considered. Instead, we assume that these effects can be corrected using wing twist.

The total weight of the wing is the sum of the wing-structure weight and the net weight of all non-structural components, i.e.

(7) \begin{equation}W=W_{s}+W_{n}\end{equation}

The net weight W _n is found from the relation

(8) \begin{equation}W_{n}=W_{r}+\int\limits^{b/2}_{z=-b/2}\widetilde{W}_{n}(z)dz\end{equation}

where W _r is the portion of the net weight carried at the wing root. The bending moments are minimised when the net weight is distributed according to the weight constraints given by⁽⁴⁷⁾

(9) \begin{equation}\widetilde{W}_{n}(z)=(W-W_{r})\frac{\widetilde{L}(z)}{L}-\widetilde{W}_{s}(z)\end{equation}

(10) \begin{equation}W_{r}=\frac{n_{g}-1}{n_{m}+n_{g}}W\end{equation}

For a rectangular wing having the weight distribution from Equation (9), Equations (5) and (6) can be evaluated analytically. Assuming that the wing loading is fixed and a single lift distribution is used at all flight phases, Phillips et al.⁽⁴⁸⁾ showed that induced drag is minimised with a lift distribution having $B_{2}=0$ , $B_{3}=-3/8+\sqrt{9/64-1/12}$ , with $B_{n}=0$ for $n>3$ for the stress-limited design and $B_{2}=0$ , $B_{3}=-3/7+\sqrt{9/49-1/21}$ , with $B_{n}=0$ for $n>3$ for the deflection-limited design.

In this paper, we extend the work of Phillips et al.^(47,48) and present a method for minimising induced drag for wings with non-rectangular planforms and weight distributions other than Equation (9). It should be remembered that the present method maintains the assumptions associated with lifting-line theory, including a planar wing with zero sweep and moderate to high aspect ratio. For other wing configurations, modifications to this method may be needed.

3.0 WING-STRUCTURE WEIGHT AND INDUCED DRAG

For the stress-limited design of a wing with a non-rectangular planform and a weight distribution other than Equation (9), the integrals in Equations (4) and (5) must often be evaluated numerically. Moreover, for any given flight condition, Equations (4) and (5) show that the wing bending moments and wing-structure weight distribution are coupled. Therefore, for a wing with any weight distribution other than Equation (9), a numerical iterative method is required to compute the wing-structure weight. The induced drag can then be found by using Equation (7) in Equation (2). An implementation of one such iterative process is given by Taylor et al.⁽⁵¹⁾ for the stress-limited design.

For deflection-limited designs, the vertical spar deflection can be found using the relation⁽⁴⁷⁾

(11) \begin{equation}\frac{d^{2}\delta}{dz^{2}}=\frac{2\sigma_{\max}}{Eh(z)}\end{equation}

where E is the modulus of elasticity of the beam material. For any spanwise-symmetric load distribution, the boundary conditions on Equation (11) are

(12) \begin{equation}\delta(0)=0, \frac{d\delta}{dz}\bigg|_{z=0}=0\end{equation}

Integrating Equation (11) subject to Equation (12), the deflection at any spanwise location $z_{0}$ becomes

(13) \begin{equation}\delta(z_{0})=\frac{2\sigma_{\max}}{E}\int\limits_{0}^{z_{0}}\int\limits_{0}^{z}\frac{1}{h(z')}dz'dz\end{equation}

If both manoeuvring and hard-landing design limits are considered, maximum deflection always occurs at the wingtips. Using Equation (13), the deflection at the wingtip is

(14) \begin{equation}\delta_{\max}=\frac{2\sigma_{\max}}{E}\int\limits_{0}^{b/2}\int\limits_{0}^{z}\frac{1}{h(z')}dz'dz\end{equation}

Because aerofoil thickness is typically a fraction of the chord length, the beam-height distribution h(z) is typically related to the chord distribution. If the beam-height or chord distribution is an arbitrary function of spanwise location. Equation (14) must be evaluated using numerical methods.

Using Equation (14) to replace $\sigma_{\max}$ in Equation (5), the wing-structure weight required to support the bending moments for the deflection-limited design can be written

(15) \begin{equation}W_{s}=2\int\limits^{b/2}_{0}\frac{\left|\widetilde{M}_{b}(z)\right|}{S_{b}(z)}dz;\quad S_{b}(z)=\frac{C_\delta E[t_{\max}(z)/c(z)]c(z)\delta_{\max}}{8\gamma\int^{b/2}_{0}\int^{z}_{0}[t_{\max}(z')/c(z')]^{-1}c(z')^{-1}dz'dz},\quad C_{\delta}=\frac{8I(h/t_{\max})^{2}}{Ah^{2}}\end{equation}

Like Equation (5), Equation (15) is coupled with the bending-moment distribution. Thus, an iterative solver is needed to compute the wing-structure weight for the deflection-limited design.

If Equation (5) predicts a wing-structure weight that is greater than that predicted by Equation (15), the design is stress limited; if Equation (15) gives a value greater than Equation (5), the design is deflection limited. Because the limiting constraint depends on the design parameters, both stress and deflection limits must be considered at each spanwise location. However, recall that in this study, the aerodynamic effects of structural bending and twist are not included.

4.0 NUMERICAL METHODOLOGY

Here, we present a method to iteratively compute the wing-structure weight and minimise induced drag. This method is similar to that given by Taylor et al.⁽⁵¹⁾, but here we will include the deflection-limited design and several additional constraints that were not considered in Ref. (51).

4.1 Solving for Wing-Structure Weight

A fixed-point iteration scheme is used to compute the wing-structure weight and bending-moment distribution. An initial guess for the wing-structure weight is used in Equation (4) to calculate the section bending-moment distribution for both the manoeuvring and hard-landing limits. At each section, the limit that produces a higher-magnitude section bending moment is the design limit. The limiting section bending moment is used in Equations (5) and (15) to predict the section wing-structure weight for the stress- and deflection-limited designs. At each section, the limiting wing-structure weight is then passed back as the guess for the next iteration. The process is repeated until the wing-structure weight converges within some specified tolerance. For the purposes of this study, an initial guess of $\widetilde{W}_{s}(z)=0$ provides good results. The process is summarised as follows:

1. 1. Input b, $\widetilde{L}(z)/L$ , W _r , $\widetilde{W}_{n}(z)$ , c(z), $t_{\max}(z)/c(z)$ , $\gamma$ , E, $\sigma_{\max}$ , $\delta_{\max}$ , n _m , n _g , $C_{\sigma}$ and $C_{\delta}$ .

2. 2. Calculate the total weight using Equation (7). For the initial guess, use $\widetilde{W}_{s}(z)=0$ , W _s = 0.

3. 3. Calculate the total net weight using Equation (8).

4. 4. Calculate the manoeuvring and hard-landing bending-moment distributions using Equation (4).

5. 5. Using the higher-magnitude section bending moment from step 4 in Equations (5) and (15), calculate the wing-structure weight distribution for the stress- and deflection-limited designs.

6. 6. Calculate the total wing-structure weight by integrating either Equation (5) or (15).

7. 7. Repeat steps 2 through 6 until the wing-structure weight has converged to within a specified tolerance.

Once the wing-structure weight is known, the induced drag is calculated using Equation (2). A schematic of the process is shown in Fig. 1. Note that after the first iteration, step 3 is only required if the net weight is a function of the wing-structure weight, as it is in Equation (9). In this paper, this special case will be used only for benchmarking the wing-structure weight solver against analytic solutions.

Figure 1. Schematic of the iterative wing-structure weight solver.

In general, any high-order integration scheme can be used to evaluate the integrals in Equations (4), (5), (8) and (15). In this study, the composite Simpson’s rule is used. The wing is discretised using the cosine clustering scheme given in Equation (1), with even spacing in $\theta$ . The resulting grid is shown in Fig. 2. Using Simpson’s rule, the wing-structure weight is evaluated as⁽⁵¹⁾

(16) \begin{align}W_{s}&=-b(\theta_{m}-\theta_{0})\nonumber\\[3pt]&\quad\times \left(\frac{\widetilde{W}_{s,0}\sin\theta_{0}+4\sum^{\rm m-1}_{\rm i=odd}\widetilde{W}_{s,i}\sin{\theta_i}+2\sum^{\rm m-2}_{\rm i=even}\widetilde{W}_{s,i}\sin\theta_{i}+\widetilde{W}_{s,m}\sin\theta_{m}}{3m}\right)\end{align}

where m is the number of nodes, and $\widetilde{W}_{s,i}$ is evaluated from Equation (5) or Equation (15), i.e.

(17) \begin{equation}\widetilde{W}_{s,i}=\frac{|\widetilde{M}_{b,i}|}{S_{b,i}};\quad S_{b,i}=\frac{C_{\sigma}[t_{\textrm{max},i}/c_{i}]c_{i}\sigma_{\max}}{\gamma}\end{equation}

(18) \begin{equation}\widetilde{W}_{s,i}=\frac{|\widetilde{M}_{b,i}|}{S_{b,i}}; S_{b,i}=\frac{C_{\delta}E[t_{\textrm{max},i}/c_{i}]c_{i}\sigma_{\max}}{4b\gamma\int_{0}^{\pi}\int_{0}^{z}[t_{\max}(\theta')/c(\theta')]^{-1}\sin\theta'\sin\theta d\theta' d\theta}\end{equation}

Figure 2. Discretisation of a tapered semispan with 40 nodes and cosine clustering near the wing tip.

The integral in the denominator of Equation (18) is also evaluated using Simpson’s rule, and $\widetilde{M}_{b,i}$ is found from Equation (4), i.e.⁽⁵¹⁾

(19) \begin{align}\widetilde{M}_{b,i}&=-\frac{b}{2}(\theta_{m}-\theta_{i})\left(\frac{f_{i}+4\sum\limits^{m-1}_{j=i+\textrm{odd}}+2\sum\limits^{m-1}_{j=i+\textrm{even}}f_{j}+f_{m}}{3(m-i)}\right);\nonumber\\[3pt] f_{j}&=\frac{b}{2}\sin\theta_{j}(\widetilde{L}_{j}-n_{a}\widetilde{W}_{n,j}-n_{a}\widetilde{W}_{s,j})(\cos\theta_{j}-\cos\theta_{i})\end{align}

Note that Simpson’s rule requires even grid spacing. Therefore, Equations (16)–(19) are written in terms of $\theta$ .

Figure 3 shows the results of a grid-resolution study for the iterative wing-structure weight solver using a wing with the parameters $R_{T}=0.5$ , $b=66.0$ ft, $S=267.3\ \mathrm{ft}^{2}$ , $t/c=0.1875$ , $C_{\sigma}=0.165$ , $C_{\delta}=0.653$ , $\sigma_{\max}=25\times10^{3}$ psi, $\delta_{\max}=3.5$ ft, $E=10.0\times10^{6}$ psi, $\gamma=0.10\ \mathrm{lbf/in}^{3}$ , $W_{r}=4500$ lbf, $W_{n}=7500$ lbf, $n_{m}=n_{g}=3.75$ and the weight distribution given by Equation (9). Results were compared using grids with node counts ranging between 10 and 1280, and Richardson extrapolation⁽⁵²⁾ was used to project a fully grid-resolved value from the results obtained with 160, 320 and 640 nodes. Above 40 nodes, the method shows second-order convergence, meaning that as the grid size is halved, the solution error is approximately reduced to one-fourth the previous value. The extrapolated value differs from the analytic solution⁽⁴⁹⁾ by only 0.001%. With as few as 160 nodes, the predicted wing-structure weight falls within 0.003% of the extrapolated value. Therefore, 160 nodes will be used for all subsequent results. With 160 nodes, the total predicted wing-structure weight matches the analytic solution to within 0.004%, and Fig. 4 shows that the predicted wing-structure weight distribution is in good agreement with the analytic solution.

Figure 3. Grid-resolution results for the iterative wing-structure weight solver.

Figure 4. Comparison of the wing-structure weight predicted by the numerical wing-structure weight solver and the analytic solution from Ref. (49).

4.2 Minimising Induced Drag in an Optimisation Framework

The induced drag from the wing-structure solver can be used as an objective function in an optimisation framework similar to that shown in Fig. 5. Any of the parameters from Equations (2), (5), (7) or (15) could be used as design variables. However, in this study, we will use only the lift distribution (B _n ) and wingspan (b). Note that in the previous sections, the lift distribution is assumed to be spanwise symmetric (B _n = 0 for all even n). Therefore, for the remainder of this study, we will assume that the even Fourier coefficients are identically zero.

Figure 5. Example optimisation framework for minimising induced drag using wingspan and lift distribution.

The optimisation process is summarised as follows: an initial guess is made for the design variables b and B _n ; the wing-structure weight and induced drag are computed using the methodology explained in the previous section; the design variables are updated using an optimisation method of the user’s choosing, subject to relevant constraints; and the updated design variables are fed back to the wing-structure weight solver. The process is repeated until the induced drag converges within some specified tolerance. Because the relationship between induced drag and the design variables is well behaved, any gradient-based method with appropriate constraints should be adequate for updating the design variables b and B _n . The method used in this study for updating b and B _n is discussed in the following section.

The choice of design constraints can have a significant impact on the minimum-induced-drag solution^(47,48) . In this study, we will consider only a few example constraints proposed by Phillips et al.^(47,48) , including an all-positive spanwise lift distribution, fixed net weight and fixed wing loading. We will also constrain spar height h and width w, as explained in Ref. (53), to ensure that the spar fits within the local aerofoil section. The optimisation problem can be summarised as follows:

(20) \begin{equation}\left\{\begin{array}{l@{\quad}l}\mathrm{minimise:} & D_{i}\\\\[-9pt]\mathrm{with\ respect\ to:} & b,B_{n}(n\ \textrm{odd})\\\\[-9pt]\mathrm{subject\ to:} & L-W=0\\\\[-9pt]& \sigma_{\max}-\sigma(z)\geq0\\\\[-9pt]&\delta_{\max}-\delta(b/2)\geq0\\\\[-9pt]&\widetilde{L}(z)\geq0\\\\[-9pt]& 1-h(z)/t_{\max}(z)\geq0\\\\[-9pt]& w_{\max}/c-w(z)/c(z)\geq0\\\\[-9pt]& W_{n_{0}}-W_{n}=0\\\\[-9pt]& \widetilde{W}_{n_{0}}(z)-\widetilde{W}_{n}(z)=0\\\\[-9pt]&(W/S)_{0}-W/S=0\end{array}\right.\end{equation}

where the subscript 0 indicates that the parameter value is prescribed. The first three constraints in Equation (20) are enforced implicitly in Equations (2), (5) and (15). The remaining constraints can be enforced as explained in Ref. (53).

5.0 RESULTS

As an example of minimising induced drag for a wing with a non-rectangular planform and a net-weight distribution other than Equation (9), consider the NASA Ikhana airframe^(54–57). Ikhana has a linearly tapered wing with a wingspan b = 66 ft, an aspect ratio R _A = 16.296 and a taper ratio R _T = 0.421. A generic instrumentation pod weighing 500 lbf⁽⁵⁷⁾ is sometimes mounted at a hard point outboard of the wing root. Assuming that all of the fuel is distributed in fuel bladders that extend to 83.1% semispan⁽⁵³⁾, the net-weight distribution can be approximated as

(21) \begin{equation}\widetilde{W}_{n}(z)=Kc(z)^{2}\end{equation}

where K is a scaling constant that depends on the length of the fuel bladder and the weight of the fuel carried in the wing. Using Equation (21) in Equation (8) gives a relationship that can be solved to find K for a fuel bladder that extends to 83.1% semispan with a given fuel weight $W_{f}$ , i.e.

(22) \begin{align}W_{f}=2K\int\limits^{0.831b/2}_{0}c(z)^{2}dz\\[-17pt]\nonumber\end{align}

For this study, we will consider two example Ikhana configurations in steady level flight at sea level with a cruise velocity of 287 ft/s⁽⁵⁵⁾. The first configuration has 3000 lbf of fuel distributed according to Equation (21) in fuel bladders spanning 83.1% semispan with no instrumentation pod. This gives a scaling constant $K=2.8212$ . The second example configuration includes a generic instrument pod mounted on each wing at hard points located at 25% semispan that each cover 1ft spanwise. To maintain the same fixed net weight as the no-pod configuration, the fuel weight is reduced to 2000 lbf, which gives a scaling constant $K=1.8808$ . The resulting net-weight distribution is shown in Fig. 6. All other parameters for both configurations are given in Table 1. Note that the values for $C_{\sigma}$ and $C_{\delta}$ correspond to a beam with a rectangular cross section, and the values for $\sigma_{\max}$ , E, and $\gamma$ were selected to be conservative. The manoeuvring and hard-landing load limits represent a typical load limit of 2.5g with a safety factor of 1.5. The maximum deflection is just over 10% of the semispan, which is reasonable for a high-aspect-ratio wing. However, it will be shown that results are sensitive to changes in this parameter.

Figure 6. Example net-weight distribution for the Ikhana wing carrying 2000 lbf of fuel and a generic instrumentation pod.

Table 1 Example specifications for the Ikhana airframe

Wings with taper ratios near $R_{T}=0.4$ produce a nearly elliptic lift distribution with no aerodynamic or geometric twist^(58,59) . Therefore, we will use the elliptic lift distribution for the baseline design. The solver described in Section 4.1 predicts a wing-structure weight of 1008.4 lbf and induced drag of 54.040 lbf for the baseline no-pod configuration. The total weight is 8508.4 lbf, and the wing loading is 31.831. For the baseline pod configuration, the solver predicts a wing-structure weight of 1080.5 lbf, giving a total weight of 8580.5l bf and a wing loading of 32.101. The induced drag is 54.959 lbf. A summary of the results for the baseline design is included in Table 2.

Table 2 Example optimisation results for the NASA Ikhana airframe

5.1 Minimising Induced Drag

The lift distribution, wingspan and wing-structure weight that minimise induced drag were found using the framework from Fig. 5, in conjunction with the SciPy¹ implementation of the Sequential Least-Squares Programming (SLSQP) method⁽⁶⁰⁾. Using SLSQP, the non-linear constrained optimisation problem is cast as an approximate linear least squares problem around the initial design variables x. This problem is solved to give an update for the design variables Δx. The original problem is then recast as a linear least squares problem around the updated point x + Δx, and the process is repeated until Δx falls below a specified tolerance. Gradients for the objectives and constraints are calculated using finite differencing. For additional details, see Ref. (60).

The wing loading is fixed at 31.831 for the no-pod Ikhana configuration and at 32.101 for the pod configuration. The net weight for both configurations is fixed at $W_{n}=7500$ lbf. A spar-width constraint of $w/c\leq0.1$ is also imposed. The wingspan b and the Fourier coefficients $B_{n}$ that define the lift distribution are the design variables. For the results shown here, the Fourier series is truncated at n = 29.

The optimum lift distribution for each configuration is shown in Fig. 7, along with five reference lift distributions labelled a, b, c, d and e. Curve a is the elliptic lift distribution. Curve b is Prandtl’s 1933 lift distribution⁽³⁹⁾. Curves c and d are the optimum lift distributions found by Phillips et al.⁽⁴⁸⁾ for the stress- and deflection-limited designs, respectively, of a rectangular wing with fixed wing loading and the weight distribution given by Equations (9) and (10). Curve e is the optimum lift distribution found by Taylor and Hunsaker⁽⁴⁹⁾ for the stress-limited design of a tapered wing with fixed wing loading, the weight distribution given by Equations (9) and (10), and a taper ratio of R _T = 0.4. Additional optimisation results are summarised in Table 2. Note that in this study, we have fixed the taper ratio to R _T = 0.421 for all configurations. Therefore, the optimum solutions presented in Table 2 have a different root and tip chord than the baseline configuration.

Figure 7. Solutions for the lift distributions that minimise induced drag for the example no-pod and pod configurations of the NASA Ikhana airframe.

Figure 7 shows that the lift distributions that minimise induced drag for the no-pod and pod configurations are nearly identical, and both lift distributions are noticeably non-elliptic. Table 2 shows that the magnitude of the Fourier coefficients decreases rapidly as n increases. The same trend is shown in Refs. (49) and (51). Both lift distributions are primarily dominated by B ₃, with $B_{3}=-0.091066$ for the no-pod configuration and $B_{3}=-0.084530$ for the pod configuration. These values fall near the theoretical optimum $B_{3}=-0.059716$ for the deflection-limited design of a rectangular wing with fixed wing loading⁽⁴⁷⁾. Indeed, both optimum Ikahana designs are deflection limited.

The reader is reminded that, to obtain any of the lift distributions in Fig. 7, the wing must be twisted. For an unswept wing with any given planform shape, the twist distribution required to produce a desired lift distribution, specified by B _n , can be computed using the method shown by Phillips and Hunsaker⁽⁵⁰⁾. However, in this study we assume that the wing is correctly twisted to produce the desired lift distribution.

From Table 2, we see that for the no-pod configuration, using the optimum lift distribution allows an increase in wingspan of 18.31% and an increase in wing-structure weight of 97.21%, and results in a reduction in induced drag of 8.93% over the baseline no-pod configuration. For the pod configuration, the optimum lift distribution allows an increase in wingspan of 16.79%, an increase in wing-structure weight of 86.32% and a reduction in induced drag of 7.95% over the baseline pod configuration.

The wing-structure weight distributions for the baseline Ikhana designs and the optimum designs are shown in Fig. 8, along with their corresponding planforms. Although Ikhana has a non-rectangular planform and a weight distribution other than Equation (9), the optimum wing-structure weight for each configuration is just over 26% of the net weight. This agrees relatively well with the theoretical optimum wing-structure weight of $W_{s}=W_{n}/4$ ⁽⁴⁸⁾ for the deflection-limited design of a rectangular wing with the weight distribution given by Equation (9).

Figure 8. Wing-structure weight distributions and corresponding planforms for the baseline design and optimum design of the example no-pod configuration and pod configuration of the NASA Ikhana airframe.

Induced-drag contours around the optimum design for each example Ikhana configuration are shown in Fig. 9 as a function of the design variables b and $B_{3}$ . In reality, the lift distribution is a function of n Fourier coefficients, and the design space is more than n-dimensional. However, because the optimum lift distribution for each Ikhana configuration is dominated by $B_{3}$ , we approximate the lift distribution using $B_{3}$ alone. Note that the induced-drag contours are not smooth at low wingspans, since the wing design transitions from stress limited to deflection limited at a low wingspan for each Ikhana configuration.

Figure 9. Induced-drag contours for the example no-pod configuration and pod configuration of the NASA Ikhana airframe.

Figure 9 gives some insight into the relative influence of the wingspan, weight and lift distribution on the induced drag at different points in the design space. For example, for both Ikhana configurations, the induced drag is much more sensitive to changes in wingspan than it is to changes in lift distribution around the baseline design. Since the wing-structure weight typically increases as B ₃ and b increase, it is more advantageous to increase the wingspan and the weight than to decrease the weight by changing the lift distribution near the baseline design. On the other hand, there are regions in the design space where reducing the weight by changing the lift distribution gives a greater reduction in induced drag than changing the wingspan.

The characteristics of the design space depend on the wing configuration, and a figure like Fig. 9 can require more than 100,000 function evaluations. However, using the methods presented in this paper, Fig. 9 was produced in seconds. Understanding of the design space during early design phases can facilitate rapid conceptual optimisation and reveal important aspects of the design that cannot be easily seen using high-fidelity methods alone.

5.2 Sensitivity of Optimum Solution to Design Parameters

To illustrate the sensitivity of the optimum solutions in this paper to changes in design parameters, Fig. 10 shows the percent change in the minimum induced drag, optimum wingspan, optimum $B_{3}$ and optimum wing-structure weight as a function of the percent change in pod location, average S, and the parameters W _r , W/S and R_T for the pod configuration of the NASA Ikhana airframe. The percent change in pod location is measured in percent semispan.

Figure 10. Percent change in minimum induced drag and optimum wingspan, B ₃, and wing-structure weight with change in pod location and the parameters S _b , W _r , W/S and R _T for the example Ikhana pod configuration.

The plots for W _r and pod location in Fig. 10 show that the optimum lift distribution, characterised by B ₃, is most sensitive to the weight distribution. As the pod is shifted away from the wing root and the root weight decreases, the value for B ₃ also decreases. This corresponds to a less-elliptic lift distribution, which results in an increase in the wingspan and lower induced drag. This supports the result found by Phillips et al.⁽⁴⁷⁾ that the optimum root weight is given by Equation (10), which predicts that the theoretical optimum root weight for Ikhana is close to $W_{r}\approx3500$ lbf. The lift distribution is not sensitive to changes in average S _b or W/S, and B ₃ only changes by about $\pm$ 1% with $\pm$ 10% changes in R _T , which agrees with the observation made by Taylor and Hunsaker⁽⁴⁹⁾ that the optimum lift distribution is relatively insensitive to the taper ratio.

Figure 10 also shows that the wing-structure weight does not change with changes in average S _b and W/S, and it changes by less than ±0.65% with ±10% changes in pod location, W _r and R _T. This supports the analytic solutions found by Phillips et al.⁽⁴⁸⁾ and Taylor and Hunsaker⁽⁴⁹⁾ that the optimum wing-structure weight is independent of all other design parameters.

Only the optimum wingspan and corresponding induced drag are affected by changes in average S _b and W/S. For S _b, this is not surprising, since increasing S _b means that less weight is required to support the bending moments. This allows for larger increases in wingspan with smaller corresponding increases in wing-structure weight. For the range of S _b shown, the optimum design is deflection limited, which means that S _b is inversely proportional to $\gamma$ and directly proportional to $C_{\delta}$ , E and $\delta_{\max}$ , as shown in Equation (15). Therefore, the sensitivities shown in Fig. 10 for S _b are also characteristic of the sensitivities for $C_{\delta}$ , E, $\delta_{\max}$ and the quantity $1/\gamma$ .

The results in this section show how the methods presented in this paper can be used for design-space exploration. Because the methods are fast, they can be used to rapidly visualise the coupled aerostructural design space and obtain solution sensitivities to various design parameters. It should be remembered that the results shown here are only valid for the two example configurations of the NASA Ikhana airframe given in Table 1. Nevertheless, the methods presented in this paper can be used for any unswept planar wing with arbitrary planform and weight distribution to rapidly iterate on possible design concepts.

6.0 Conclusions

Low-fidelity methods are valuable for rapid aerostructural optimisation during the conceptual and preliminary design phases. However, most modern aerostructural methods use mid- and high-fidelity solvers, which are better suited for later design phases. The majority of analytic and low-fidelity aerostructural optimisation methods are limited in application to wings with specific planforms and weight distributions. Here, a low-fidelity numerical method has been presented that includes more general approximations corresponding to arbitrary planforms and weight distributions. The method uses an iterative solver to determine the wing-structure weight and induced drag for a given lift distribution and wingspan. The solver is used within an optimisation framework for rapid design-space exploration and optimisation.

Section 5.0 shows an example application of the method presented in this paper to two configurations of the NASA Ikhana airframe. A summary of the optimisation results, including the optimum wingspans, wing-structure weights and lift distributions, is given in Table 2. The optimum lift distributions for both Ikhana configurations are shown in Fig. 7. It has been shown that the optimum lift distributions for the Ikhana configurations are very similar to the analytic optimum lift distribution for a rectangular wing, with the ideal weight distribution given in Equation (9). The optimum wing-structure weight for each Ikhana configuration is also in good agreement with theoretical solutions.

A visualisation of the design space for each Ikhana configuration is shown in Fig. 9. The relative influence of the wingspan, lift distribution and wing-structure weight depends on the location of the design in the design space. Figure 10 shows the sensitivities of the design values around the optimum solution to changes in pod location, proportionality coefficient, root weight, wing loading, and taper ratio for the pod configuration of the Ikhana airframe. The optimum wingspan is most sensitive to the proportionality coefficient and wing loading, and the optimum lift distribution is most sensitive to the weight distribution. The optimum wing-structure weight is nearly independent of all other parameters. For the Ikhana configurations considered here, the optimum design allows a wingspan increase of up to 18.31%, an increase in wing-structure weight of up to 97.21% and a reduction in induced drag of up to 8.93% over the baseline Ikhana configuration. All results were obtained in a matter of seconds.

It should be remembered that the methods presented here were derived using the assumptions associated with lifting-line theory, including wing planarity, zero sweep and moderate to high aspect ratio. For other wing designs, modifications to these methods may be needed. However, the methods presented here are useful for many practical aircraft configurations. In early design phases, these methods can be used for rapid conceptual optimisation and visualisation of the design space. These results can provide important insight into the effects of the wing aerodynamic and structural properties and the wing weight distribution on the minimum-induced-drag design.