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

2.0 ANALYTICAL FOUNDATION

Using Prandtl’s classical lifting-line theory^{(Reference Prandtl9,Reference Prandtl10)} , the dimensionless spanwise section-lift distribution on a finite wing with no dihedral or sweep immersed in a uniform flow can be written as^{(Reference Phillips, Hunsaker and Joo47)}

(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^{(Reference Phillips and Hunsaker50)}. 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^{(Reference Prandtl39)} 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^{(Reference Prandtl39)}. 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^{(Reference Prandtl39)}.

Phillips et al.^{(Reference Phillips, Hunsaker and Joo47,Reference Phillips, Hunsaker and Taylor48)} 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^{(Reference Phillips, Hunsaker and Joo47)}

(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^{(Reference Phillips, Hunsaker and Joo47)}

(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. (Reference Phillips, Hunsaker and Joo47). For deflection-limited designs, Equation (5) can be rewritten as^{(Reference Phillips, Hunsaker and Joo47)}

(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^{(Reference Phillips, Hunsaker and Joo47)}

(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.^{(Reference Phillips, Hunsaker and Taylor48)} 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.^{(Reference Phillips, Hunsaker and Joo47,Reference Phillips, Hunsaker and Taylor48)} 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.^{(Reference Taylor, Hunsaker and Joo51)} for the stress-limited design.

For deflection-limited designs, the vertical spar deflection can be found using the relation^{(Reference Phillips, Hunsaker and Joo47)}

(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.^{(Reference Taylor, Hunsaker and Joo51)}, but here we will include the deflection-limited design and several additional constraints that were not considered in Ref. (Reference Taylor, Hunsaker and Joo51).

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^{(Reference Taylor, Hunsaker and Joo51)}

(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.^{(Reference Taylor, Hunsaker and Joo51)}

(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^{(Reference Richardson52)} 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^{(Reference Taylor and Hunsaker49)} 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. (Reference Taylor and Hunsaker49).

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^{(Reference Phillips, Hunsaker and Joo47,Reference Phillips, Hunsaker and Taylor48)} . In this study, we will consider only a few example constraints proposed by Phillips et al.^{(Reference Phillips, Hunsaker and Joo47,Reference Phillips, Hunsaker and Taylor48)} , 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. (Reference Taylor and Hunsaker53), 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. (Reference Taylor and Hunsaker53).