On the teaching of the principles of wing

flexure-torsion flutter
G. J. HANCOCK and J. R. WRIGHT
Departmentof Aeronautical Engineering,
Queen Mary College, University of London
and
A. SIMPSON
Departmentof Aeronautical Engineering
University of Bristol

because there is insufficient time to discuss unsteady aero-

SUMMARY dynamics, although compatibility with standard quasi-steady
aerodynamics should be noted. More formally, flutter is the
This note attempts to clarify some basic ideas, and correct some condition when an aircraft component (eg wing or tailplane)
common misunderstandings, on wing flexure-torsion flutter. exhibits self-sustained oscillations at .a certain forward speed.
Topics include: That forward speed is the critical flutter speed. At speeds
below the critical flutter speed any initial dynamic structural
(i) an understanding of the physics of wing flexure-torsion vibrations will be damped whereas at speed above the critical
flutter through a number of specific examples, flutter speed any initial dynamic structural disturbance will
(ii) an examination of the energy interpretation of the physics grow, leading, if not limited by non-linear aerodynamic or
of flutter, indicating its limited usefulness, structural behaviour, to structural failure. It is mandatory that
(iii) a discussion of the significant differences between wing the critical flutter speed be outside the operational flight
flexure-torsion flutter and the Duncan flutter engine, envelope of an aircraft with a sufficient safety margin (usually
(iv) a graphical representation to assess the contribution of at least 15%), demonstrated by a combination of estimation
various parameter to flutter onset, and flight testing procedures. Ground resonance testing (also
(v) subcritical response below the critical flutter speed, mandatory in UK) provides a preliminary check on the mathe-
(vi) some practical wing flutter considerations. matical model of the structure.
Broadbent wrote an excellent series of articles on elemen-
tary concepts in aeroelasticity in 1954 which were published in
Aircraft Engineering. These articles were also published
1. INTRODUCTION collectively as a monograph which unfortunately has been out
of print for many years. Nevertheless, these articles are
The aim of this paper is to outline a few basic ideas on wing essential reading for teachers of aeroelasticity. (In passing, the
flutter for incorporation into an undergraduate degree course, past volumes of Aircraft Engineering are well worth browsing
through, many articles, particularly on aircraft design, are still
recognising that the amount of time for the teaching of aero- relevant today and make stimulating reading.)
elasticity is strictly limited. Text books are not altogether
helpful in clarifying the physics of flutter. Furthermore, as will This note is complementary to the articles by Broadbent.
be pointed out, superficial reading of text books has lead to a In Sections 3, 4 and 5 of this paper the critical flutter
number of misconceptions. conditions for a number of examples of idealised two degree of
Before introducing students to the concepts of flutter, a freedom systems (ie wing bending and wing torsion) are
background knowledge of conservative structural dynamics is described. It is explained how the aerodynamic loads couple
required. In particular, for a continuous aircraft structure, it is the two degrees of freedom to give instability. Although the
necessary to understand the mathematical modelling, associ- analysis employed involves simple mathematical modelling,
ated approximations and numerical procedures which lead to the principal aim is to build up physical insight.
estimates of structural normal modes and natural frequencies. The arguments depend only on the understanding of the
To complement the mathematical modelling, some knowledge characteristics of the standard linear differential equation for
of vibration testing and the measurement of normal modes is a single response variable,
desirable.
Flutter is a complex phenomenon where, in the classical
sense of the term, two or more structural normal modes are
coupled and excited through time dependent aerodynamic m^+cft + kz = F(t) . . . (1)
loads. The form and character of time dependent aerodynamic
loadings have to be taken on trust at undergraduate level
where m, the inertia or mass coefficient, c the viscous damping
Paper No 1312. Reprints available — see page 305. coefficient, and k, the stiffness coefficient, are all constant,

Aeronautical Journal October 7985 Hancocketal 285

 independent of time. In equation (1) when all of the coeffi- EI(y) spanwise bending stiffness distribution
cients are positive, ie F(t) generalised force input
f„(x/c, yls) mode shape corresponding to generalised
m>0, d>0, k>0 co-ordinate
GJ(y) spanwise torsional stiffness distribution
the response z(t) is stable in the sense that an initial distur- generalised displacement
bance, following an impulse (F(t) = Fa8(t)), decays to zero moment of inertia in bending
asymptotically with time. When c is positive and small (ie c2 < moment of inertia in torsion
4mk), the response is a damped oscillation, whereas when d is Iy» product of inertia
positive and large (ie c2 > 4mk) the response is a pair of k stiffness
non-oscillatory subsidences. ky bending moment stiffness
When the mass and stiffness coefficients are positive, but the ke torsion moment stiffness
damping coefficient is negative, ie Kmn generalised stiffness
Ku K2 constants
m > 0, c < 0, k > 0, Lz(v), ^i(v) X two dimensional lift oscillatory derivatives for
the system is unstable; the response due to an initial distur- Wv), W?) ) arbitraryflexuralaxis
bance grows indefinitely with time either as a divergent LM, W)\ two dimensional lift oscillatory derivatives
L,(y), i4r)\ whenflexuralaxis is on VA chord axis
oscillation when (c2 < 4mk) or as a pure non-oscillatory m j
mass
divergence when (c2 > 4mk). Mmn generalised mass
When the mass and damping coefficients are positive, but non-dimensional aerodynamic torsional
the stiffness coefficient is negative, ie Me
damping derivative
My(v) M-(v) } aerodynamic torsional oscillatory derivatives
m>0,c>0,*<0, Me(v) Me{v) ( about flexural axis
the system is unstable; the response due to an initial distur- *fr(y) M-(v) aerodynamic torsional derivative where
bance grows with time as a pure non-oscillatory divergence. Me(v) Aft(v) flexural axis is on VA chord axis
Critical conditions arise therefore when P power output
PL power load demand
Q Vzp^Sca
0 1n(t) generalised co-ordinates
In complex amplitude when qn(t) = ?„£""'
which implies an undamped simple harmonic oscillation at S
frequency o>, where semi span
S (half) wing area (= sc)
T torsional moment
w2 = k/m; V free stream velocity
vF critical free streamflutterspeed
and when w(x,y (downward) displacement function over a
'')
wing surface
k = 0, work done
ft distance of flexural axis from wing leading
Xf(=x- f)
which implies a state of neutral equilibrium in which z(t) is edge
independent of time, retaining any initial value, z(o), given to distance of centre of mass from wing leading
it. *cm \ ScmC)
edge
Section 6 deals with the energy approach of understanding X co-ordinate along chord from leading edge to
flutter, it is argued that this approach has limited usefulness. trailing edge
In Section 7 the difference between wing flexure-torsion spanwise co-ordinate
flutter and the well known Duncan flutter engine is explained. y
z{i) displacement variable
In Section 8 a graphical method is outlined which shows the 7 downward non-dimensional displacement of
parametric contributions to the flutter condition. flexural axis
In Section 9 a description is given of the response characteris- Zy„ downward non-dimensional displacement of
tics of a binary wingflexure-torsionsystem at speeds below the V4 chord axis
flutter speed (so called sub-critical response). P air density
An outline is given of the formulation and solution of more y bending angle (+ve wing tip down)
general binary flutter equations in Section 10. y complex amplitude when y(t) = ye*"*
Finally, in Section 11 some aspects of practical wing fluttei ya amplitude taking y(f) = y0 sin tot
calculations are summarised. 0 torsional angle (+ve nose up)
e0o complex amplitude when 0(t) = He*"
amplitude j w h e n e{f) = e s i n (fttf + ^
2. NOTATION phase difference (
e(0n generalised torsional displacements
a two dimensional lift curve slope frequency, rad/s
A generalised aerodynamic loads °>y natural bending frequency = {kylly)yi
B bending moment <0e natural torsion frequency = (kg/Ig)''2
By{v), By{v) aerodynamic bending oscillatory derivatives <i)
F frequency at critical flutter speed
B»(v), Be(y) V (= a)c/v), non-dimensional frequency
c chord parameter
c structural damping M damping
ec distance of flexural axis aft of aerodynamic € strain field
centre <T stress field

286 Hancock eta/ Aeronautical Journal October 1985

 motion in an air stream is always highly damped with the
3. SINGLE DEGREE OF FREEDOM MOTIONS damping increasing in proportion to V, the stream velocity, and
decreasing with altitude. At high speed, the (aerodynamic)
3.1 Wing Bending damping can be over-critical.
Consider a rigid non-swept, rectangular wing (chord c, semi-
span s) hinged at its root to allow a bending degree of freedom 3.2 Wing torsion
only, as shown in Fig. 1. This arrangement literally represents a
helicopter rotor blade attachment, but it is intended to idealise Consider the same rectangular wing hinged about a flexural
the bending characteristics of a finite wing. The spring axis to allow the single degree of freedom of torsion, as shown
constant, ky, represents the flexural stiffness of die wing in its in Fig. 2. The spring constant ke represents the fundamental
fundamental mode, and the inertia, Iy, the associated torsion-mode stiffness and Ie the corresponding moment of
(generalised) inertia. inertia about the flexural axis.
In still air, when structural damping is ignored, the equation In still air, again ignoring structural damping, the equation of
of bending motion following an initial disturbance is torsional motion is

ly>y + kyy = 0, . . . (2) Ied + ked = 0, (5)

where the bending moment of inertia Iy is given by where

Iy= I 8my2, I„ = J8m(x-xf)2 (6)

wing

and ky is the bending moment stiffness, as indicated in Fig. 1. and ke is the torsional stiffness. The still air natural torsional
The still air natural bending frequency wy is given by frequency o), is given by

ov ' ky/Iy. (3)

(7)
In a relative air stream of velocity V, the equation of bending
motion, following an initial disturbance becomes, on the basis
of quasi-steady aerodynamic strip theory, In a relative air stream of velocity V, a quasi-steady
aerodynamic lift force acts at the aerodynamic centre which is
Iyy + k/y = aerodynamic bending moment in direction of y located a distance ec ahead of the flexural axis, as shown in Fig.
increasing 3. A symmetric aerofoil profile is assumed so there is no steady
moment about the aerodynamic centre. The equation of
torsion motion becomes, by reference to Fig. 3,
= -Jv2pV*cdya(j£\)
O

(4) y
7\ " K
T y FLEXURAL AXIS
where TORSION SPRINS
STIFFNESS k .

S = wing area = sc, *•& (POSITIVE NOSE UP)

a = two-dimensional sectional lift curve slope. X
,—yC^^FLEXUR AL AXIS

The inertia, damping and stiffness coefficients in equation > **7 ^^

(4) are all positive. The aerodynamic loads in equation (4)
contribute to the damping only. Any oscillatory bending

Figure 2. Wing with single degree of freedom in wing torsion.

BENDING
MOMENT
STIFFNESS k BENDING ANGLE
(+VE WING TIP DOWN)
FIXED PIVOT

Figure 1. Wing with single degree of freedom in wing bending. Figure 3. Notation for aerodynamic loading.

Aeronautical Journal October 1985 Hancock et a I 287

 s

Ie0 + ke9 = jVip^cdy a dec 05 Ml 1-5 —— V

= VipV2Scae9 . . . (8) -0-5- y^

At a forward speed V, there is no aerodynamic damping, the / ROTATION ABOUT MID-CHORD

aerodynamic loads contribute to the effective stiffness only. -10-

/
/
MOMENT ABOUT MID-CHORD

The wing therefore oscillates in torsion in simple harmonic

motion at frequency o> where

Figure 4. Variation of M$ with frequency parameter.

h
So the correct quasi-steady limit on the right hand side of
It is seen that the frequency decreases with dynamic pressure
equation ( l l ) a s i ' — » 0 i s identical to the right hand side of
VipV2 when e is positive and that the motion becomes unstable equation (8).
as a pure divergence when
For a two dimensional aerofoil oscillating in pitch in simple
1 harmonic motion in an incompressible stream the variation of
/2pV2Scae>ke.
Mg with v, according to oscillatory aerofoil theory, when the
flexural axis is at the mid-chord location, is shown in Fig. 4.
The well known wing divergence speed VD is the critical
Now the 'standard' quasi-steady value of Me as quoted in text
condition of zero effective stiffness given by
books for a flexural axis at mid-chord is actually zero. It is seen
that this quoted quasi-steady value for Me bears no relation to
VipVDzScae = k e. . . . (10) the true value of Mg.
Confusion is introduced because the aeroelasticity derivative
Me is not the same as the aircraft stability and control derivative
In reality unsteady aerodynamic effects introduce positive Mg which does have a quasi-steady limit (given by an effective
aerodynamic damping so equation (8) becomes quasi-steady camber) effectively
—
(™e)aeroelasticity l ™ # "*" -™d)aircraft stability
Ifi + ke0 = Vip^Sc Led + M-J^\ . . . (11)
and it is the M& term in aircraft stability which has no (two
dimensional) quasi-steady limit.
By dimensional analysis it is seen that Me is dimensionless; Me is
To solve equation (11), since Me is a function of the
usually negative so equation (11) implies damped oscillations
frequency of the motion, an iterative procedure is required. To
below the wing divergence speed.
start the calculations, the effect of the damping term (Me) can
It should be pointed out that the right hand side of equation
be neglected, thus at a given speed V a first estimate of the
(11) is not strictly a correct representation of unsteady aero-
frequency parameter is
dynamic loads, however equation (11) has sufficient validity
for the present qualitative purposes. _o>c _ c {ke — Vip^Scae) Vz
The understanding of Me is one of the common misconcep-
tions referred to in the first paragraph of the Introduction. Text
books suggest that there is a quasi-steady value for Me based on
Then the value of Me corresponding to the above value of v can
the argument that 9 introduces an effective quasi-steady
be obtained from Fig. 4. The solution of equation (11) with this
camber which can be calculated by steady aerodynamic theory.
value of Mg will change the frequency; the procedure can be
This idea is totally erroneous.
repeated until convergence of v.
Quasi-steady aerodynamics are correctly defined as the limit
It is important to realise that the above converged answer is
of oscillatory aerodynamics as the frequency tends to zero, (ie
approximate, not 'exact'. The actual motion will be a damped
as frequency parameter v (= utcjV) —» 0). Now in two
oscillation, but the value of Me taken from Fig. 4 has been
dimensional oscillatory aerofoil theory for small v
derived on the assumption that the motion is an undamped
simple harmonic motion thus the effect of the damping on Me
Me(v) = 0(inv) . . . (12) has been neglected. For the present qualitative purposes, this
effect can be ignored but for practical quantitative calculations
so such effects need to be incorporated, especially for highly
damped motions.
Afj(i>)->°°asi'->0. . . . (13)

Oscillatory tests in a wind tunnel confirm this behaviour. Thus

a quasi-steady limit for M0 does not exist in two dimensions. An
4. COUPLED MOTIONS
explanation of where the text book explanation goes wrong is given
in Appendix 1. In three dimensions with finite wing theory Me
does in fact tend to a finite value as v —> 0, but, this value is not Consider the rectangular wing with the two degrees of freedom
related in any way to a camber effect. of bending and torsion, combining the arrangements shown in
Figs. 1 and 2.
The coupled equations of motion in still air are
It should be appreciated that in the limit as v—> 0, although in
two dimensions Me(y) —»°°, the product 9Me{v) -* 0, since Iyy + Iye9 + kyy = 0)
. . . (15)
&Mi(v)~e(v€nv)^>0asv-+0. . . . (14) Iyey + Ifi + ke0 = 0>

288 Hancocketal Aeronautical Journal October 1985

 where The aim of the following sections is to build up a qualitative
understanding of flutter through inspection of the equations
Iyg = product of inertia = / 8m(x — xf)y for a number of special cases.
wing

= m{xcm-xf)ycm • (16) Example 1

where m is the mass of the wing and (JCOT, yOT) are the co- Consider the example where the wing is mass balanced about
ordinates of wing centre of mass. Note that IyS is positive when its flexural axis (ie, Iy6 = 0) with the flexural axis passing
the centre of mass lies aft of the flexural axis. through the aerodynamic centre (ie e = 0); this example was
The still air coupled natural frequencies are given by close to the case on Spitfire Mklll. Equations (18) become

(- iy + ky){- w + k^ - IVM = o. (17)

Iyy + k^y = - Vip^Ss • (19(0)
One natural frequency will be somewhat higher than the
torsional uncoupled natural frequency (keIIg)Vl (the mode will
be primarily torsion with some inertial cross coupling with
bending) while the other natural frequency will be somewhat I ee + k e 6 = Vzp^Sc (19(H))
lower than the bending uncoupled natural frequency (ky/Iy)V2
(the mode is primarily bending with some inertial cross
coupling with torsion). The torsion equation is now uncoupled from the bending
In a relative air stream of velocity V, equation. The torsional motion is always damped because of
the Mg term, thus as time increases, 0—* 0 asymptotically. As $
Iyy + Iye6 + kyy = aerodynamic bending moment —* 0 in equation (19(i)) then the bending motion will decay
because of the positive aerodynamic damping, associated with
s
y. Sustained simple harmonic oscillations cannot therefore
= - \VzpV2 cdya Up1 + ojy occur, so flutter cannot occur. Furthermore, Hi this example,
the wing divergence speed is infinite (since e = 0; see equation
O (10)).

= - VipVSsa ( £ + ! ) • • • (18(i))
Example 2
Consider next the example where the wing is mass balanced, (ie
ly9y + Ie0 + ke0 = aerodynamic torsional moment
Iye = 0), the torsional damping term Me(0c/V) is neglected, but
e > 0 (ie flexural axis lies aft of the aerodynamic centre).
The equations of motion become
= I VipV1 c dy a (^ + 0jec
O

Iyy+kyy=-V2PV^Ssa^+0fy . . . (20(i))
2
+ VipV Sc Mg6c/V

IgO +kg0 = VipV^Scae (^ + d\. . . . (20(H))

= %p^5c[«(^+fl>)+^j.
For self-sustained simple harmonic motion to occur it is
(18(H)) essential, when taking into account the fact that equations
(20(i)) and (20(H)) are intrinsically coupled that,
(i) in the bending equation (20(i)), there shall be no
5. PHYSICS OF FLUTTER damping, that is, there shall be no effective y term,
(ii) in the torsion equation (20(H)), there shall be no
The equations of motion (equations (18)) represent the damping, that is, there shall be no effective 6 term,
coupling of the degrees of freedom of bending and torsion. At (Hi) the frequencies of bending and torsion motions shall be
any dynamic pressure (V-zpV2) there will be two independent identical.
coupled modes; each mode will have a motion of the form Consider the torsion equation (20(H)). For zero damping the y
gOi+K")^ with a frequency <o and damping p (p. < 0 gives positive term must not introduce a contribution to an effective 0 term,
damping), an amplitude ratio \y/0\ and a phase angle arg(y/0). It this can only occur if,
is important to realise that for each of these modes, equations
18(i) and 18(ii) are simultaneously satisfied with the same
the bending angular velocity y is in phase (or in antiphase)
values of co, pi, \y/d\ and aig(y/0).
with the torsional displacement 8 . . . (21)
The critical flutter condition occurs when the damping in one
of the modes becomes zero; this occurs at a critical flutter speed But if condition (21) is satisfied then 6 is Hi phase with y, so the
VF, at the coupled flutter frequency <oF, in a flutter mode, expression ys/3V + 6/2) in the bending equation (20(i)), is in
defined by \y/6\ and arg(-y/0). phase with y. Thus there is an effective damping Hi bending
Thus at the critical flutter condition, the effective dampings unless
in each of the separate motions in equations (18) must be
simultaneously zero and the frequencies of both motions (ie
equation (18(i)) and equation (18(h))) must be identical.
^ +4=0
3V 2 " U • (22)

Aeronautical Journal October 1985 Hancock etal 289

 When equation (22) is satisfied the equations of motion, As seen from equation (26), this type of flutter can only
equations (20(i)), (ii)), reduce to occur if the uncoupled frequencies are not too far apart. If a
purely two-dimensional problem had been considered in heave
Ly + kyy = 0 (23(0) and pitch, the uncoupled bending and torsion frequencies
would have to be identical for this type of flutter; this is the
reason that the above model is three-dimensional rather than
Ie6 +ke0 = V2pV2Scae- • • (23(ii)) the more common two-dimensional model.
However this kind of flutter cannot occur in practice
For the bending and torsion modes, as expressed by because the unsteady aerodynamic torsional damping term
equations (23(i)), (23(H)), to be coupled together at a common (Mg) completely changes the character of the flutter
frequency (in accordance with (ii) above), then mechanism.
Doubt is sown that all is not right with the above type of
VipVtfScaelA flutter when the reader is informed that at speeds below the
2 \^_ _ * above critical speed VF as given by equation (25) the system is
(Op = • (24)
unstable, whereas the speeds above the critical speed VF the
system is stable (this point is demonstrated later). So the
where VF (the critical flutter speed) is the unique value of V at example described in this section, and which is often high-
which the common frequency wF (the associated flutter lighted in undergraduate lectures, although of mathematical
frequency) occurs. Both are given by equation (24). Thus: interest, has no practical significance since the actual flutter
speed is zero airspeed.
the flutter frequency <uF is the same as the natural bending
frequency (from the first of equations (24))
Example 3
the critical flutter speed is given by the second equation of
(24) namely Consider the previous example with the wing mass balanced (ie
IyB = 0), with the damping term (— M^) now retained, but
kyle assuming that (— M-e) is small.
VipV/Scae/A V, (25) The equations are

bending is 90° out of phase (ie in quadrature) with torsion Vip^Ss a

(condition (21)).
Ly +kyy= -
[W+2) (27(i))

the amplitude ratio, from (equation (22)), is

+e
\y/6\ = 3/2(o>,*/V).
Ied + ke0 = Vip^Sc
[ae{w ) w AC

These are the four conditions which must be satisfied if flutter is • • • (27(fi))
to occur in the present case. Applying the same arguments as those applied to equations
For the flutter speed to be less than the wing divergence (20) in the previous example to equations (27), the (M-e) term
speed in this example would remain in the torsion equation damping out any
torsional vibrations and eliminating any possibility of sustained
oscillations. So that type of flutter is not possible. It is back to
(26) square one.
Consider the torsional equation (27(ii)). The damping in the
torsional motion, equation (27(ii)), would be eliminated if
This type of flutter is the one usually described in under-
graduate courses, depicted pictorially in Fig. 5.
yy + Mi)—=0. . . . (28)
In the motion shown in Fig. 5, the quasi-steady aerodynamic
bending moments due to y and 0 cancel out at each instant of
time, leaving the wing to oscillate in bending under its bending Equation (28) implies that
inertia and stiffness at its natural bending frequency; the air (i) y is in phase (or antiphase) with 0 . . . (29(i))
speed is such that the effective torsional stiffness, (ie structural
stiffness plus aerodynamic stiffness) in the torsion mode gives a (ii) |y/0| is small, assuming \2Mgc/aes\ to be small when \M-e\
torsional frequency equal to the natural bending frequency. is small . . . (29(ii))

. L I F T DUE : DUE
TO 8

L I F T DUE LIFT DUE

TO '& TO 8

Figure 5. An oscillatory motion (but not flutter onset).

290 Hancock eta/ Aeronautical Journal October 1985

 Turning to the bending equation,(27(i)), if 0 is in phase with y Example 4
(condition 29(i)), then it is literally not possible to eliminate the
y term in equation (27(i)). But if 0 is large compared with y, The next example is an extension of the previous example
from condition (29(H)), then the 0 term on the right hand side where again the wing is mass balanced about the flexural axis
of the bending equation (27(i)) is dominant over the y term and (ie Igy, = 0) but now the damping derivative (—Me) is not
so only a small phase difference between y and 0 would assumed to be small).
eliminate the damping, (this conjecture is proved in the next
example). The equations of motion are
Thus conditions (28) and (29) imply that equations (27)
reduce to,
Iyy + kyy=-V2pV>Ssa^ + <?) (33(i))

Iyy + kyy « - Vip^Ss • (30(0)

A{rMb)c
Ie0 + ke0 = VipVtSc \ae (^ + 6j + M^i
1
Ie6 + ke0 = /2pV2Scaed (30(ii))
(33(H))
noting (—Mg) is positive.
At the critical flutter condition, the two frequencies of At the critical flutter condition the effective dampings in
equation (30) must be identical, so each of equations (33(i)), (33(H)) must be zero. The damping in
the bending equation (33{i)), will be eliminated if (ys/3V + 0/2)
ky + VzpVF2Ss is in phase with y, ie
4(-Me VzpVj-Scae)
Oif
-{ (34)
. . . (31)

where Ky is a constant. And the damping in the torsion

It is seen from equation (31) that the flutter frequency is greater equation, (33(H)), will be eliminated if (ae ys/2V + Me6c/V) is
than the bending frequency <oy( = ky/Iyy/2) but less than the in phase with 6, ie
torsional frequency <ue( = (kg/Ie)'/2).
From equation (31), the flutter speed is given by a e ds 0c
j . \4 - v a . . . (35)

V2pVF2Ssa =
(-Me)k(
('-3) (32) where K2 is also a constant.
(lsle +(-Mi)c\
Now writing
\4 c Iy as )
0 = d0 sin cot
This example describes a practical type of wing flutter, which is
y = y„ sin (cot + <f>) • (36)
primarily in torsion, although the flexure-torsion coupling is
essential. This type of flutter cannot be estimated on the basis
of quasi-aerodynamics since the term Me is crucial; the value of where 0O, y0 are the amplitudes of motion and (f> is the phase
Mg to be consistent with the flutter frequency, (ie v = <i>pC/VF). angle, it follows from substitution of equations (36) in
equations (34), (35) that
A stability analysis shows that the system is stable for V < VF
and unstable for V > VF, as required.
Equation (32) highlights other main features of practical sin Zm = — (37)
flutter (remembering Iye = 0): ae 3
0o _ v 2s
(i) flutter speed increases as ec (ie the distance between the (38)
centre of pressure and flexural axis) decreases, y0 sin$ 3c

(ii) flutter speed decreases as axy —* o>s; it is therefore

important to keep the still air natural frequencies well Kx = v cot <f> (39)
separated, 3c
(iii) in practice (coy/t^)2 is small compared with unity, so the
flutter speed increases with ke; thus ke is the primary
structural design variable for increasing the flutter speed. K2 = — 7 ae sin2<t> (40)

It is noted that as (— Me) tends to zero the flutter speed tends to where v is the frequency parameter at the critical flutter
zero. This result is consistent with the statements of Example 2. condition,
It is an interesting observation to note from this example that
if additional structural damping were added to the torsional
mode, it would complement (—Mg) and therefore increase the _ oiFVF
flutter speed. But additional structural damping in the bending
mode would have little effect, in fact it would decrease the
flutter speed slightly. and (— Me(v)) is given from Fig. 4.

Aeronautical Journal October 1985 Hancock etal 291

 Flutter is possible when equations (37-40) are satisfied and (iii) For the flutter frequency to be between the natural wing
when the bending and torsional frequencies in equations frequencies of bending and torsion, then from equation
(33(i)), (33(H)) are the same, that is, by including equations (41)
(34), (35),
Kx > 0 and ae + K2 > 0.
_ ky + VipV^SsaK-Ay)
u>f From equations (39), Kx is positive with either

0 < 0 < 77/2 Or 77 < <f> < 37T/2.

2
= ke-V2PVF Sc(ae + K2(v))
• • (41) Both of these ranges for <j> are also compatible with the
equation of the phase angle <\>, namely equation (37).
The condition for K2 is automatically satisfied from
The complete solution can be obtained from a simple iterative equation (40).
process between equations (37)-(41).
The above equations provide further insight into the flutter
phenomenon for a mass balanced system. Example 5

According to equation (37) the phase angle depends on All of the previous examples have assumed the system to be
(i)
the flutter frequency parameter v, M-e{v) and e. It is of mass balanced. So finally, consider an interesting example
interest to note, as shown in Fig. 6, that in the range of where the wing is not mass balanced but the flexural axis is
interest of v the combination (— M-e{v))v for a flexural assumed to pass through the aerodynamic centre (so e = 0),
axis at the mid-chord, varies little, with a value about 0-4. and the aerodynamic term {M-JdcIV) is neglected.
Thus, from equation (37), The basic equations (18) become

s i n 2<j) •
018
e
(42) Iyy + Iyee + kyy = = -%,**.(£ + |). • (44(i))

Hence with e 0-25 the phase angle between bending iyey + iee + kee = o (44(H))
and torsion is
Consider the torsion equation (44(ii)); for there to be no
0-23° (43) effective damping in the torsion, y must be in-phase with 6. But
turning to the bending equation if 0 is in-phase with y then it is
It is of interest to note that according to equation (42) not possible to eliminate the y damping term, unless
flutter cannot occur in this example if e is less than 0-18.
But this conclusion should be treated with caution since •y = 0. . . . (45)
M-e also depends on e because that value M-e depends on
the position of the flexural axis {see Appendix 2). And then equations (44(i)), (44(H)), become

Iye6 + V2pV2Ssa/20 =0 . . . (46(i))

(ii) The amplitude ratio from equation (38), with the phase
angle from equation (43) is Iee + kee = 0 . . . (46(H))

Flutter can only occur if both frequencies are the same, ie

He. 6v
7o V2PVF2Ssa/2
O>F =T (47)
he
taking s/c ~ 3. Thus for typical value of v between 0-5
and 1 -0 the torsion motion is the dominant one.
Thus there is a finite flutter speed when Iy0 is positive (ie when
the centre of mass lies aft of flexural axis).
In this example the flutter frequency is the uncoupled
torsional natural frequency (remember that the still-air
frequencies are coupled). Furthermore, there is no bending
motion but the bending degree of freedom is vital. This type of
flutter depends on the total bending moment due to 0 (inertial
and aerodynamic) being zero. It can be shown that the
addition of a small (—Mg) term does not significantly change
(-M,M)v
the above result for reasonable magnitudes of Iy6.
0-4-
Once again it is seen that ke is the primary structural variable
0-3-
which determines the critical flutter speed.
0.2-

0'1-

1
6. ENERGY CONSIDERATIONS
0-5 I'd
1-5
y .

A common approach to the understanding of the physics of

flutter is to consider the net flow of energy into the structure
Figure 6. Variation of (-M^(v))v. from the airstream. If there is a net flow of energy into a

292 Hancock eta/ Aeronautical Journal October 1985

 structural vibration, the structure will be unstable and flutter parameter is important and the value of the frequency para-
will occur; if there is a net flow of energy out of a structural meter depends on the inertial and stiffness characteristics.
vibration, the structure will be stable and any oscillations will In equation (51)
decay. At the critical flutter condition there will be a zero flow
of energy.
Now from equations (18(i)), (18(ii)) the bending and ( r A = * ) '•ViirpVtSsa
torsional moments are

B = - Iyy - Iye6 - VipW-Ssa (ys/3V + 0/2) represents the energy associated with the aerodynamic torsion
moment driving the bending motion, it is seen that this energy
is a maximum when <t> = IT12;
T=-Lty-I$ + VipV1Sc
(„[fv+.]-<-<) (yo0oev cos $/2 ] Vb np^Ssa
• (48)
represents the energy associated with the aerodynamic bending
where (— Mg) is positive. The work done W on the structure moment driving the torsion motion, it is seen that this energy is
per cycle is a maximum at <j> = 0;

fy02ps/3c\ x
/V2irpVzSsa
W = J (Bdy + Tdff)
represents the energy required to overcome the aerodynamic
cycle
bending damping;

fa\- Mi,(v))p/a) 1/2TTpV2Ssa

4/ (By + T0)d(a)t). • • (49)
represents the energy required to overcome the aerodynamic
torsional damping.
For flutter instability, W must be positive, so the first
Writing, term, which represents the driving mechanisms, must over-
come the two damping mechanisms. Energy extraction for an
0 = 0o sin bit, y = T0 (sin att + <f>) (50) airstream cannot occur in y, or 0 motions perse, it only occurs
when both are present with an appropriate phase angle. As
at the critical flutter condition, where 0o, y0 are the amplitudes far as helping towards the understanding of flutter, equation
of motion and <f> is a phase angle, then the work done is (51) is uninformative, it does not explain in any way how modes
couple together to give an oscillatory motion at a specific
frequency and at a specific flight condition.
* = fcirplVS,. faffi + J cos*) - y.*% It is left to the reader to check that the flutter condition in
each of the five examples in Section 4 when substituted into
equation (51) satisfy the condition W equal to zero.
-e0H-Mb{v)vla)\ =Q (51)

at the critical flutter speed where

7. THE FLUTTER ENGINE
v = (oFVF/c.
It is often said that the famous Duncan 'flutter engine'
It is seen that there are no direct terms involving the struc- demonstrates wing flutter, this is not necessarily so. A photo-
tural inertias or stiffness (all of these terms are conservative and graph of the flutter engine is shown in Fig. 7.
so do not contribute to the work done over the cycle). Thus it is The flutter engine is a device by which a mechanism with two
said that the work done only depends on the aerodynamic degrees of freedom of y and 6 (as previously defined) extracts
terms, but this statement is superficial since the frequency energy from a uniform airstream and then via gears and
linkages, drives a shaft which may be connected to a generator.
In the flutter engine, there are no counterparts of the wing
structural stiffnesses ky and ke; these are effectively zero.
Furthermore, the amplitude ratio and the phase angle between
y and 0 use adjusting screws in the mechanism in order to
optimise useful work output at any airspeed V. The details of
the mechanism need not concern us here.
The formula for the work done on the flutter engine when it
is oscillating in simple harmonic motion is exactly the same
formula as that derived for the wing (see equations (47)-(50)
namely,

fy = VnrpVtSsa \yoeo ^ ? + y060 e v cos 0/2 - y02w/3c

-
0«o2 ( " Mi»("))"/flJ (52)
Figure 7. Flutter engine.

Aeronautical Journal October 1985 Hancocketal 293

 In practice, flutter engines are designed with e close to zero:
hence the second term in the brackets of equation (52) can be 8. GRAPHICAL SOLUTION OF BINARY FLUTTER EQUATIONS
ignored.
By adjustment of the phase angle 4>, and the angular The physics of flutter, as indicated by the examples in Section
amplitudes y0 and d0, a useful work output may be obtained. 5, is not an easy phenomenon to explain in qualitative terms.
On the basis of equation (52) putting e equal to zero the best From a design point of view what is required is a feel for how
choice of phase is the flutter speed depends on the main variables; position of
0 = IT12 centre of mass, position of flexural axis, the natural torsional
leading to the average power output and bending frequencies, and torsional aerodynamic damping.
A graphical method to illustrate the influence of these
cu n, _ pai^Ssa fyo0o Jo2™
P = eJi-Miiv)) v/a\ variables on flutter was first introduced in Ref. 2 and further
2TT 3c developed by Niblett (3) . This method is described here for the
coupled equations derived in Section 4. Those equations
(equations (18(i)) and (18(ii)) are

. . . (53)
Iyy+IyeO + kyy=-Q^+^
It is noted that if either y0 or 0o is zero, there is no possibility of
positive work. The type of motion (with <f> = IT12) is precisely (57)
that shown in Fig. 5, but as explained earlier, such a motion is
not indicative of wing torsion-flexure flutter. Reference to the
flutter engine as an aid for understanding practical torsion-
flexure flutter phenomena on aircraft wings is therefore mis-
leading. where
It is nevertheless of interest to pursue the treatment of the Q = Vip^Sca.
flutter engine a little further. The output shaft is driven usually
by a slider-crank-chain output from the flapping motion, y, at a At the critical flutter condition
1:1 gear ratio. To obtain a reasonable output, 0o has to be
sizeable since any output depends principally upon the product y = y-e><»< a n ( j 0 = 0 ^ (58)
0oyo. It follows that if stall is not to occur, the throw of the crank
in the slider crank chain cannot be too large. If the power
demanded by the load is PL (including any friction within the
engine itself), then P is equal to PL and equation (53) yields the On substitution of equation (58), equations (57) become
running speeds o>(rad/s), in terms of the frequency parameter v
( = (oc/V), via the v roots of the quadratic equation,

[ay02(s/c) + 30 O 2 (- M-0)V -
2 2
(54)

. . . (59)
where
ScPL
P- 'power/wind speed ratio'. . . (55)
pV^Ssa
For a non-trivial solution of equation (59) the determinant must
be zero. Since the determinant is complex the real part and
the imaginary part of that determinant are separately zero.
The engine will stall if After some manipulation

p, ^ 3/8 a (yoeofl[ ay02(s/c) + 30O2 ( - Mb)]. real part:

The no-load running speed, assuming zero friction, is given

by /x = 0 in equation (54), ie ['-#]-[-Hi-s Iy 2c k9

pSc s2 M-e
v = 3al& [ay02(s/c) + 3d2 ( - Mi)] (56)
6L ^M'-Sh • • • «•»
The fact that there are two running speeds (given by the
solutions of equation (54) is a consequence of the specification where o>#, u>y are the natural (still air) frequencies in torsion and
of a power output, PL. This may be achieved by a large torque bending respectively, (note when
at low o> (the smaller v solution of equation (54)) or a small
torque at large o> (the larger v solution). If the engine is to be
used to generate electrical power, o> will be specified along with e> 0 0 2 5 ,
pS
^f Me
' <<e and negligible).
PL; V and yo/0o must then be adjusted to give this condition. 6/,

294 Hancock etal Aeronautical Journal October 1985

 imaginary part: In Fig. 8(ix) when the wing is mass balanced about the
mid-chord the straight line for a given value of (— Me)
intersects the 'conic' twice, once at the flutter onset con-
dition, and later when (1/2pV25ca/fc9) is about 150, which
(61) represents the speed at which flutter ceases (above this speed
the wing does not flutter). This cessation condition, which
occurs when the wing oscillating in its natural bending mode, is
the condition identified in Example 2. The fact that in Fig. 8(ix)
all the straight dotted lines for different values of (— M-e) go
Now the 'real part' equation (equation (60)) is a quadratic through the same flutter cessation point is a consequence of the
relation between co2 and Q(= VipV2), while the 'imaginary part' neglect of the small relative density term in equation (60).
equation (equation (61)) is a linear relation between a>2 and When IyB is negative (ie centre of mass ahead of flexural axis)
Q(=V2pV2). Plotting the two graphs of to2 against V2QV2 for as shown in Figs. 8(ii, iii, vi), the wing has aflutterinstability for
each of these relationships, the flutter condition (caF, VF) is small positive values of ( - Me) at all speeds. At larger values of
given by the intersection of the two graphs. (— Me) the straight dotted lines cut the conies at two points, the
Defining lower value of {VjpV^Scalkg) is the flutter onset condition
positioning of centre of mass = xcmc aft of leading edge while the higher value is the flutter cessation condition. At high
position of flexural axis = xf aft of leading edge enough values of (—Me) there is no flutter condition since there
aerodynamic centre = 0-25 aft of leading edge is no intersection. In all of the cases shown, for values of (-Mj)
m = mass of wing of 0-5 and 1-0, wing divergence occurs before flutter, this
conclusion is not valid for smaller values of (— M0).
and taking When Iy6 is positive (ie centre of mass aft offlexuralaxis), as
_ ms2 shown in Figs. 8(iv, vii, viii), flutter occurs at non-zero speeds
y for all values of (— Me), and the flutter onset speed increases
T with the magnitude of (— Me). Keeping the flexural axis fixed,
the movement of the centre of mass aft reduces the flutter
Ie (about axis through centre of mass) = mc2xcm2fl speed. In all of these cases flutter occurs before wing diverg-
Ie (about flexural axis) = mc2xcm2l3+ m (xm — x^c2 ence.
Figure 8 gives a bird's-eye view of the effects of the influence
Iy» = m (JOT - xf) 0-45 sc of the various design parameters or flutter. The lecturer is
encouraged to write his own mini-program and to look a little
K =(47T)2/r deeper into the interaction between the variables in the sense
of relating ke, Xf, Iyg etc to real stiffness characteristics of an
k
« =(20T7-)2OTC2/12 aerofoil, looking at spar and skin thickness dimensions and to
real inertial characteristics, including fuel in internal tanks.
equations (60, 61) have been solved for

= 0-3,0-4,05,

*f = 0-3, 0-4,0-5,
9. SUBCRITICAL RESPONSE
-Mj, = 00, 0-5,10,
What has been described so far is the system behaviour at the
and the results are shown in the nine diagrams which make up critical flutter condition itself. It is also of fundamental
Fig. 8. The plots are of importance to know about the response characteristics of the
system at flight speeds below the criticalflutterspeed (so-called
subcritical response) because flight flutter tests, which are
/a>V
againstQ Vip^Sca (62) mandatory, cannot approach the critical flutter condition too
(s) *r ^ — closely. Comparisons have to be made between the measured
and predicted subcritical response in flight for V < VF.
Results have been calculated for the two degree of freedom
model used so far as described by the equations.
The solid lines, given by equation (60), are conies, the straight
dotted lines are equations (61). It should be noted from Iyy + Ij + kyy=- VzpVSsc^ + y)
equation (60) that the wing divergence condition is given when
one of the conic lines intersects the abscissa axis (ie when co is
zero).
lyeJ + h ,6 + kee = VipWSc U f e + o\ + MJ£ .
In the diagonal diagrams, Figs. 8(i), (v), (ix), the wing is mass . . (63)
balanced about the flexural axis, as the flexural axis is moved
aft. Here the conies degenerate into two straight lines. It is seen In this worked example,
that in all three of these diagrams for (— Me) zero the flutter
speed is the zero speed condition, as stated in Section 5 in the (i) the uncoupled natural frequency in bending is taken to
discussion of Example 2. In these cases of mass balance the be about 2Hz
flutter speed increases with increase in (— M-e) the aerodynamic (ii) the uncoupled natural frequency in torsion is taken to
torsional damping. Furthermore the flutter speed decreases as be about 8Hz
the flexural axis moves aft. In these cases flutter always occurs (iii) the centre of mass is taken to be 0- \c aft of the flexural
before wing divergence. axis, at 0-5 s from the wing root

Aeronautical Journal October 1985 Hancock etal 295

 300-

""--^ Xp -0-4
TF)

300 !„-»3
K, .0-3 200- ^v \
*.N
Hr----..^ i„.o 200-
*.f ^v N
- N X ( -0'5

(00
\^ • FLUTTER ONSET \ N
S
• FLUTTER ONSET
\ N • FLUTTER CEASES
100-
100-
too —-_^_ (-Mil-O-5 N.
~~ - - J,-M*0* \ . ^v
f.*." ~^V
- tMJrt° ~* ~ *~ ^ * ^ ~ - ~ -
tn-i» ^ ^
^ - W I N G DIVERGENCE
a 2 i, 6 S 10 12 . U 0 2 A 6 B 10 12 U 16

)ipv'Soo / k t Sf'S'H.

Figure 8(i) Figure 8(ii). Figure 8(iii).

Kern""-'
X, - 0 3
I„>0
--_J-Mil-0 -J-M.1-0

FLUTTER ONSET
• FLUTTER ONSET
FLUTTER ONSET
• FLUTTER CEA5ES

WING DIVERGENCE

^pVsc V k . HPVscH# V**«n,,

Figure 8(iv). Figure 8(v). Figure 8(vi).

WO-5
Xt -0-3
I|.»0

(-Mj)-O

• FLUTTER FLUTTER ONSET

ONSET FLUTTER CEASES
(-MJ-0-5
0-5
Me)"1-0

WING DIVERGENCE WING DIVERGENCE

10 12 U 16
>2PV*Sca
%.
Figure 8(vii). Figure 8(ix).

Figure 8. Graphical solutions of binary flutter.

296 Hancock eta/ Aeronautical Journal October 1985

 (iv) e = 0-2
(v) s/c = 3-0 10. FLUTTER CALCULATIONS
(vi) M-e = -0-30
(vii) relative mass parameter ml1hpSwc = 200 The purpose of a flutter calculation for a coupled wing flexure-
(viii) the results are plotted as a function of 1/2pVzSca/k( torsion system is to determine the critical condition and hence
to determine the flutter speed, the flutter frequency and, if
required, the modulus and phase of the modal relationship
between bending and torsion.
At any speed V, because equations (63) combine two second Since only the critical condition of oscillatory motion is
order systems, there will be four characteristic roots (or sought, oscillatory aerodynamics can be used. As the name
eigenvalues) in the form of two complex roots, when V is well implies if a wing is oscillating in simple harmonic motion with
below the wing divergence speed VD. Normally there will be
two (complex) modes with an amplitude and phase relationship y = yeiM and 0 = 0?°" . . . (64)
between y and 0 for each of the two modes.
At any speed there will be Mode 1 and Mode 2. In each where y/0 is complex, then the aerodynamic bending and
of these Modes the torsion and bending degrees of freedom torsion moments can be expressed as
are coupled together at a particular amplitude \y/0\ and phase
angle arg (y/0), moving together in a motion of the form aerodynamic bending moment
e'M+"°'',withthe same rate of damping p. and same frequency &>.
For both Mode 1 and Mode 2 variables (\y/0\, arg y/0, p., co) will V2pV2Ss[(By(v) + ivBy{v))ye^ + {B0{v) + ivBi,(i>)jdeu»']
change with V. As V —* VF then for one of the Modes p. —* 0,
w -> o)F, and \y/0\ and arg(y/0) tend to be the flutter mode. The aerodynamic torsional moment aboutflexuralaxis
other Mode still exists at the flutter frequency but will be
heavily damped. = VipV^Sc [(My(v) + ivMyiv^ye^ + (Me(v)
Results for the present example are shown in Figs. 9.
Figure 9(i) shows the frequencies of Mode 1 and Mode 2 as a + ivMl)(v))eeiu'] . . . (65)
function of dynamic pressure. At V equal to zero (ie zero
airspeed) the system oscillates in its natural modes; Mode 1 is where
identified with the higher natural frequency of about 10Hz
(which is primarily torsion with a little bending, because of the v = frequency parameter = coc/V
cross inertia), whereas Mode 2 is identified with the lower
natural frequency about 2Hz (which is primarily bending with a (in Europe and the USA, the reduced frequency k(= <DC/2V) is
little torsion). It is seen that the frequency of Mode 2 decreases used instead of v (so v = 2k)), and
markedly with increase in dynamic pressure whereas the
frequency of Mode 1 increases only slightly. By(v),Be(v),My(v),Me{v)
Figure 9(ii) shows the variation in the damping of the two
Modes with dynamic pressure. At zero air velocity, there is no are known as the in-phase derivatives,
damping since any structural damping is assumed to be zero,
and the aerodynamic loads are zero. At low dynamic pressure, By{v), Bb(v), M-(v), Mfc)
the dampings in both Modes increase with increase in dynamic
pressure. It is seen that Mode 2 goes unstable at (1/2pV2Sca/fc„) are known as the quadrature or out-of-phase derivatives.
equal to 1-4, so this is the critical flutter condition. Meanwhile
the damping in Mode 1 continues to increase, becoming large. In a dynamic motion, these oscillatory derivatives should not
The rate of loss of damping in the neighbourhood of flutter is be referred to as stiffness or damping terms because their effect
severe; this type of flutter is known as a 'hard' flutter, typical of depends on the phase relationship between y and 0.
wing flutter, (a 'soft' type of flutter is where the rate of loss of These oscillatory aerodynamic derivatives should not be
damping through flutter is small, as might occur with store confused with aircraft stability derivatives which are defined in
flutter). a totally different manner; an illustration was given in the
Figure 9(iii) shows the variation of the amplitude ratio of discussion of M-e(v) in Section 3.2.
}y/d\ for Mode 1 with variation of dynamic pressure. At zero Two dimensional oscillatory derivatives are tabulated as
airspeed \y/0\ is about 7-0, confirming that in still air Mode 1 is functions of v and free steam Mach number and can be used in
the 'bending' mode. As the dynamic pressure increases the strip theory analyses to derive the above wing derivatives. The
amplitude ratio \y/0\ decreases, indicating that the pitch values can be modified to take into account experimental
amplitude becomes more significant. measurements. Some typical graphs are shown in Appendix 2.
Figure 9(iv) shows the variation of the amplitude ratio \y/0\ Full 3-dimensional derivatives can be calculated from oscilla-
for Mode 2 with dynamic pressure. At zero airspeed \y/0\ is tory finite wing theory.
about 1/6, confirming that the mode is primarily torsion. As the In incompressible fluid dynamics certain terms which are
dynamic pressure increases the bending contribution actually proportional to IP- can be identified in the oscillatory deriva-
decreases, the motion becoming more or less pure torsion even tives; the so-called apparent mass terms. These i*2 terms
up to the flutter condition where \y/0\ is about 1/15. This effectively give the asymptotic values of the fluid dynamic
behaviour tallies with that of the examples discussed earlier. forces as v —* °° (ie either as o> —» °° or as V—> 0, since v = wc/V).
Figure 9(v) shows the phase angle of Mode 1. At zero Since the fluid dynamic forces are proportional to VipV2 the
airspeed y and 0 are in phase; the phase angle then increases product of VipV2 and the v2 terms gives a finite limit as V—> 0
with dynamic pressure. and so the apparent mass terms are the fluid dynamic forces
Figure 9(vi) shows the phase angle of Mode 2. At zero acting on an oscillatory profile in a stream of zero velocity.
airspeed y and 6 are in antiphase (ie phase angle of 180°). But Apparent mass is associated with the fact that in an incom-
with increase in dynamic pressure there is a dramatic change in pressible fluid, there is an instantaneous effect throughout the
the phase angle, there is a rapid decrease to zero, becoming entire flow field due to a motion of a profile because the speed
slightly negative (about -10°) at the flutter condition. of sound is infinite.

Aeronauticel Journal October 1985 Hancock etal 297

 70-
Q ^ ^ FREQUENCY
60-

50-

40-

Mode 2
30-

20-
•Model

10-
IpVSca/
n1 l 1 T
0'5 VO V5

2 /k
(i) e

AMPLITUDE RATIO

Figure 9. Subcritical response.

But in aerodynamics, because air is compressible, distur- where vF = (oFVFlc. These two complex equations (66) are to be
bances take time to radiate through the flow field, conse- solved for coF, Vp, \y/d\ and arg(y/<?).
quently, there are no v1 terms in the derivatives, and apparent
mass terms cannot be identified. There has been considerable An iterative procedure can be followed. An initial value of vF
confusion about apparent mass, only because some people can be assumed so that values can be given to the aerodynamic
needlessly wish to retain the concept. Tabulated derivatives derivatives. The simultaneous equations (66) can then be
give the totality of the aerodynamics loadings; there is no need solved, following the method outlined in Section 8, for <oF and
to bother about or attempt to identify the so-called 'apparent VF, which gives a new value of vF. And then the process can be
mass terms'; whatever they are, they are automatically repeated until convergence of vF.
included. Any attempt to identify apparent mass terms in, for An alternative, indirect, method can be used to solve
example, wind tunnel testing is completely artificial. equations (66). Initial values of <oF and V^ (and hence vF) are
The full equations, equations (18), (65) for simple harmonic assumed and y is taken to be I/O. Then 0 can be found from
nxotion at the critical flutter condition are equation (66(i)). Substitution of 0 into equation (66(ii)) gives a
(complex) error between the left hand side of the equation and
'.'v-V^y - o&yfi + Ky = Vipv/ss [{By(vF) the right hand side. The method updates the values of <nF and
VF in a Newton-Raphson manner to reduce the error to zero to
+ ivFBy(pF)}y + {Be(vF) + iv^iv^O] . . . (66(i)) give the correct critical values of wF, VF and y/0. A code for a
binary flutter calculation is given in Ref. 4.
- <oF2Iyey - vF2le6 + ke0 = V2PVF2Sc[{My(vF) +wFMy^F^j
The methods described here can be extended directly to
+ {Me(vF) + ivF Me(vF)}e] . . . (66(H)) calculations involving more realistic wing modes'5'.

298 Hancock etal Aeronautical Journal October 1985

 TABLE 1

Mode No. 1 2 3 4 5 6 7 8 9 10
Mode Type B1 B2 T1 B3 0A IA T2 B4 T3 B5
Hz 207 6-48 10-36 13-85 17-43 18-75 21-61 24-64 29-71 35-14

B = Mode dominated by wing bending

T = Mode dominated by wing torsion
OA = Mode dominated by outbound aileron deflexion
IA = Mode dominated by inboard aileron deflexion

by increasing the number of modes; taking only modes 1 and 3

11. SOME PRACTICAL FLUTTER CHARACTERISTICS gives an overestimate of the flutter speed by 1-25% while
taking only modes 1, 2, 3 gives an underestimate of the flutter
In practice for the design of aircraft it is necessary to consider speed by 3%.
more than two degrees of freedom to obtain an accurate In practice theflexiblefuselage connection between the two
estimate of the critical flutter condition. For example, as half wings needs to be considered, so both symmetric and
indicated below, a typical transport wing flutter calculation asymmetric wing flutter needs to be calculated.
could involve up to ten normal wing modes. Physical insight into extensivefluttercalculations is well nigh
The derivation of the basic equations is via the principle of impossible, so that if insufficient flutter margin is predicted
virtual work either directly 6T indirectly through Lagrange during the design phase, it is by no means obvious how to
equations. These two formulative procedures are essentially redesign the wing structure to increase theflutterspeed. One of
the same, the advantage of working from first principles rather the most imaginative and major contributions to contemporary
than applying the Lagrange equations is that the methodologies flutter is due to Baldock(6); he developed a technique of
are then seen to be compatible with the Rayleigh-Ritz method systematically reducing a system comprising a large number of
and the finite element displacement method in structural normal modes to an effective binary system. This binary system
analysis. The approach is summarised in Appendix 3. can be understood in qualitative and graphical terms similar to
the graphs in Fig. 8) and appropriate structural changes can
A typical set of data for the normal modes for a wing of a then be identified.
transport aircraft (wing aspect ratio 8, leading edge sweep 30°)
is indicated in Table 1.

CONCLUDING REMARKS

TABLE2 The aim of this note has been to outline some of the underlying
ideas on wing flutter which could form the basis of an intro-
Flutter Speed Flutter Frequency ductory course of instruction. Because of limitations in space
Modes Retained m/s Hz (and time), other forms of flutter, which are important in
practice, have not been mentioned, but the underlying ideas
1.3 802-6 6-67
are similar to those presented here. For example, control
2.3 1500
1.2.3 769-2 6-30
surface flutter involves coupling of control surface rotation
1.2.3.4 777-5 6.33 with wing degrees of freedom, and store flutter involves
1.2.3.4.5 794-8 6-63 coupling of store, pylon and wing degrees of freedom. It is
1.2.3.4.5.6 803-5 6-66 possible for a wing to flutter in the one single degree of pitch if
1.2.3.4.5.6.7 790-3 6-71 the wing goes in and out of stall during the oscillatory pitch
1.2.3.4.5.6.7.8 791-0 6-69 motion; this phenomenon is known as stall flutter. Further-
1.2.3.4.5.6.7.8.9 7920 6-69 more, at supersonic speeds, different types of behaviour can
1.2.3.4.5.6.7.8.9.10 792-2 6-68 occur, for example, at low supersonic Mach numbers the aero-
dynamic pitch damping derivative (— M-e) can be negative so a
single degree of freedom motion in pitch can be unstable.
Finally, it should be mentioned that the transonic region is one
Predicted flutter speeds and flutter frequencies for various where flutter is likely but the aerodynamic loadings are most
combinations of normal modes are indicated in Table 2. complicated, involving non-linear aerodynamics, which means
It can be deduced from Table 2 that the two important that the approach used in this note is no longer valid, other
degrees of freedom are the fundamental bending mode (Bl) techniques are needed.
and fundamental torsion mode (Tl), with the first bending Only flutter has been discussed in this note. The important
harmonic mode (B.2) playing an important secondary role. complementary branch of static aeroelasticity has not been
The need to include a significant number of modes in the flutter discussed. In fact static aeroelasticity can only be understood in
calculation to obtain an accurate estimate of theflutterspeed is the context of aircraft stability and control as will be explained
seen. There is no clear trend towards a converged flutter speed in a future note.

Aeronautical Journal October 1985 Hancock etal 299

 4. COLLAR, A. R. and SIMPSON, A. Matrices and Engineering
REFERENCES Dynamics: Programs. Ellis-Horwood (to be published).
5. SIMPSON, A. The Solution of Large Flutter Problems on Small
Computers. The Aeronautical Journal, 88, April 1984.
1. BROADBENT, E. G. Elements of Aero-elasticity. Aeronautical 6. BALDOCK, J. C. A. The Identification of the Flutter Mechanism
Engineering, 26 1954. from a Large-Order Flutter Calculation. RAE TR 78017,
2. FRAZER. R. A.. DUNCAN, W. J. and COLLAR, A. R. Elementary MoD(PE), 1978.
Matrices, Cambridge University Press. 1965. 7. ACUM, W. E. A. The Comparison of Theory and Experiment for
3. NIBLETT, T. A Graphical Representation of the Binary Flutter Oscillating Wings. Manual on Aeroelasticity Vol II, AGARD
Equations in Normal Co-ordinates. ARC R&M No 3496, 1968. (Circa 1960).

APPENDIX 1
Thus the coefficient of iv, which is equal to Le, behaves like tnv
at small v. So Le tends to infinity as v —* 0. The concept of
ON THE BEHAVIOUR OF Le AT LOW FREQUENCIES quasi-steady two dimensional aerodynamics does not apply for
Le, or for any 0 derivatives.
Although the above analysis is presented for incompressible
The following discussion is given for the derivative Le rather
flow, the same behaviour occurs for all subsonic Mach numbers
than Mg, for convenience; the arguments apply equally well to
up to the transonic range.
all two dimensional 8 derivatives.

The formula given in text books for the lift on a two dimen-
sional aerofoil oscillating in pitch in simple harmonic motion
about a mid chord axis in an incompressible stream is
APPENDIX 2
Lift = irP(c/2)2V0 + 2TTPV(C/2)C(I>/2)[V0 + cO/4]

. . . (Al) OSCILLATORY AERODYNAMIC DERIVATIVES

Two-dimensional derivatives are defined by reference to Fig.

where v = a>c/V and C(v/2) is the standard (complex) Theodor- A2.1:
sen function. The first term in equation (A.l) is the so-called
apparent mass term while the second term is referred to as the Lift = 1 /2pV 2 c[(4 + JVLJ)(Z>/4 + (Le + ivtiie]^
circulatory lift.
. . . (A2.1)
In text books, it is argued that as v—> 0, since C(iV2)—»• 1 then
Nose-up pitching moment about the quarter chord axis

lift-* W pV*c{£+(» + £ ) } . . . (A2) = MVi = V 2 pV 2 c 2 [(M z + ivM-z)(z)Vi + (Me + wlH^eY"

. . . (A2.2)
which is interpreted as the circulatory lift being given by the
incidence at the 3A chord point (remembering that in this where
example, the axis at V2 chord). This argument is fallacious.
Zy4c is the amplitude of heave (z is measured downward) of
Now equation (A.l) is not strictly correct since it only applies V* chord axis,
to simple harmonic motion so it should be more correctly 6 is the amplitude of pitch angle (pitch angle is positive
written in the form nose up) with axis of rotation about xk chord axis,
Lift = 7T p V>c h - + C(W2) [ l + ' f ] k ^ < and so the derivatives Lz, Me etc, are given relative to the axis of
twist at the V4 chord. Graphs of these derivatives are shown in
. . . (A3) Figs. A2.

The oscillatory derivative Le is the coefficient of iv inside the

bracket of equation (A3), (see for example equations (65)). / / / / / / / / / / /
Now the limiting value of C(v/2) is, at small values of v, of the
form

C(v/2) = 1 + K(iv(nv) + 0{Mnv) . . . (A.4) V Ycr

1 ^^~
where K is a constant. C
A
Substituting equation (A4) in equation (A3) gives

Lift~2TTty2pV2c){l + iv(Ktnv + V2))eoei"' . . . (A5) Figure A2.1. Notation.

300 Hancock eta/ Aeronautical Journal October 1985

 0-5-

0-8
0 & = ^ ^
0-5 N
\\° i!s V

-0-5-

-1-0-
\ ^0-5

-1-5- \M=O

Figure A2.2. Figure A2.3.

10-0- \

t. \
«-0-
\

frO-
\ ^
\\. 0-8
— 0-5
4-0- ^ ^
^"""--^M-O

2-0-

0' ' I I
0-5 1-0 v's V

Figure A2.4. Figure A2.5.

Oscillatory derivatives.

AeronauticalJoumal October 1985 Hancock etal 301

 0-4- /0-8

7o 1T5 V (-Mi)
(-M2j
0-3-
-0-2-

0-2"
-0-4-

0-1- ^0-S
-0-6-

M=0
0 I 1
0-5 1-0 i!s V

Figure A2.6. Figure A2.7.

25-

1-2- (-Me) 1
I
^0-8
r B,
1-0- 20-

0-8-

06- 1-5-

0-4-
^
X^O'8
0-2- \ _9JL—-^S
1-0-
1
0
/ V i°
1.5 V M=0
-0-2-
V^J""" 0-5
0-5-
-0-4- ^^M=0

-0-6-
0 1 1 1
0-5 1-0 1-5 V

Figure A2.8. Figure A2.9.

Oscillatory derivatives.

302 Hancock etal AeronauticalJournal October 1985

 If the flexural axis is located at
APPENDIX 3
x = Xf
DERIVATION OF FLUTTER EQUATIONS
aft of the leading edge the two dimensional oscillatory deriva-
tives can be written

Lift = VipV2^^ + ivL-z)z + OU + ivL-e)W' Assume an aircraft displacement w(x, y, t) downward given by
the superimposition of N modes
(A2.3)

Moment (+ ve nose up) (*,:M) =!>(^>f) «.« . (A3.1)

= V2PV2c2[(Mz + ivM-z)z + (Me + wM^ey*'

(A2.4)
where where fn(x/c, y/s) is an assumed mode shape for the n,h mode
zc is the amplitude of heave (displacement is measured (n = 1, N), q„(t) is the n th generalised co-ordinate, and c is the
downward) of the flexural axis, mean chord.
In simple harmonic motion at frequency a>
6 is the amplitude of pitch angle (pitch angle is positive
qn(t) = qne« (A3.2)
nose-up) about the flexural axis.

In terms of the 'V4 chord' derivatives

Lift = 1 /2pV 2 c[(4 + ivLz) (z - (Xf - y 4 )0) The aerodynamic surface loading distribution (upward) can
be expressed as
ia
+ (Le + ivLe)'e]e '

= VipV2<i(Lz + ivLz)z + ^Vtfn^ne (A3.3)

{(Le- (xf - VA)LZ) + HU - (xf- v*)L$ny
(A2.5) where
Moment = VzpW[(M z + ivM-z) (z - (Xf - V*)0)

+ (lke + ivMeJBW'*' + (Xf - y 4 )c (Lift)

= 1/2PV1c2[{(Mz
H^H
is the oscillatory loading distribution induced by the oscillatory
+ (Xf - V*)Lz) + ip(Mz+(Xf-V4)Lz)}z + {[Me- (xf- V*)MZ displacement in the n th mode

+ (Xf - V.) (Le - (Xf - y 4 ) 4 ) ] + utfti - (Xf- V<)MZ

+ (xf - y4) (Le - (xf - v^LjjW

. . . (A2.6)
Mv-M
Note that {€n(x/c, y/s)} will be complex with in-phase and
The reader could check the variation of the derivatives with the quadrature components. It is implied by equation (A3.3) that
position of the flexural axis, in particular to verify Fig. 4. the aerodynamic loadings are linear and superimposable.
The above theoretical two-dimensional values tie up reason- The dynamic problem can be replaced by an equivalent static
aoiy well with experimental measurements below the transonic problem by applying D'Alembert's principle, where the 'static
region (see Ref. 7). inertial' load (downward) on an element of massrfm situated at
Crude values of finite wing derivatives can be obtained from
application in equations (66) by assuming aerodynamic strip
theory with the variation of z(y)/c and 0(y) specified, and then d2w
the expressions (A2.3,4), are integrated in the spanwise = -dm
direction. For better values it is necessary to turn to finite wing
theory.

^dm(P(l 7) f) (A3.4)

The displacement function cf„(x/c, y/s) as defined in equation

(A3.1) implies a strain tensor field throughout the structure

N
. (A3.5)

Aeronautical Journal October 1985 Hancocketal

 which in turn implies, through the stress-strain relationships a Once the still air normal modes have been determined then
stress tensor field throughout the structure these normal modes could be used in the above analysis as the
shape functions f„(x/c, y/s) themselves, rather than arbitrary
assumed mode shapes. In this case, since the normal modes are
? = > .g-n^neial. . (A3.6) orthogonal, then
' &
Mmn = Km„ = 0,mi= n. (A3.14)
Now apply a virtual change 8qm, which implies a virtual
displacement The dynamic behaviour at a finite speed would then be
determined once the aerodynamic matrix has been evaluated,
^ = cfm(-/-\bqme^. . . (A3.7) using the normal mode shapes.
& The above analysis can be applied as a first approximation to
Thus the incremental or virtual work done by the inertial and the analysis of a transport aircraft type wing; the wing can be
aerodynamic loads is divided into N spanwise strips as shown in Fig. A3.1.

. . . (A3.8) j-l,N

which is equal to the incremental or virtual strain energy

{sjWs-.'^voolfJ 8y m « a - . . (A3.9)

Figure A3.1. Notation.

by the principle of virtual work where

+ (<Txy)n (exy)m + K * ) „ ( ^ +(o"5C)„ fej^} The overall displacement of the wing is represented by the
displacements of each strip defined non-dimensionally as
Thus cftje"*1 = downward displacement offlexuralaxis at stripy,
?jC"" = nose up pitch rotation of strip; about its flexural axis.
XiU. +
I)*«A =
-V2pV2^^OT,?^m = l,A0
. (A3.10)
These displacements are taken as the generalised co-
where ordinates q„ei'"' (n = 1,2N) where

Mmn = mass matrix

-#4^)*" . . . (A3.11)
{qu q2... q2N) = {&!, \\... TtN\ \,...eN).
The generalised co-ordinates, and the displacement
functions fn{x/c, y/s) as introduced in equation (A3.1) are

/£„,„ = stiffness matrix =Jjjq'„emd (vol) . . . (A3.12)

A-mn = aerodynamic matrix =// £n F , H/ m (-, jjds n = l,N

(A3.13)
. (A3.15)

<ln - Qn-N
The stiffness matrix can be determined in simphfied analyses
from EI and GJ spanwise distributions, otherwise from finite,
element displacement analyses.
Note that rotary inertia terms have been neglected.
Equations (A3.10) can be used to find the normal modes. where *
In this case the aerodynamic terms are put equal to zero (ie
V = 0). The resulting eigenvalue problem can be solved for
natural frequencies and eigenvectors which define the
proportion of each arbitrary assumed shape fn(x/c, y/s) in each
normal mode shape. = 0 otherwise.

304 Hancock etal Aeronautical Journal October 1985

 The (2N x2N) mass matrix, as defined by equation (A3.11), when 0m-n(y) is the continuous spanwise displacement
can be evaluated so function for an applied angle of twist 0m_N at y = ym„N (ie on
strip m — N) and 0n-N(y) is the continuous displacement
function for an applied angle twist 0„-N aty= yn-N (ie on strip n
:
M„ c2 (mass of strip n) n = 1, N -N).
And

Mnn
•'HI
=<?!!!
strip (n-N)
dmixlcf n=N+l,2N in>N
Kmn = 0 for \n<N
m<N
m>N
. . (A3.19)

. . . (A3.16)
The aerodynamic loadings are taken to be two dimensional
on each strip so the loadings are given by the two dimensional
oscillatory derivatives as listed in Appendix 2.
( n = l, N
Mm = c2fffdm(x/c)
\m=n+N
strip n In terms of the notation in equations A2.3 and A2.4, the
aerodynamic matrix then becomes
= product of inertia of strip n
Ann = S£LZ + ivL;] n = l,N
n = N+l,2N
M„, = Sn[Me + ivMi\
m =n—N n = N+l,2N

Mmn = 0 otherwise. (n = l,N

Anm = S„[MZ + ivMz]
\m=n+N
The structural characteristics can be expressed in terms of
bending stiffness EI(y) and torsional stiffness GJ(y) in = N+l,2N
distributions, modified to take into account the built-in root = Sn[Le + ivLg]
{m—n—N
effects. In this case the stiffness matrix is
= 0 otherwise,
d2hn(y)&hm(y)rlln=l,N
Kmn-f EI(y)
dy2 dy2 dy\
m=l,N
(A3.17)
S„ = area of strip n.

where hm(y) is the continuous spanwise bending displacement

function along the wing flexural axis for an applied
displacement K, aX y = ym (ie on strip m), and hjy) is the
continuous bending displacement along the wing flexural axis
for an applied bending displacement hn aty = y„ (ie on strip n).
And

in =N+1, IN
Kmn=fGJ(y)d0"-^Ue^dy \m = N+l, IN REPRINTS
dy dy
Reprints of this paper may be obtainedfromthe Society, price
(A3.18) £2.00, including postage and packing.

Aeronautical Journal October 1985 Hancocketal 305