Category Flight Vehicle Aerodynamics

Lifting Surface Theory

Lifting surface theory is an extension of thin airfoil theory to 3D. It models the flow about the wings and tails of a general 3D aircraft configuration using vortex sheets y(s,£), or the equivalent normal-doublet sheets p(s, e). The objective is to represent the lift, sideforce, moments, and induced drag of the configuration using only the camber surface shapes, with the volume effects of the various components being ignored.

6.4.1 Vortex/doublet sheet geometry

The assumed geometry of the vortex or doublet sheets is shown in Figure 6.4. The sheets are assumed to be everywhere parallel to the x axis, with the camber-surface shapes of the actual geometry represented only by their normal vector distribution n(s, e). The sheet strengths 7(s, e) or p(s, e) are unknown only over the extent of the actual surface. On the trailing wake portions of the sheets, the strengths are constant in x, and equal to their trailing-edge values.

Подпись: Y(s) = Y(S,ITE) ■ x (on wake) or p,(s) = P(S,£TE) (on wake) (6.16)

These are also the sheet strengths in the Trefftz plane, as shown in Figure 6.4.

6.4.2 Lifting-surface problem formulation

Подпись: v-'"> = і Lifting Surface Theory Подпись: (6.17)

The perturbation velocity field of the vortex sheet distribution is given by the superposition integral (2.16),

Lifting Surface Theory

where r'(s,£) is the assumed vortex sheet geometry. The integral is over both the surface and wake vortex sheets. With an airmass which is still in the Earth frame, the total fluid velocity observed by a point r fixed in the body frame is then obtained by subtracting that point’s velocity Up as given by (6.3).

V(r) = V7 – (U + Oxr) (6.18)

Подпись: V(r) ■ n Подпись: VY Подпись: (U + О x r) Подпись: n = 0 Подпись: (6.19)

The flow-tangency boundary condition is then

for each r(s, e) surface point. The Kutta condition

Y(s, eTE) x X = 0 (6.20)

is also applied all along the trailing edge of each surface. Equations (6.17), (6.19), (6.20) together constitute an integral-equation problem for the unknown 7(s, e) distribution. If instead p(s, e) is chosen as the unknown variable, then 7 = П x Vp would be substituted into all the above expressions. This lifting-surface problem can be solved by the Vortex Lattice Method, described later in Section 6.5.

Quasi-steady force and moment parameterization

Dimensional analysis indicates that all the force and moment coefficients have the following parametric dependence for steady flows. The same dependencies also approximately hold for quasi-steady flows such as an aircraft in slow maneuver.

Cd = Cd(а, в, p, q,r, а, в,5t,5a,5e, Sr, M^,Beao)

Cy = Су (а, в, p, q, Г, а ,в, St, Sa, Se, Sr, MTO, Beoa)

Cn = Cn(a, в, p, q,r, а, в, St,5a,5e, Sr, MTO, Be00)

The parameters St, Sa, Se, Sr are throttle, aileron, elevator, rudder control parameters (there may be more). These represent thrust settings or control-surface deflection angles, which influence the overall force and moment on the aircraft.

Подпись: Figure 6.3: Aircraft а, в flow parameter space (remaining parameters p,r ,r ... not shown). CL (a) variation is shown for parameter slice A, and Cn (в) variation is shown for slices B, C. The linearized force or moment representations, indicated by dashed lines, are valid in a sufficiently small region near the chosen baseline operating point. Two baseline points 1,2 at two different a0 values are indicated, with the stability derivative CLa changing dramatically between the two points. Quasi-steady force and moment parameterization

Figure 6.3 shows typical flows over the entire possible а, в range, most of which involve large-scale separa­tion and flow reversal. The result is that the force and moment coefficients have complicated dependencies on the operating parameters, as indicated in the sample CL(a) and Cn(e) line plots in Figure 6.3.

The overall parameter space is enormously larger than the two-axis slice shown in Figure 6.3, since there are additional axes for p, q, r, a,3 ,ST…, MTO, Re. For many applications, such as mission performance estimates, stability and control analysis and design, etc., the force and moment coefficients only need to be defined within a small region of the parameter space, described by small deviations about some operating point or trim state, denoted by the ()0 subscript. Here the force and moment coefficient functions can be approximated by their linearized forms, or equivalently then – first-order Taylor series approximations.


^ Clo

+ CLa Aa

+ CLq



CLa A«

+ CLSe Д5е



– Cno

+ Cne Дв

+ Cnp

, ap +

Cnr A3

+ CnSr A$r



Clo =

– CL(ao, eo,.



= dCL/да





CLq =

: dCL/dq (ао, во,…)

Cno =

_ Cn(a0,e0v



= dCn/дв





Cnp =

: dCn/dp (ао, во,…)

■ ■ ■ (6.14)

The series variables are the aerodynamic parameter perturbations from the trim state.

Да = a — a0 Aft = в — во Д p = p — p0 … (6.15)

The series coefficients CLa, CLq… are stability derivatives, and Cls, Cngr … are control derivatives. These play a crucial role in aircraft flight dynamics and stability and control, as outlined in Chapter 9. Note that these coefficients can substantially depend on the baseline trim state values, and some may have their signs reversed between different trim states, as for example CLa shown in Figure 6.3.

Non-Dimensionalization and Parameterization

6.3.1 Dimensionless variables

Aerodynamic characteristics are almost invariably defined, described, or provided in terms of the following dimensionless force coefficients, moment coefficients, and rotation rates. The latter two are commonly used in either the body or the stability axes.

Подпись: CD = D G = c P bref Q/ Sref q/ Sref bref p = 21/ Gy = Y Cm = M а = Q Cref Q/ Sref Q/ Sref cref 21/ CL = L Gn = N r = Г bref Q/ Sref Q/ Sref bref 21/ (6.9)

All the reference quantities are arbitrary. The traditional choice for the reference area Sref is the projected wing area S. This typically includes any hidden “carry-through” wing area inside a fuselage, but may exclude root fairings or fillets. The choice for bref is the actual projected wingspan, which may or may not include tip devices such as winglets. The traditional choice for cref is the wing’s mean aerodynamic chord,

1 Ґ2 2

Cref — cmac = — / C(y) d у (6.10)

S – b/2

which is in effect a root-mean-square chord. Workable alternatives are the average chord cref = S/b, or simply the root chord cref = c(0).

Подпись: a Cref Подпись: а в bref в = Ж Подпись: (6.11)

Theoretically, additional important parameters are the dimensionless flow-angle rates.

This antifies the strength and influence of the wing’s shed vorticity, which is present in unsteady airfoil flows and is discussed in more detail in Section 7.4.2. In brief, a determines the time delay in the wing’s downwash seen by the horizontal tail, and therefore during pitching maneuvers it influences the time evolu­tion of the overall Cm, and to a lesser extent of the Cl also. In most aircraft в has relatively little influence and is usually ignored, although it may be significant for unusual aircraft configurations.

Axis Systems

In computational methods, the aircraft motion vectors U, О, and the aerodynamic force and moment vectors F, M are most easily specified or calculated in the same xyz axes which are used to specify the geometry itself, shown in Figure 6.1. To apply the results to aircraft performance, stability and control, and other related disciplines it is necessary to provide these vector quantities in other more relevant axes.

6.2.1 Stability axes

Axis Systems Axis Systems Axis Systems

The drag, sideforce, and lift force components are most commonly defined in the stability axes, which are rotated from the geometry axes by only the angle of attack a (not by sideslip в), as shown in Figure 6.2.

Figure 6.2: Orientation of stability-axis (freestream-aligned) aerodynamic force components D, Y,L, moment components Ls, Ms, Ns, and rotation rate components ps, qs, rs.

All the standard vector components in stability axes are defined from their components in geometry axes using the T rotation matrix.

Подпись: T s Подпись: Fx Fy Fz Подпись: ? Подпись: Ts Axis Systems Подпись: s Подпись: -Q2 Подпись: (6.5)
Axis Systems


Note that the moment and rotation-rate components, Ls, Ns and ps, rs, have reverse signs compared to the force components D, L. In effect, the stability axes used for the moments and rates are rotated by 180° about the y axis relative to the stability axes used for the forces. Note also that the T matrix leaves all the vector y components unchanged. Hence we have Y = Fy, Ms = My, qs = Qy.

6.2.2 Wind axes

The stability axes are not quite appropriate when examining the drag of an aircraft in sideslipping flight with

/1/0, since the drag I) as defined by the T matrix in (6.5) is not the hue stream wise drag force. In this

= W

Подпись: w Подпись: (6.6)
Axis Systems

situation we can invoke the wind axes, which are implemented by the rotation matrix T, which consists of a and в rotations, applied in that order as shown in Figure 6.1.

Подпись: Fx Fy Fz cos/ -sin/ 0 cosa 0 sina cos/ cosa —sin/ cos/ sina sin/ cos/ 0 0 1 0 = sin/ cosa cos/ sin/ sina (6.7) 0 0 1 _ —sina 0 cosa — sina 0 cosa w


Note that the drag as defined by the T matrix in (6.6) is exactly equivalent to the dot product of the total force and the unit freestream.

D = F ■ V//!/ (6.8)

= W = s

Furthermore, the sideforce Y produced by T in (6.6) is almost the same as that produced by T in (6.5), and the lift L is identical. Because the simple relation (6.8) is available to define the exact D when needed, and the exact Y is of relatively little importance, wind axes see little use in practice.

Aerodynamics of Aircraft in Maneuver

This chapter will examine the aerodynamics of thin wings of arbitrary planform and of slender bodies in arbitrary translation and rotation. Quasi-steady flow will be assumed.

6.1 Aircraft Motion Definition

Chapter 9 will derive in detail the Earth and body axis systems used for describing aircraft motion. Here, a few of those key relations will be simply stated without derivation. Unless otherwise indicated, all vector components will be assumed to be in the geometry axes shown in Figure 6.1, which have x and z reversed from the body axes given in Chapter 9. The other axis systems will be discussed where appropriate.

6.1.1 Aircraft velocity and rotation

Aerodynamics of Aircraft in Maneuver

The aircraft motion is defined by the velocity U of its axis-origin point, and by its rotation rate П. Both are shown in Figure 6.1. These are defined relative to the Earth frame, and hence they are also the velocity and rotation rate of the aircraft relative to a still airmass.

The aerodynamic “freestream” velocity V, is directly opposite to U, and is conventionally specified by the two aerodynamic flow angles a and в, applied in that order as shown in Figure 6.1.


Ux 1 ( — cos a cos в 1

Uy = —V, = V, I sin в > (6.1)

Uz — sin a cos в

Ко = /ux + Ц? + U.’} , a = arctan—K – , f3 = arct. an — ^4 (6.2)

V ‘T y ~ ’ – Ux ’ ^

Given these reciprocal relations, {К,, а, в} and {Ux, Uy, Uz} are equivalent alternative parameter sets. In practice, а, в are chosen as the independent parameters. These define the three components of the normal­ized aircraft velocity U/V, via (6.1), which are needed to compute the aerodynamic forces and moments.

6.1.2 Body-point velocity

The Earth-frame velocity of any point rp fixed on the body is given by

Up = U + Oxrp (6.3)

as shown in Figure 6.1. If the airmass is still (without wind or gusts), then the apparent airmass velocity seen by this point is —Up, which is in effect a “local freestream.” This will be used to formulate flow-tangency boundary conditions in computational methods.

Non-planar wake

Another interesting case is a non-planar wake, such as that produced by a wing with winglets, which were originally developed by Whitcomb [49]. A winglet acts much like a span extension in that both spread out the shed vorticity, which reduces the velocities and kinetic energy in the Trefftz plane, and thus reduce Di. A span increase does this more effectively than a winglet, but on the other hand a winglet produces a smaller increase in the root bending moment.

Figures 5.18 and 5.19 show two possible ways to parameterize the geometry of the wing+winglet combi­nation, and the resulting Di relative to the no-winglet case value Dil. Results both without and with the root bending moment constraint are shown. The bending moment constraint is seen to put a floor on the

Di/Di-і ratio at about 0.84, which is comparable to the minimum value of the best flat-wing case shown in Figures 5.16 and 5.17. The only apparent advantage of the winglet is that much of this benefit can be obtained with a smaller overall span. The conclusion is that winglets are effective mainly in cases where the overall span is limited by other than structural constraints. The relative merits of winglets and various other types of non-planar lifting surface systems is discussed by Kroo [47].

Non-planar wake

Non-planar wake

Lift fixed


Non-planar wake

0 10 20 30 40 50 60 70 80 90

^ winglet і de9 ]

Figure 5.18: Induced drag of wing + winglet with fixed inner span, versus winglet height and angle above horizontal. Lift is the same for all cases. Plot on right in addition has a fixed root bending moment. The no-winglet case with elliptical loading provides the reference value Dil, and also the fixed lift and bending moment values.


0 10 20 30 40 50 60 70 80 90

^winglet і de9 ]




Fixed total span


z winglet


Non-planar wake

Non-planar wakeNon-planar wakeNon-planar wake

Example optimum load distributions Planar wake

As shown in Section 5.8.1, the elliptical load distribution is optimum for the case of a planar wake with a fixed span. However, in many aircraft applications a more relevant constraint is not on the span but on the root bending moment, since this dominates the wing’s structural weight which offsets induced drag reductions.

To illustrate this tradeoff, Figure 5.16 shows three load distributions for three different specified spans, each having the same lift and root bending moment. Increasing the span reduces Di significantly, even though the resulting load distributions are very “sub-optimal” in a fixed-span sense. Referring to the Fourier Di expression (5.71), the increase in the span b more than overcomes the increased parameter S which measures how much the loading deviates from elliptical.

Figure 5.17 shows the relative Di for a continuous range of spans for the flat-wing case with constrained bending moment. With no bending moment constraint, the optimum loading is elliptical for any span, and the induced drag then scales simply as Di ~ 1/b2, indicated by the thin line in Figure 5.17.

Example optimum load distributions Planar wake

0 0.2 0.4 0.6 0.8 1 1.2 1.4

2y / b1

Figure 5.16: Optimum Ap(y) load distributions and corresponding w(y) normal velocity distribu­tions for three specified spans, all with the same specified lift and root bending moment. The baseline elliptic-loading case has span bi and induced drag Dil.

Example optimum load distributions Planar wake

Figure 5.17: Induced drag versus span, with fixed lift and root bending moment, relative to base­line case. The three symbols correspond to the three cases in Figure 5.16. Thin line is the result Di/Dii = (b/b1)-2 for elliptical loading for all spans, which is optimal in the absence of abending moment constraint.

Additional constraints

Any number of other constraints can be added in addition to the lift, such as the root bending moment mentioned earlier. An effective general solution technique here is define a Lagrangian function L, which is the objective function plus all the constraints,

L(a^(s) , Л1 , Л2 …) = Di + Л1 (L — Lspec) + Л2(М0 — M0spec) + … (5.102)

where the Лі,2… coefficients are Lagrange multipliers associated with the constraints, and are solved as part of the problem. The constant Л in the previous section was in fact the same as Л1 here.

The constrained-optimum solution Ap(s), Л1, Л2 is defined by the requirement that L is stationary.

5C = 5Di + 5Лі (L-LSpec) + ^2 (Mo-Mospec) + Лі 5L + Л2 SMo = 0 (5.103)

P* A Sp

y cos в + z sin в P^V* Ap cos в ds – Lspec

Подпись: + Ai cos в + Л2 —К dn 2
Подпись: + SЛl + SЛ 2
Подпись: I Ap
Подпись: y cos в + z sin в
Подпись: ds - M0
Подпись: ds
Подпись: spec
Подпись: 0
Подпись: (5.104)

Substituting for 5Di, SL, 5M0, L, M0, and collecting terms having the same A Sp, 5Л1,5Л2 factors gives


У Aij Apj – Л1 V* cos ві j=1

P*V*y A pi cos ві Asi

Подпись: A2-hL Additional constraints Подпись: 0
Подпись: (i = 1 ...N) (5.105)
Подпись: L
Additional constraints

which for optimality must be zero for any A6p(s), dA , dA2. This requirement is met by setting all the quantities in the brackets to zero, using the wake panel method to discretize the integrals. Since the first bracket is inside the integral (and inside the equivalent discrete sum), it must be set to zero at each of the N discrete panel points. In contrast, the second and third brackets set to zero are single equations. The result is the following (N+2) x (N+2) linear system for A pi, Л1, Л2.

Optimum potential jump calculation

The above application of Calculus of Variations produced the optimum normal velocity distribution dp/dn(s). The one remaining step is to determine the corresponding Ap(s). A suitable numerical approach is to use the 2D panel method sketched in Figure 5.12. Using the AIC matrix Aj defined by (5.81), condition (5.98) is imposed at each panel control point.


^ Aij Apj — Л cos ві = 0 (i = 1 …N) (5.100)


The constant Л is one additional unknown in the problem. The appropriate additional equation is the specified-lift constraint, written by using the discrete lift expression (5.78).


P^VL^2 a pi cos ві Аві = Lspec (5.101)

i= 1

Equations (5.100) and (5.101) together constitute a (N+1) x(N+1) linear system for the unknowns Api, Л. The corresponding induced drag can then be computed by re-using the AIC matrix Aij to obtain dp/dni, and then using this in the discrete induced drag expression (5.82).

Fuselage wake contraction effect

For a configuration with a fuselage of significant size, such as the one shown in Figure 5.5, the average sheet velocity Va can no longer be assumed to be parallel to the freestream, so the wake does not trail straight back from the wing trailing edge. The actual velocities and streamline trajectories can be determined from a panel or slender-body model of the fuselage (see Section 6.6). Nikolski [48] used instead a simple axisymmetric fuselage flow model shown in Figure 5.13, where Va is assumed to be parallel to axisymmetric streamtubes. Conservation of mass between the streamtube cross-section at the wing and in the Trefftz Plane gives

m(y) = pTOV* n (y2 – (d/2)2) = pTO V* ny2 (5.83)

y(y) = sjy2 + (d/2)2 (5.84)

which assumes that the mass flux pV magnitudes adjacent to the fuselage of maximum diameter d are nearly the same as in the freestream.

Fuselage wake contraction effect

Figure 5.13: Axisymmetric streamtube flow around fuselage determines wake contraction from wing to the Trefftz plane.

Equation (5.84) is the correspondence function which specifies the wing location y which is connected to wake location y by an average streamline. For a given wing circulation distribution Г(у), the potential jump in the wake is then given in terms of the correspondence function.

A p(y) = r(y(y)) (5.85)

As an example, consider the case of an elliptical spanwise loading in the wake.

A <p(y) = A P’0 J 1 – (2 у/b)2 = Ap0 sin d (5.86)

b2 = b2 – d2 (5.87)

Fuselage wake contraction effect Fuselage wake contraction effect
Подпись: 2
Подпись: 1
Подпись: ds
Подпись: MQ
Подпись: (5.92) (5.93)

Only the lift constraint will be assumed first to simplify the initial discussion of the concepts. Adding other constraints will then be considered.

which holds for any two fields f, g which satisfy V2/ = 0 and V2g = 0. For our case we choose / = p and g = 5p, in which case the identity shows that the two terms in the equation (5.96) integrand are actually equal. Hence, omitting the second term and doubling the first term will not change the result.

SDi = – p^J ASip^ ds (5.97)

If Di is to be a minimum it’s necessary that it be stationary, specifically that 5Di = 0, for any admissible A 5p(s) distribution along the sheet. This is satisfied by a normal velocity distribution which is everywhere proportional to the local cos в,

—^-(s) = Л cos 9(s) (5.98)


where Л is some constant. This solution can be verified by putting it into the £Di expression (5.97), to give

5Di = —poo f А-Sip A cos 6 ds = —SL = 0 (5.99)

as required. The conclusion is that a normal velocity which is given by (5.98) results in the smallest possible induced drag for a given lift and a given wake shape. This result is exactly consistent with the result (5.76) obtained via the Fourier series approach for the flat wake case. The great advantage of (5.98) is that it applies to a wake of any shape.