THE AERONAUTICAL JOURNAL 1

Page 1 of 34.
c The Author(s), 2020. Published by Cambridge University Press on behalf of Royal Aeronautical
Society
doi:10.1017/aer.2020.84

A comparative assessment
of flutter prediction techniques
U.P.V. Sudha
[email protected]
Aeronautical Development Agency
Bangalore
India 560017

G.S. Deodhare
Programme Director (Combat Aircraft) & Director
Aeronautical Development Agency
Bangalore
India 560017

K. Venkatraman
Aerospace Engineering
Department Indian Institute of Science
Bangalore
India 560012

ABSTRACT
To establish flutter onset boundaries on the flight envelope, it is required to determine the flut-
ter onset dynamic pressure. Proper selection of a flight flutter prediction technique is vital to
flutter onset speed prediction. Several methods are available in literature, starting with those
based on velocity damping, envelope functions, flutter margin, discrete-time Autoregressive
Moving Average (ARMA) modelling, flutterometer and the Houbolt–Rainey algorithm. Each
approach has its capabilities and limitations. To choose a robust and efficient flutter prediction
technique from among the velocity damping, envelope function, Houbolt–Rainey, flutter mar-
gin and auto-regressive techniques, an example problem is chosen for their evaluation. Hence,
in this paper, a three-degree-of-freedom model representing the aerodynamics, stiffness and
inertia of a typical wing section is used(1) . The aerodynamic, stiffness and inertia properties
in the example problem are kept the same when each of the above techniques is used to pre-
dict the flutter speed of this aeroelastic system. This three-degree-of-freedom model is used
to generate data at speeds before initiation of flutter, during flutter and after occurrence of
flutter. Using these data, the above-mentioned flutter prediction methods are evaluated and
the results are presented.

Keywords: Flight flutter test prediction technique; ARMA; Flutter Margin; Three-degree-
of-freedom model

Received 23 January 2020; revised 26 July 2020; accepted 20 August 2020.

2 THE AERONAUTICAL JOURNAL

NOMENCLATURE
Roman symbols
xα distance from centre of gravity to elastic axis
xβ distance from hinge to elastic axis
rα radius of gyration of the aerofoil
rβ radius of gyration of the aileron
b semi-chord of the aerofoil
a non-dimensional distance between elastic axis and mid-chord of the aerofoil
c non-dimensional distance between the hinge point to the mid-chord of the aerofoil
Kh spring stiffness in heave
Kα spring stiffness in pitch
Kβ spring stiffness of the control surface
h translational displacement
ḣ heaving velocity
ḧ heaving acceleration
nc non-circulatory component
C(k) Theodorsen function
H1 , H0 Hankel function
V Free-stream airspeed
s Non-dimensionalised Laplace variable

Greek symbols
ρ air density
α̇ pitching velocity
α̈ pitching acceleration
β aileron rolling displacement
β̇ aileron rolling velocity
β̈ aileron rolling acceleration
μ mass ratio
ωh heaving frequency
ωθ pitching frequency
ωβ aileron pitching frequency
α rotational displacement

Acronyms
FM flutter margin
AR auto-regressive
ARMA auto-regressive moving average
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 3

MBA moving block approach

LSCFM least-squares curve-fitting method
FPE final prediction error
AIC Akaike information criterion
PSD power spectral density
FFT fast Fourier transformation
Fz flutter prediction parameter for binary flutter system
Fn flutter prediction parameter for multi-mode flutter system

1.0 INTRODUCTION
Flight flutter tests are an indispensable part of the certification process of any new aircraft
programme and are mandatory for envelope expansion. It must be demonstrated that the
aircraft is aeroelastically stable throughout its design envelope. The approach adopted for
flutter clearance involves a combination of pre-flight flutter analysis and flight flutter test-
ing. Flutter prediction techniques lie at the core of all flight flutter tests, and proper use of
these techniques to accurately predict flutter is a challenging task. The commonly used flutter
prediction techniques available in literature are velocity damping, envelope functions, flutter
margin, discrete-time ARMA modelling, flutterometer and the Houbolt–Rainey algorithm. In
this work, using a canonical three-degree-of-freedom aeroelastic system, some of these dif-
ferent flutter prediction algorithms are examined in terms of their accuracy, efficiency and
robustness.

2.0 LITERATURE REVIEW

Flight flutter testing is a mandatory part of an aircraft’s certification process and thus must
be undertaken to demonstrate that the aircraft is flutter free throughout its desired flight
envelope(2) . To establish the flutter onset boundaries on the flight envelope, one has to deter-
mine the flutter onset dynamic pressure. The method that is most commonly used to predict
the onset of flutter is to extrapolate the trends of modal damping(3,4) . This method consists
of noting the variation in the modal damping with airspeed and extrapolating those varia-
tions to an airspeed at which the damping should become zero. The resulting airspeed is
considered as the predicted flutter speed. Desmarais and Bennett(5) developed an automated
procedure for implementing the traditional V –g method of flutter solution. Koenig(6) was
the first to emphasize the problems with using damping estimation as an indicator of flut-
ter onset. Adverse effects such as noise and poor excitation will cause a scattered damping
estimation(6) . Another area of difficulty is the extrapolation of the damping. Damping can be
a highly nonlinear function of airspeed, thus its extrapolation must be performed carefully
to ensure that any higher-order nonlinearity is accounted for. Velocity damping trends are
shown in Figure 1(7) for three different configurations that are prone to mild flutter (on the
left), exhibit moderate flutter (in the centre), and prone to explosive or sudden and violent
flutter (on the right). In all these cases, the damping initially increases until some speed Va
is reached, above which the damping decreases until flutter occurs. From Figure 1, it is seen
that, in the case of mild flutter, the damping degradation above Va is gradual, while in the case
of moderate flutter there is sufficient warning before the onset of flutter. However, in the case
of explosive flutter, the damping degradation above Va is abrupt, with no prior warning. Thus,
4 THE AERONAUTICAL JOURNAL

Figure 1. Damping as an indication of flutter.(7)

the real problem here is to use other flutter prediction techniques to avoid the deficiencies of
the velocity damping technique.
There are several techniques to extract the frequency and damping and determine the stabil-
ity of an aircraft from data obtained during flight flutter testing. Cooper(8) and Kayran(9) have
reviewed various techniques that have been used to estimate modal parameters from flight
flutter test data. These methods can be classified into frequency-domain methods and time-
domain methods, depending upon whether the curve fitting is performed in the frequency
or time domain. Frequency-domain methods—the traditional approach to modal parameter
estimation—use a parameter estimation procedure based on the calculation of the frequency
response function; random noise is removed through averaging. Any errors in the frequency
response function estimation are passed onto the estimated modal parameters.
Time-domain methods, on the other hand, use time response data by curve fitting the
impulse response function. The impulse response function is either calculated from the inverse
Fourier transform of the frequency response function or the time-domain system response to
an impulse excitation, or generated from the system response to an unknown random input.
Time-domain modal analysis offers several advantages over frequency-domain approaches(9) .
It does not depend on excitation or equipment needed to obtain specially designed excita-
tion forces. Also, in general, it needs fewer response data. As the analysis is not based on
frequency response function data, it can represent an alternative for analysing close vibration
modes. However, time-domain modal analysis requires measured data with a very high signal-
to-noise ratio. Data inaccuracy and the presence of noise can lead to computational modes in
the calculation. Cooper(8) does not recommend any particular technique, but refers to the
wide practice of linear modal analysis technique used for estimating frequency and damping.
Dawson and Maxwell(10) used amplitude and velocity inputs, wherein the onset of flutter is
predicted by recording the maximum amplitude data for a structure at various test points.
The maximum amplitude is obtained from a sharp, adequate and precisely repeatable step
excitation or by monitoring the maximum modal amplitude of a structure under the influence
of a frequency sweep excitation. This technique is applied on an F-16 carrying a heavy load.
Cooper(11) and Dimitriadis and Cooper(12) used the envelope function technique, which unlike
modal damping extrapolation, makes use of the envelope bounding an impulse response as
obtained after a sharp impulse excitation of the structure. The size and shape of the response
envelope change, reflecting the loss of damping with varying speed. Flutter instability can
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 5

be predicted by extrapolating the envelope shape parameter, which reaches a certain value as
flutter approaches.
Houbolt and Rainey(13) proposed the Houbolt–Rainey technique, which is essentially a sub-
critical response-based flutter onset prediction that can be used with either a sinusoidally
varying forced excitation or a randomly varying forced excitation. The analytical founda-
tion of this approach is straightforward and simply requires a measurement of the amplitude
of the response in the structural mode that is important to flutter. The reciprocal of these
amplitude measurements are plotted against a flow parameter such as density or dynamic
pressure, as flutter is approached, and the resulting curve is extrapolated to the flutter con-
dition, which occurs when the reciprocal of the amplitude becomes zero. Sandford et al.(14)
presented some results using a subcritical response method that they called the peak-hold
method. Foughner(15) used this method during flutter testing at the NASA Langley tran-
sonic dynamics tunnel. Doggett(16) examined the relationship between the Houbolt–Rainey
and peak-hold method and stated that these two methods are the same, suggesting that the
method be referred to in the future as the Houbolt–Rainey method.
The μ flutterometer method used by Brenner et al.(17) , Lind(18,19,20) and Lind and
Brenner(21, 22,23) provides a methodology for the computation of a worst-case estimate of flut-
ter speed that is conservative with respect to the true flutter speed by considering the errors
and uncertainty in a mathematical model. The flutterometer approach uses flight flutter test
data to identify the errors and uncertainty in a model. The resulting worst-case flutter mar-
gins predicted using the μ flutterometer directly account for data that describe the aircraft, as
distinguished from a mathematical model of the aircraft. However, this method is more con-
servative in predicting the flutter onset speed owing to the incorporation of uncertainty into
the model.
Zimmerman and Weissenburger(7) introduced the flutter margin concept based on Routh’s
stability criteria. The flutter margin (FM) is computed using frequencies and modal damping
values. It was developed based on a binary flutter system. The flutter margin decreases much
more gradually than the damping of the critical mode. The flutter margin is also expressed
as a quadratic function of the dynamic pressure. Hence, if the flutter margin is known at
three different airspeeds, it can be fitted using a second-order polynomial and subsequently
extrapolated to the zero flutter margin where flutter occurs. Price and Lee(24) developed a new
form of flutter margin for trinary flutter that also varies linearly with dynamic pressure; it is
also relatively insensitive to errors in damping, but is very sensitive to errors in frequency.
However, they applied this technique to a specific case of a known trinary flutter problem,
which may not always be the case. Lind(25) introduced an extension to the Zimmerman and
Weissenburger(7) approach that accounts for multi-mode coupling in the flutter instability.
Both the classical two-mode coupling and a multi-mode coupling approach showed similar
accuracy and confidence level when the critical flutter mode was considered in the modal
pairing. Torii and Matsuzaki(26) proposed a new flutter prediction parameter based on the
Jury stability criteria in the sampled time domain, which is analogous to the continuous time
flutter margin defined by Zimmerman and Weissenburger(7) .
Kayran(9) reviewed various flight flutter test techniques utilised in flight testing of aircraft.
The flight test procedures typically followed in flight flutter testing are described for different
airspeed regimes. Kayran(9) mentions that the flutter margin method can be used for systems
having more than two degrees of freedom, if it is known beforehand which two modes will
cause flutter. In that case, the flutter margin can be applied successfully; otherwise, all possible
pairs of modes must be examined. Flutter margin-based flutter prediction was applied to the
F-15 flight flutter test data by Katz et al.(27) . They defined the critical flutter boundary using
6 THE AERONAUTICAL JOURNAL

the flutter margin for all the important modal combinations. When using the flutter margin,
the frequency and damping must be known beforehand.
System identification(28,29) is the process of determining an adequate mathematical model,
usually containing differential or finite-difference equations, with unknown parameters to be
determined indirectly from measured data. These techniques can be used for the estimation of
the frequency, damping and stability from flight flutter test measurements. The stability can
be estimated directly using flight flutter test response data without estimating the frequency
and damping by using such system identification techniques.
Sundaramurthy et al.(30) used a technique based on the flow unsteadiness as the pri-
mary excitation for dynamic stability measurements in a trisonic wind tunnel. Time series
auto-regressive modelling was used to derive a digital spectrum of the unsteadiness excited
model response, and the system damping was evaluated from the half-power bandwidth
of the spectrum. Perangelo and Waisanen(31) examined the capability of the maximum-
likelihood technique to determine the modal frequency and damping estimates from noisy
flutter response signals by comparison with least-squares flutter analysis procedures. Batill
et al.(32) , Pinkelman et al.(33,34) used parameter identification time series models for linear
dynamic structural systems. Total least squares provides an approach that appears to signifi-
cantly reduce the bias error in the parameter estimates. These methods are applied to a simple,
simulated system, and then to flight flutter test data with a particular emphasis on accu-
rate modal damping estimates. McNamara and Friedmann(35) compared three different flutter
identification methods, namely the Moving Block Approach (MBA), Least-Squares Curve-
Fitting method (LSCFM) and a system identification technique using an Auto-regressive
Moving-Average (ARMA) model. Using analytical data, they found the ARMA model-based
prediction tool to be superior to both MBA and LSCFM. Jategaonkar and Balakrishna(36) used
analytical data and a known number of modes to estimate the frequency and damping using
an AR model.
Kay(37) classified spectral estimation techniques into non-parametric and parametric meth-
ods. Non-parametric methods are referred to as classical techniques, and parametric methods
as modern spectral estimation techniques. Fast Fourier Transform (FFT) methods, which are
non-parametric, do not require a model or parameter estimation for the spectral estimation.
The Fourier transformation operations are performed on either windowed data or windowed
Auto-correlation Function (ACF) estimates. Such windowing of data or ACF values makes
an implicit assumption that the unobserved data or ACF values outside the window are zero,
which is normally an unrealistic assumption. A smeared spectral estimate is a consequence
of windowing. AR, MA and ARMA models are parametric methods. Parametric modelling
for spectral estimation consists of choosing an appropriate model, estimating the parameters
of the model and then substituting these estimated values into theoretical expressions for the
power spectral density (PSD). One of the promising aspects of the model-based approach to
spectral estimation is that one can make more realistic assumptions concerning the nature
of the measured process outside the measurement interval, beyond assuming it to be zero or
cyclic. Thus, the need for window functions can be eliminated, together with their impact.
Kay(38) stated that FFT methods are inferior to the AR modelling method, as very long signal
records are necessary for the spectral estimation, highlighting the advantages of modern spec-
tral estimation techniques over classical techniques. Jategaonkar and Balakrishna(36) mention
that, in blow-down tunnels, long-period records are highly uneconomical and FFT methods
are not preferred for analysis.
Raol et al.(39) investigated the problem of model order determination, which is an integral
part of system identification for dynamical systems using auto-regressive and least-squares
techniques. Twelve model order testing criteria available in literature were critically evaluated
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 7

using simulated data. The study revealed that, among all the tested criteria, only a subset are
adequate to establish the model order reliably. Two criteria proposed by Akaike, namely the
Final Prediction Error (FPE) and the Akaike Information Criterion (AIC)(40,41,37) , are found
to be more reliable. These model order estimators are based on the estimated prediction error
power. Pan(42) used the Akaike information criterion in his biomedical studies. Bozdogan(43)
studied the basic underlying idea of the Akaike information criterion and introduced a new
measure of complexity—the information complexity (ICOMP) criterion—as a decision rule
for model selection and evaluation.
Matsuzaki and Ando(44) suggested a system identification analysis procedure using the
ARMA representation of an aeroelastic system as a means of predicting the flutter veloc-
ity. The basis for this approach is the Jury stability criteria(45) for evaluating the stability of a
discrete system, which are very similar to Routh’s stability criterion for continuous-time sys-
tems. This was applied to binary flutter systems. The flutter prediction parameter Fz proposed
by(26) was updated to include multiple modes by(46) , thus defining a new parameter Fn . This
technique has been applied to subsonic wind-tunnel test data. Recently, Torii(47) described a
procedure for three-mode discrete-time systems.
In the discussion that follows, we present an analysis based on simulations using five
different flutter prediction methods, starting with the traditional velocity damping-based
flutter prediction, followed by the envelope function-based technique, the Houbolt–Rainey
approach, the flutter margin-based flutter prediction, and AR-based prediction models. The
same aeroelastic system is simulated and used for flutter prediction by each of the five differ-
ent techniques. The aeroelastic system has three degrees of freedom, and the flutter speed is
computed for this system using the V –g method with unsteady potential flow aerodynamics
resulting from a rational function approximation of Theodorsen’s function.

3.0 AEROELASTIC RESPONSE OF A WING SECTION

This section describes the flutter instability of a wing section with a control surface hinged to
it. The three-degree-of-freedom model of this wing section is shown in Figure 2, where b is the
semi-chord of the aerofoil section, a is the non-dimensional distance between the elastic axis
and the mid-chord, and c is the non-dimensional distance between the hinge point and the mid-
chord, both expressed as a ratio of the semi-chord. We first develop the equations of motion of
this three-degree-of-freedom vibration model, then describe the unsteady aerodynamic terms
and the flutter instability problem.

3.1 Equations of motion

The wing section model with heave, pitch and control surface degrees of freedom having
linear stiffness constraints is shown in Figure 2. The heave stiffness is represented by Kh ,
the torsional stiffness of the spring is represented by Kα and the torsional stiffness of the
control surface is represented by Kβ . h is the vertical displacement of the elastic axis, positive
downward; α is the angle-of-attack, positive clockwise, of the entire aerofoil, about its elastic
axis; β is the angular deflection, positive downward, of the control surface relative to its
undisturbed position. For the three-degree-of-freedom wing section, the differential equations
of motions can be written using Newton’s second law of motion(1) as follows:

mḧ + Sα α̈ + Sβ β̈ + Kh h = F, · · · (1)

Sα ḧ + Iα α̈ + [Iβ + b(c − a)Sβ ]β̈ + Kα α = M · · · (2)

8 THE AERONAUTICAL JOURNAL

Figure 2. The three-degree-of-freedom aerofoil.

and

Sβ ḧ + [Iβ + b(c − a)Sβ ]α̈ + Iβ β̈ + Kβ β = Mβ , · · · (3)

where Sα = mxα , Sβ = mxβ , Iα = mrα2 , Iβ = mrβ2 .

In matrix form, these three equations can be represented as
⎡ ⎤⎧ ⎫
m mxα mxβ ⎨ ḧ ⎪
⎪ ⎬
⎢ ⎥
⎣mxα mrα2
mrα + mxβ (bc − ba)⎦ α̈
2
⎩ ⎪
⎪ ⎭
mxβ mrα2 + mxβ (bc − ba) mrβ2 β̈
⎡ ⎤ ⎧ ⎫ ⎧ ⎫
Kh 0 0 ⎪
⎨h⎪ ⎬ ⎪⎨F⎪ ⎬
⎢ ⎥
+⎣ 0 Kα 0 ⎦ α = M . · · · (4)
⎪
⎩ ⎪ ⎭ ⎪⎩ ⎪ ⎭
0 0 Kβ β Mβ

3.2 Aerodynamic forces

The unsteady aerodynamic forces are calculated based on linearised thin-aerofoil theory.
In this study, Theodorsen’s approach(48) is used to calculate the aerodynamic forces and
moments. Here, we follow the representation of this theory given in(1) .
The total force and moment resulting from the non-circulatory and circulatory flows are
expressed as

V b
F = −πρb ḧ + V α̇ − baα̈ − T4 β̇ − T1 β̈ − 2πρVbQC(k),
2
· · · (5)
π π

1 1 V2
M = πρb2 baḧ − Vb − a α̇ − b2 + a2 α̈ − (T4 + T10 )β
2 8 π

Vb 1 b2
+ −T1 + T8 + (c − a)T4 − T11 β̇ + (T7 + (c − a)T1 )β̈ · · · (6)
π 2 π

1
+ 2πρVb a +
2
QC(k)
2
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 9

b Vb 1
Mβ = πρb 2
T1 ḧ + 2T9 + T1 − a − T4 α̇
π π 2
2
2b2 V
− T13 α̈ − (T5 − T4 T10 )β · · · (7)
π π
2
Vb b
+ T4 T11 β̇ + T3 β̈ − ρVb2 T12 QC(k)
2π 2 π

In the above set of equations, V is the free-stream airspeed, ρ is the air density and C(k)
is Theodorsen’s function as a function of the reduced frequency k, and is associated with the
circulatory component of the unsteady forces. Q in the above is defined as

1 V b
Q = V α + ḣ + α̇b − a + T10 β + T11 β̇, · · · (8)
2 π 2π

and the Ti functions are given in(1) pp.16–20

Q can be rewritten in the following form by using S1 and S2 :

Q = V [S1 ]{xs } + b[S2 ]{x˙s } · · · (9)

where

T10
S1 = 0 1 · · · (10)
π

T11
S2 = 1 (0.5 − a) · · · (11)
2π
⎧ ⎫
⎨h/b⎬
xs = α . · · · (12)
⎩ ⎭
β

In matrix form, the aerodynamic forces and moments will be

⎧ ⎫
⎪
⎨ F/b ⎪ ⎬ 2b
M/b2 = q (2C(k){R}[S1 ]{xs } + C(k){R}[S2 ]{x˙s }
⎪
⎩ ⎪
⎭ V
Mβ /b2

2b2 2b
+ 2 [Mnc ]{x¨s } + [Bnc ]{x˙s } + 2[Knc ]{xs } · · · (13)
V V

where
⎧ ⎫
⎨ −2π ⎬
R = 2π (a + 0.5) · · · (14)
⎩ ⎭
−T12
10 THE AERONAUTICAL JOURNAL

In the above set of equations, the subscript nc denotes the non-circulatory contribution.
⎡ ⎤
−π πa T1
⎢ ⎥
⎢ 1 ⎥
Mnc = ⎢ π a −π + a2 −2T13 ⎥ · · · (15)
⎢ 8 ⎥
⎣ T3 ⎦
T1 −2T13
π
⎡ ⎤
0 −π T4
⎢ ⎥
⎢ 1 ⎥
⎢
Bnc = ⎢ 0 π a − −T 16 ⎥ · · · (16)
2 ⎥
⎣ T19 ⎦
0 −T17 −
π
⎡ ⎤
0 0 0
⎢0 0 −T ⎥
Knc = ⎢⎣
15 ⎥
· · · (17)
T18 ⎦
0 0 −
π
Equation (4) can now be presented as
⎡ ⎤⎧ ⎫ ⎡ 2 ⎤⎧ ⎫
1 xθ xβ ⎨ ḧ/b ⎪
⎪ ⎬ ωh 0 0 ⎪
⎨ h/b ⎪
⎬
⎢ ⎥ ⎢ ⎥
⎣ xθ rθ2
rθ + xβ (c − a) ⎦
2
α̈ +⎣ 0 ωα2 r2θ 0 ⎦ α
⎪
⎩ ⎪
⎭ ⎪
⎩ ⎪
⎭
xβ rθ + xβ (c − a)
2 2
rβ β̈ 0 0 ωβ2 r2β β
⎡ ⎧⎧ ⎫ ⎧ ⎫
q⎢
⎪
⎨⎪ ⎨ −2π ⎪
⎬ ⎨h/b⎪
⎪ ⎬
T10
= ⎣2C(k) 2π (a + 0.5) 0 1 α
m ⎪
⎩⎪ ⎩ ⎪
⎭ π ⎪ ⎩ ⎪ ⎭
−T12 β
⎧ ⎫ ⎧ ⎫⎫
⎪
⎨ −2π ⎪
⎬ ⎪
⎨ḣ/b⎪ ⎬⎪
⎬
2b T11
+ 2π (a + 0.5) 1 (0.5 − a) α̇
V ⎪⎩ ⎪
⎭ 2π ⎪ ⎩ ⎪ ⎭⎪
⎭
−T12 β̇
⎡ ⎤
−π πa T1 ⎧ ⎫
⎢ ⎥ ⎪ḧ/b⎪
2⎢ 1 ⎥ ⎨ ⎬
2b ⎢ ⎥
+ 2 ⎢ π a −π 8 + a −2T13 ⎥ α̈ · · · (18)
2

V ⎢ ⎥⎩⎪ ⎭ ⎪
⎣ T3 ⎦ β̈
T1 −2T13
π
⎡ ⎤
0 −π T4 ⎧ ⎫
⎢ ⎥ ⎪ḣ/b⎪
⎢
2b ⎢0 π a − 1 ⎥ ⎨ ⎬
+ ⎢ −T16 ⎥ ⎥ α̇
V ⎢ 2 ⎥⎪
⎣ ⎩ ⎪ ⎭
T19 ⎦ β̇
0 −T17 −
π
⎡ ⎤ ⎤
0 0 0 ⎧ ⎫
⎢ ⎥⎪ ⎨h/b⎪ ⎬⎥
⎢0 0 −T15 ⎥ ⎥
+ 2⎢ ⎥ α ⎥
⎣ T18 ⎩ ⎭⎦
⎦ ⎪ ⎪
0 0 − β
π
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 11

or more compactly
q 2b
[M]{x¨s } + [K]{xs } = [2C(k){R}[S1 ]{xs } + C(k){R}[S2 ]{x˙s }
m V
2b2 2b
+ 2 [Mnc ]{x¨s } + [Bnc ]{x˙s } + 2[Knc ]{xs }]. · · · (19)
V V

3.3 The rational function approximation of the unsteady aerodynamic

forces
To transform the equations of motion into state-space form, the frequency-domain unsteady
aerodynamic forces must be approximated in terms of rational functions of the Laplace vari-
able. However, when such Rational Function Approximation (RFA) is conducted, it results
in an increase in the number of augmented aerodynamic states required to accurately repre-
sent the unsteady aerodynamic forces. There are several approaches for RFA, such as the
matrix Pad approximant technique(49,50) , minimum-state method(51-53) , corrected/modified
minimum-state method(54,55) , mixed method(56) , Chebyshev polynomials method(57,58) and
Rogers method(59) .
In this paper, the RFA approach based on Roger’s method is used. This method is very
simple and accurately transforms the unsteady aerodynamic forces from the frequency to
time domain(1) . The unsteady aerodynamic forces given by Equation (13) in the frequency
domain can be recast as
⎧ ⎫
⎨ F/b ⎪
⎪ ⎬
M/b2 = q 2C(k){R}[S1 ] + i2kC(k){R}[S2 ] − 2k 2 [Mnc ]
⎪
⎩ ⎪
⎭
Mβ /b2
+ i2k[Bnc ] + 2[Knc ] {xs }
· · · (20)
= q 2C(k){R}[S1 ] + is C(k){R}[S2 ] − 2s 2 [Mnc ]

+ i2s [Bnc ] + 2[Knc ] {xs }
= q[A(s )]{xs }

where s is the non-dimensionalised Laplace variable s = sb/V . Let the calculated aerody-
namic forces be represented as [A(s )] = [F(s )] + i[G(s )]. The real and imaginary part of the
aerodynamic matrix will then be
N
[Pj ]s
A(s ) = [P0 ] + [P1 ]s + [P2 ]s 2 + ;
j=3
s + γj−2

N
[Pj ](−s 2 )
[A(s )] [F(s )] = [P0 ] + [P2 ]s 2 + , · · · (21)
j−3 (−s 2 ) + γj−2
2

N
[Pj ]γj−2 s
[A(s )] [G(s )] = [P1 ]s + .
j−3
(−s 2 ) + γj−2
2

γj−2 are the aerodynamic poles, which are usually preselected in the range of reduced fre-
quencies of interest. The real matrices [Pj ] are determined using the least-squares technique
for a term by term fitting of the aerodynamic matrix [A(s )]. The augmented aerodynamic
state is defined as
12 THE AERONAUTICAL JOURNAL

s s
{xaj } = {xs } = · · · (22)
s + γj−2 s + b γj−2
V

V
{ẋaj } = {ẋs } − γj−2 {xaj } · · · (23)
b
One can now represent Equation (19), together with the unsteady aerodynamic states,
because the augmented state-space equation of motion now represents a linear time-invariant
system.
⎡ ⎤
⎧ ⎫ 0 I 0 ··· 0 ⎧ ⎫
⎪
⎪ ẋ ⎪
s ⎪ ⎢ ⎥⎪⎪ xs ⎪⎪
⎪
⎪ ⎪
⎪ ⎢−M −1 K −M −1 B qM −1 P3 −1
⎥⎪⎪ ẋs ⎪
⎪
⎨ ẍs ⎪
⎪ ⎬ ⎢ ··· qM PN ⎥⎪⎨ ⎪ ⎬
ẋa3 = ⎢ ⎥ xa3
⎢ 0 I − Vb γ1 I ··· 0 ⎥ · · · (24)
⎪
⎪ .. ⎪⎪ ⎢ ⎥⎪ . ⎪
⎪
⎪ . ⎪⎪ ⎢ ⎥⎪⎪ .. ⎪⎪
⎩ ⎪
⎪ ⎭ ⎣ ··· ··· ··· ··· ··· ⎦⎪⎪ ⎪
⎩ ⎪
⎭
ẋaN V xaN
0 I 0 0 − b γN−2 I

where
2
b
M = M − qP2 · · · (25)
V

b
B = B − qP1 · · · (26)
V

K = K − qP0 · · · (27)

4.0 ANALYTICAL FLUTTER PREDICTION

We assume that the aeroelastic system is vibrating harmonically as xs (t) = xs eiωt . Substituting
into Equation 19 and re-arranging terms, we get

1 1 1 1
[K]{xs } = [M]{xs } + C(k){R}[S1 ]{xs } + i C(k){R}[S2 ]{xs }
ω2 μ k2 k

1 1
− [Mnc ]{xs } + i [Bnc ]{xs } + 2 [Knc ]{xs } , · · · (28)
k k

where μ is the mass ratio and k = ωb/V is the reduced frequency. In the V –g(1) or k method(60)
of flutter analysis, Equation 28 is recast as a complex eigenvalue problem by introducing
an artificial structural damping term g. This structural damping force is proportional to the
structural stiffness but in quadrature with the displacement. Note that this structural damping
is frequency independent.

1 + ig 1 1 1
[K]{xs } = [M]{xs } + C(k){R}[S1 ]{xs } + i C(k){R}[S2 ]{xs }
ω 2 μ k 2 k

1 1
− [Mnc ]{xs } + i [Bnc ]{xs } + 2 [Knc ]{xs } . · · · (29)
k k
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 13

Table 1
Typical section parameters for generating the response data(1)

Parameter Values
Heaving frequency, ωh 50.0rad/s
Pitching frequency, ωα 100.0rad/s
Aileron pitching frequency, ωβ 300.0rad/s
Distance to rotation axis, a −0.4
Distance to aileron hinge axis, c 0.6
Semi chord, b 1
Distance from centre of gravity to elastic axis, xα 0.2
Distance from hinge to elastic axis, xβ 0.0125
Square of radius of gyration of the aerofoil, r̄α2 0.25
Square of radius of gyration of the aileron, r̄β2 0.00625
Mass ratio, μ 40.0
Air density, ρ 0.002378slugs/ft3

For a given value of reduced frequency k = ωb V

, the eigenvalue λ λRe + iλim = ω12 + i ωg2
is determined, and we will have λ1Re = ω2 and the structural damping factor g = λλRe im
.
A MATLAB code based on the theory developed above was written to determine the
response data and to obtain the frequency and damping values for different airspeeds. The
parameters for this simulation are presented in Table 1(1) . The different k values were obtained
by varying the aircraft speed V . The response obtained at a speed of 301.68ft/s corresponds
to the flutter speed, and the dynamic pressure at this point is 0.75psi(1) . The acceleration
response data at flutter is shown in Figure 3. Similarly, the responses obtained at a speed of
200ft/s, which is less than the flutter speed, and at a speed of 305ft/s, which is more than the
flutter speed, are shown in Figure 3.
Using this program, it becomes possible to generate the aeroelastic response at different
airspeeds, which is equivalent to obtaining data from flight flutter tests. Therefore, these
response simulations are assumed to provide flight flutter test data and will be used to obtain
the common aeroelastic response data sets for evaluating the flutter prediction techniques
such as the velocity damping, envelope function, Zimmerman–Weissenburger flutter margin,
discrete-time auto-regressive modelling and modified Houbolt–Rainey approaches.
The aeroelastic response is categorised into four sets. The first set consists of the responses
obtained at airspeeds of 200, 225, 250 and 275ft/s, which can be considered as representing
tests at subcritical speeds. The remaining sets are around the critical flutter speed, at speeds
varying from 275 to 300ft/s in steps of 5ft/s. More specifically, the second set consists of
the responses obtained at airspeeds of 275, 280, 285 and 290ft/s, the third set at airspeeds of
275, 280, 285, 290 and 295ft/s and the fourth set at airspeeds of 275, 280, 285, 290, 295 and
300ft/s. This response segregation is done to study the sensitivity of the methods to the data
in terms of correctly predicting the flutter speed.
The test points of all four sets are presented in Table 2. Using the response at this combi-
nation of velocities, the flutter speed prediction capability of the above-mentioned methods is
demonstrated.
14 THE AERONAUTICAL JOURNAL

(a)

(b)

(c)

Figure 3. Acceleration response data. (a) Before flutter (V = 200ft/s) (b) At flutter (V = 301:68ft/s)
(c) Beyond flutter (V = 305ft/s).

5.0 FLUTTER PREDICTION USING VELOCITY DAMPING

The technique that is most commonly used to predict the onset of flutter is to extrapolate
trends of modal damping. These are data-based tools, since they rely entirely on the analysis of
flight measurements of the modal damping ratio with no assumption regarding the theoretical
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 15

Table 2
Velocity and dynamic pressure data sets

Velocity, Dynamic pressure, Velocity, Dynamic pressure,

V (ft/s) q (psi) V (ft/s) q (psi)
Set 1 Set 2
200 0.33 275 0.62
225 0.42 280 0.65
250 0.52 285 0.67
275 0.62 290 0.69
Set 3 Set 4
275 0.62 275 0.62
280 0.65 280 0.65
285 0.67 285 0.67
290 0.69 290 0.69
295 0.72 295 0.72
300 0.74

model of the specific system being tested. This approach has been used by many authors to
predict flutter(6,3,8,61,19) . The underlying idea is that the modal damping of at least one mode
of vibration becomes zero at the onset of flutter. Therefore, this flutter prediction technique
consists in noting the variation in the modal damping with airspeed, and extrapolating those
variations to an airspeed at which the damping would become zero. This resulting airspeed is
considered as the predicted flutter speed. Although the principle is quite simple, one encoun-
ters difficulty in practice. One area of difficulty is the extraction of the modal damping. Flight
test data often have low signal-to-noise ratios, hence sophisticated techniques such as param-
eter estimation or modal filtering may be required. Another area of difficulty is extrapolation.
The damping can be a highly nonlinear function of airspeed, thus the extrapolation must be
performed carefully to ensure that it accounts for any higher-order nonlinearity. Below, we
outline the approach.
The equations of motion for a two-degree-of-freedom aerofoil are given as

h h
Kij b
=ω 2
Aij + Mij b
, · · · (30)
α α

where Kij is the stiffness matrix, Mij is the mass matrix and Aij is the aerodynamic matrix.
The terms in the aerodynamic matrix are a function of the reduced frequency k. The velocity
damping or V –g technique for determining
flutter(62) assumes
a structural or hysteretic damp-
ing parameter g, such that Kij is replaced by (1 + ig) Kij . For a given reduced frequency
k = ωb
V
, Equation (30) becomes a complex eigenvalue problem

ωα2 hb hb
(1 + ig) 2 Kij = Aij + Mij . · · · (31)
ω α α

The matrices Aij , Mij and Kij are functions of the reduced frequency k. The eigenvalue,
2
λ = (1 + ig) ωω2 , for a given value of k is computed, and the ratio of the real to the imaginary
α
16 THE AERONAUTICAL JOURNAL

Table 3
Frequency and damping values at different speeds

Mode 1 Mode 2 Mode 3

Velocity Frequency Damping g Frequency Damping g Frequency Damping g
ft/s rad/s rad/s rad/s
200 53.65 −0.2037 100.50 −0.1110 342.80 −0.0154
225 55.49 −0.2669 97.22 −0.1242 342.20 −0.0168
250 57.30 −0.3616 92.89 −0.1333 341.50 −0.0180
275 58.59 −0.4900 86.58 −0.1316 340.90 −0.0191
280 58.80 −0.5234 84.86 −0.1280 340.70 −0.0193
285 58.93 −0.5401 82.91 −0.1221 340.60 −0.0195
290 59.03 −0.5694 80.39 −0.1113 340.40 −0.0197
295 59.13 −0.5954 77.00 −0.0892 340.30 −0.0198
300 59.22 −0.6240 71.97 −0.0283 340.10 −0.0199

value of the eigenvalue, g = λλRI , is varied as a function of λR = ωω2 for each reduced frequency.
2

α
The plot of g as the ordinate and V as the abscissa, for each eigenvalue, is known as the V –g
plot. The value of V at which g = 0 is the flutter speed, because the effective damping at this
airspeed is zero and the system is in a state of neutral dynamic equilibrium.
For each airspeed, the effective damping g is determined; the closer this damping is to
zero, the closer the aeroelastic system is to flutter instability. One can extrapolate the flutter
speed from off-flutter damping values and airspeeds. This is the basis of the velocity-damping
technique for predicting flutter speed.
Now, referring back to Equation 29, the complex eigenvalue problem is solved by start-
ing with large values of k; then the unsteady aerodynamic terms for the assumed value of k
are looked up, and finally the eigenvalues are solved. The reciprocal of the real part of the
eigenvalue gives the frequency of oscillation ω2 , while the imaginary part gives the damping
parameter g. Once the frequency of oscillation ω is determined, the corresponding airspeed
V is known. One then plots ω versus V and g versus V . This procedure is then repeated with
decreasing values of k until the imaginary part of the eigenvalue, or g, becomes zero. This
corresponds to the flutter condition. The corresponding real part of the eigenvalue gives the
flutter frequency ωf2 , and this together with k gives the flutter velocity Vf . The V –g method
derives its name from the graph of the total aeroelastic system damping g versus airspeed,
and from the fact that the flutter condition corresponds to the point on this graph where the
aeroelastic system damping crosses the V axis.
Flutter prediction is carried out using the different data sets, and the results are tabulated
below. The frequency and damping for various airspeeds corresponding to three modes are
presented in Table 3. From this table, it can be observed that the artificial damping values are
negative for all three modes, implying that the aeroelastic system is stable at these airspeeds.
The procedure to predict the flutter speed using the different data sets is explained below.
The damping values estimated at these data points are extrapolated to the zero damp-
ing value, and the speed corresponding to zero damping is taken as the flutter speed. The
quadratic extrapolation method is used for the extrapolation. Figure 4(a) shows the variation
of the damping with velocity for this case. From this plot, it is seen that the flutter velocity
is 410.20ft/s. This predicted flutter speed is very high compared with the theoretical value of
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 17

Figure 4. Flutter prediction with the V − g method. (a) Damping g versus velocity—data set 1 (b) Damping
g versus velocity—data set 2 (c) Damping g versus velocity—data set 3 (d) Damping g versus velocity—
data set 4.

301.68ft/s(1) . Figure 4(b) shows the variation of the damping with velocity corresponding to
data set 2. From this figure, it is seen that the flutter velocity is 316.0ft/s. This predicted flutter
speed is closer to the theoretical value of 301.68ft/s. Similar procedures are adopted for data
set 3 and data set 4 for the flutter prediction and are shown in Figure 4(c) and Figure 4(d),
respectively. Data set 3 predicts a flutter speed of 309.5ft/s, whereas data set 4 predicts a flut-
ter speed of 303.2ft/s. Therefore, the data sets at airspeeds closer to the flutter speed predict
the flutter onset speed more accurately than those that are further away.

6.0 FLUTTER PREDICTION WITH THE ENVELOPE

FUNCTION
The envelope function(11) was originally proposed as a tool to provide a quick assessment of
the overall aeroelastic stability to complement standard flutter prediction analysis. The basis
of this technique is that the impulse response of any stable damped system decays, with the
shape of the decay in the time domain being described by the decay envelope. As the damping
in a given aeroelastic system decreases, the decay envelope grows wider, eventually becoming
a rectangle as the damping becomes zero. Evaluation of the position of the centroid of the
decay envelope and the way that it shifts on the time axis as the damping decreases enables an
assessment of the stability of the system. For an aeroelastic system with an impulse response
of y(t), the decay envelope, or envelope function(11) , is given by
18 THE AERONAUTICAL JOURNAL

env(t) = (y(t))2 + (yH (t))2 , · · · (32)

where yH (t) is the Hilbert transform of the impulse response defined as

yH (t) = F −1 {Im[Y (ω)] − jRe[Y (ω)]} , · · · (33)

where Y (ω) is the Fourier transform of y(t), Im and Re are the imaginary and real part, respec-
tively, and F −1 is the inverse Fourier transform. The time centroid of the decay envelope is
given by
! tmax
env(t)tdt
t̄ = !0tmax . · · · (34)
0 env(t)dt

The upper limit of the integration, tmax , serves to define the rectangle within which the
integration takes place. For a single-degree-of-freedom system, when the damping is zero,
the time centroid lies at t = tmax /2. For a multi-degree-of-freedom system, it is suggested that
t ≈ tmax /2 is an adequate approximation for the position of the time centroid. Since t̄ tends
to increase as the damping drops, its inverse is usually employed as the significant shape
parameter, that is

1
S= . · · · (35)
t̄
In that case, the value of S at the flutter condition is S = 2/tmax for a single-degree-of-
freedom system. For a multi-degree-of-freedom system,

2
S≈ . · · · (36)
tmax

The flutter testing procedure based on the envelope function evaluates S at a number of
sub-critical airspeeds. The variation of S with airspeed is then curve-fit using a polynomial,
as is done in the case of the velocity damping technique for flutter prediction, and extrapolated
to the point where S = 2/tmax , thus yielding the flutter velocity.
Note that this tool is not intended to replace any flutter prediction methods, but rather to
provide a quick indication of whether there has been any significant change in stability since
the previous test points. In this study, the shape parameter corresponding to the four data sets
is computed and presented in Table 4.
The response from data set 1 is used to estimate the shape parameters. Usually, the shape
parameter values estimated at these data points are extrapolated to the shape parameter S =
2/tmax , and the speed corresponding to this value is taken as the flutter speed. However, in
this case, it is observed that the shape parameter values are increasing with velocity and do
not show any tendency for flutter. This is shown in Figure 5(a). Similarly for data set 2, the
shape parameters are computed and plotted with respect to velocity in Figure 5(b). From this
figure, it is seen that the shape parameters show a tendency towards flutter. The predicted
flutter speed is 303.15ft/s. The variation of the shape parameter values with velocity using
data set 3 is shown in Figure 5(c). The plot indicates a decrease in the shape parameter as
the velocity increases. This trend when extrapolated to the shape parameter S = 2/tmax thus
yields the flutter speed. The predicted flutter speed is 300.16ft/s. The response from data set 4
is used to estimate the shape parameters. These shape parameters are plotted against velocity
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 19

Table 4
Velocity and shape parameter values

Velocity, V (ft/s) Shape parameter Velocity, V (ft/s) Shape parameter

Set 1 Set 2
200 2.77 275 3.55
225 3.23 280 3.40
250 3.51 285 2.98
275 3.55 290 2.48
Set 3 Set 4
275 3.55 275 3.55
280 3.40 280 3.40
285 2.98 285 2.98
290 2.48 290 2.48
295 1.41 295 1.41
300 0.41

in Figure 5(d), showing a decreasing trend of the shape parameter with increasing velocity.
The predicted flutter speed is 300.71ft/s, and this value is close to the actual flutter speed.
Therefore, the envelope function technique can predict a flutter speed very close to the
actual flutter speed, except for data set 1.

7.0 THE HOUBOLT–RAINEY METHOD

The Houbolt–Rainey method(13,16) is a subcritical response method for flutter onset prediction
that can be used with either sinusoidally or randomly varying forced excitation. The analyt-
ical foundation of this technique is straightforward and simply requires the measurement of
the amplitude of the response in the structural mode that is important for flutter. Houbolt
and Rainey used a classical autospectrum to determine the amplitude. The amplitude values
needed by the Houbolt–Rainey method are determined by taking the r.m.s. value of the spec-
trum of the autocorrelation of the response. The auto-regressive spectrum of the system is
given by Equation (37)(36) .
tσ 2

AR (ω) = . · · · (37)

n
1+ ap eipωt

p=1

The reciprocal of these amplitude measurements is plotted against velocity, and the result-
ing curve is extrapolated to the flutter condition, which occurs when the reciprocal of the
amplitude becomes zero.
To study the effectiveness of this technique for predicting the flutter speed, the absolute
value of the reciprocal of the amplitude is computed for the different data sets and shown in
Table 2. The reciprocal of the amplitude values estimated at these data points is extrapolated
to a zero value of the reciprocal of amplitude, and the speed corresponding to this point is
taken as the flutter speed. Quadratic extrapolation is used.
The response from data set 1 is used to estimate the absolute values of the reciprocal of
amplitude. These values are plotted against velocity in Figure 6(a). From this figure, it is
observed that the amplitude values do not show any tendency towards flutter. Similarly, the
20 THE AERONAUTICAL JOURNAL

(a) (b)

(c) (d)

Figure 5. Flutter prediction with the envelope function method. (a) Envelope function method—data set 1
(b) Envelope function method—data set 2 (c) Envelope function method—data set 3 (d) Envelope function
method—data set 4.

reciprocal of the amplitude corresponding to data set 2 is estimated and plotted with respect
to velocity in Figure 6(b). From this figure, it is observed that the amplitude values show a
tendency towards flutter. The flutter speed predicted using this approach is 321ft/s. This value
shows a large deviation from the theoretical estimation of the flutter speed. For data set 3, the
reciprocal of the amplitude is estimated and plotted with respect to velocity in Figure 6(c).
From this figure, it is observed that the amplitude values show a tendency towards flutter. The
predicted flutter speed is 312.5ft/s, showing a large deviation from the theoretical estimation
of the flutter speed. Data set 4 is used to estimate the absolute value of the reciprocal of the
amplitude. The absolute value of the reciprocal of the amplitude is plotted against velocity
in Figure 6(d). The predicted flutter speed after extrapolation is 304.89ft/s. These values tend
towards the actual flutter speed of 301.68ft/s.
The Houbolt–Rainey method does not predict the exact flutter speed, even from the
response data close to the theoretical flutter speed. However, it predicts the tendency for flutter
for response data close to the actual or theoretical flutter speed.

8.0 FLUTTER PREDICTION USING THE FLUTTER

MARGIN
The flutter margin (FM) is a technique to predict the flutter onset speed(7) . This method makes
use of both modal frequency and damping information for two pre-identified critical modes
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 21

(a) (b)

(c) (d)

Figure 6. Flutter prediction with the Houbolt–Rainey method. (a) Houbolt–Rainey method—data set 1
(b) Houbolt–Rainey method—data set 2 (c) Houbolt–Rainey method—data set 3 (d) Houbolt–Rainey
method—data set 4.

to calculate the flutter margin. This technique is based on Routh’s stability criteria. Using the
flutter margin information computed for off-flutter flight conditions, it is possible to predict
the flutter onset. The basic concepts and analytical development of the flutter margin tech-
nique presented here are based on a two-degree-of-freedom aeroelastic system that represents
a bending–torsion idealisation of an aircraft wing(63) . Since most cases of flutter, even for
multi-degree-of-freedom systems, occur as a result of the interaction between two dominant
degrees of freedom, this technique is valid for multi-degree-of-freedom systems too. The sys-
tem possesses a translational inertia represented by m and a rotational inertia represented by
Ic , and is subjected to a translational elastic restraint kh and a rotational elastic restraint kα
about the elastic axis. The aerodynamic forces acting on the system are with reference to the
centre of mass of the aerofoil and consist of the lift force L and pitching moment Mc . The
aerodynamic forces are assumed to be of the form
" #
L = CLα q a1 α + a2 Vḣ + a3 Vα̇ ,
" # · · · (38)
Mc = CLα q b1 α + b2 Vḣ + b3 Vα̇ .

The lift force and pitching moment are proportional to the dynamic pressure q, angle-
of-attack α, heave velocity ḣ and pitching velocity α̇. Substituting the expression for the
aerodynamics in Equation (38) into the equations of motion in Equation (39) results in a set
of simultaneous linear differential equations in the two unknowns α and h with constant
22 THE AERONAUTICAL JOURNAL

coefficients. The constants ai and bi represent quasi-steady aerodynamic terms. The equations
of motion become

ḣ α̇
mḧ + CLα qa2 + kh h + CLα qa3 + (CLα qa1 − kh x)α = 0,
V V
· · · (39)
α̇ ḣ
Ic α̈ − CLα qb3 + (kα + kh x − CLα qb1 )α − CLα q2
2
− (kh x)h = 0.
V V

The solution to such a set can be represented in the form

h = h0 est and α = α0 est , · · · (40)

where h0 and α0 are arbitrary integration constants and s represents the eigenvalues deter-
mined by substituting Equation (40) into the simultaneous differential Equation (39). This
results in a quartic equation in s

(s4 + A3 s3 + A2 s2 + A1 s + A0 ) = 0. · · · (41)

The coefficients A0 − A3 take the form

A0 = K01 (CLα q) + K02 ,

CLα q 2 CLα q
A1 = K11 + K12 ,
V V

CLα q 2 · · · (42)
A2 = K21 + K22 CLα q + K23 ,
V

CLα q
A3 = K31 .
V

The Kij are defined as

kh (a1 x − b1 )
K01 = ,
mIc
kh kα
K02 = ,
mIc
a1 b2 − a2 b1
K11 = ,
mIc
a2 (kα + kh x2 ) + kh x(a3 − b2 ) − kh b3
K12 = ,
mIc
· · · (43)
a2 b3 + a3 b2
K21 = ,
mIc
− b1
K22 = ,
Ic

kh kα + kh x2
K23 = + ,
m Ic
a2 b3
K31 = − .
m Ic
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 23

The flutter margin is derived by applying the Routh stability criterion(45) to the two-mode
system. This criterion is algebraic and provides information on the absolute stability of a
linear time-invariant system that has a characteristic equation with constant coefficients. The
Routh–Hurwitz criterion represents an algorithm for determining the location of zeros of a
polynomial with constant real coefficients with respect to the left and right half of the s-plane
without actually solving for the zeros. In a rather direct manner, Routh’s criterion for a system
to be stable requires that all the coefficients of the characteristic equation be positive, and,
furthermore, for the specific case of a quartic equation, the following relation exist between
the coefficients:

A2 (A1 /A3 ) − (A1 /A3 )2 − A0 > 0, · · · (44)

or in a modified form
2 2
A2 A2 A1
− A0 − − > 0. · · · (45)
2 2 A3

For the system to be stable, all the A’s must be positive and the inequality of Equation (45)
must be satisfied. Instability occurs when the inequality is reversed. We can define F as a
flutter stability parameter
2 2
A2 A2 A1
F= − A0 − − . · · · (46)
2 2 A3

F must be positive for stability to exist; the more positive it is, the greater the degree of
stability. F as defined here is a measure of the system stability as expressed in the Routh
criterion. Since the A’s vary with airspeed, the value of F will also depend upon airspeed. Its
numerical value at any given speed can be considered as the flutter stability margin, and we
refer to F as the flutter margin.
We designate the four roots of the characteristic Equation (41) as s1 , s2 , s3 and s4
such that

s1,2 = β1 ± iω1 ,
s3,4 = β2 ± iω2 , · · · (47)

where the ωi ’s represent the system frequencies and the βi ’s represent the negative of the
system decay rates. The A’s can be represented in terms of the βi ’s and ωi ’s, and these are
introduced into Equation (46) to yield the flutter margin
2 2
ω2 2 − ω1 2 β2 − β1 2 ω2 2 + ω1 2 β2 + β1 2
F= + + 4β1 β2 +2
2 2 2 2
2 2
β2 − β1 ω2 − ω1 2 β2 + β1 2
− +2 . · · · (48)
β2 + β1 2 2

Using this equation, it is possible to evaluate the flutter margin corresponding to any air-
speed by measuring frequencies and damping coefficients corresponding to that speed. The
24 THE AERONAUTICAL JOURNAL

Table 5
Details of the combination of the modes

Sl. no. Mode combination

1 M1 (heaving mode) & M2 (pitching mode)
2 M1 (heaving mode) & M3 (aileron pitching mode)
3 M2 (pitching mode) & M3 (aileron pitching mode)

flutter margin equation can be used to evaluate the flutter margin at test speeds at which the
aircraft has already flown. It makes no statement about what lies beyond. Substituting the A’s
defined by Equation (42) into Equation (46), the resulting equation takes the form

F = B2 (CLα q)2 + B1 (CLα q) + B0 , · · · (49)

where the coefficients are

K11 K11 ρCLα K11 K21
B2 = + K22 + ,
K31 K31 2 K31

1 K11 K12 ρCLα K12 K21
B1 = K12 K22 + K11 K23 − 2 − K01 + , · · · (50)
K31 K31 2 K31

K12 K12
B0 = K23 − − K02 .
K31 K31

Equation (49), expressing the parabolic variation of F with CLα q, is referred to as the
flutter prediction equation. The coefficients B0 , B1 and B2 are functions of aeroelastic param-
eters. The effect of CLα is captured by using flight test data corresponding to the subsonic or
supersonic flight test regimes. The Bi ’s are evaluated after inverting Equation (49) from the
measured off-flutter test points for which the flutter margin F, CLα and q are known. Once the
coefficients are determined, the CLα q corresponding to F = 0 is determined by extrapolation.
In the case of the flutter margin technique, the modal combination plays a critical role in
predicting the flutter speed. By selecting the combination of modes carefully, one can pre-
dict the onset of flutter speed accurately. The flutter margin is estimated using the frequency
and damping values for the first three data sets as mentioned for the other flutter prediction
techniques. However, the flutter margin is not estimated for the fourth set, as the other three
sets are able to demonstrate the flutter prediction efficiently. To determine the flutter onset
speed, the modes must be combined in pairs. The modal pair showing the lowest flutter mar-
gin is considered as to be the critical modal combination. The different modal combinations
considered in this study are presented in Table 5.
The frequency and damping estimations corresponding to data set 1 are used to compute
the flutter margin. The flutter margin estimated using data set 1 is plotted against velocity
in Figure 7 for all three mode combinations. It is observed that the pitch and heave mode
combination is the critical mode combination, and it exhibits a classical flutter phenomenon.
The second modal combination, where the wing pitching mode interacts with the control
surface mode, also follows the familiar control surface flutter mode combination. The third
modal combination of heaving with the control surface mode does not produce any flutter
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 25

(a) (b)

(c) (d)

(e) (f)

Figure 7. Flutter prediction with the flutter margin method. (a) Variation of flutter margin against velocity
(b) Flutter onset velocity for a critical mode combination (c) Variation of flutter margin against velocity (d)
Flutter onset velocity for a critical mode combination (e) Variation of flutter margin against velocity (f) Flutter
onset velocity for a critical mode combination.

instability. It is seen that the flutter margin shows an instability trend from the starting velocity.
The critical mode combination which leads to flutter is shown in Figure 7(b). The flutter onset
speed is estimated to be 299.50ft/s.
Similarly, the frequency and damping estimations corresponding to data set 2 are used to
compute the flutter margin, which is plotted against velocity in Figure 7(c). As predicted in
the previous case, the pitch and heave mode combination is the critical mode combination. It
is observed that the flutter margin shows instability from the beginning. However, the velocity
26 THE AERONAUTICAL JOURNAL

damping approach does not show any instability or reduction in damping until a velocity of
280ft/s is reached, for the wing pitching mode. The critical mode combination which leads
to flutter is shown in Figure 7(d). The flutter onset speed is estimated to be 302.50ft/s. The
flutter margin estimated using data set 3 is shown in Figure 7(e). As in the previous two cases,
it is found that the pitch and heave mode combination is the critical mode combination. The
flutter onset speed (Figure 7(f)) of 301.4ft/s estimated using this data set is very close to the
theoretical flutter speed. The flutter margin method is thus able to predict the onset of flutter
even when using the data set lying far below the flutter onset speed. Hence, the data set closer
to the flutter onset is redundant and data set 4 is not used.

9.0 FLUTTER PREDICTION USING THE

AUTO-REGRESSIVE MODEL
Auto-regressive models are finite-difference models. As the flutter instability prediction relies
on time-domain measurements of the system response, auto-regressive models are suitable as
they are constructed from discrete-time response sequences. The stability in the constructed
auto-regressive model is based on Jury’s stability criteria. The stability boundary is estimated
using Jury’s stability determinants, which are expressed in terms of auto-regressive coeffi-
cients. The auto-regressive method for flutter prediction consists of two parts. The first part is
the identification of the aeroelastic system as an auto-regressive model from the sampled time
response and estimation of parameters. The second part consists of an estimation of the flut-
ter prediction parameter based on Jury’s stability criteria using the estimated auto-regressive
coefficients.
The auto-regressive model identification consists of three steps, namely the selection of a
suitable model for the data, the selection of a suitable model order for the model, and the
estimation of parameters. An AR process of order p can be represented as

p
x[n] = − a[k]x[n − k] + u[n] · · · (51)
k=1

This process is termed auto-regression as the sequence x [n] is a linear regression on

itself with u[n] representing the error. Using this model, the present value of the process
is expressed as a weighted sum of past values plus a noise term. The parameters a[k] are
termed the auto-regressive coefficients.
The z transform of the auto-regressive model can be represented as

p
A(z) = a[k]z−k = 1 + a[1]z−1 + a[2]z−2 + ... + a[p]z−p . · · · (52)
k=0

The selection of the model order in the AR spectral estimation is a crucial step. If the
selected model order is too low, it results in a smoothed estimate, whereas if it is too large, it
causes spurious peaks and general statistical instability. All model order estimators are based
on the estimated prediction error power. The estimated prediction error power is guaranteed
to decrease or remain the same as the model order increases for all AR parameter estimation
methods. Two methods are the Final Prediction Error (FPE) and Akaike Information Criterion
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 27

(AIC)(40,41,37,30,39) . The optimised model order is the one with the minimum FPE, defined as

N +K
FPE(K) = ρk , · · · (53)
N −K

where ρk is the estimate of the white noise variance, or the prediction error power, for the
k th -order auto-regressive model and N is the number of data points. It is seen that, although ρk
decreases with k, the term N+KN−K
increases with k. The final prediction error is an estimate of
the prediction error power when the prediction coefficients must be estimated from the data.
The term N+K N−K
accounts for the increase in the variance of the prediction power estimator due
to inaccuracies in the prediction coefficient estimates.
A second criterion is the Akaike Information Criterion (AIC), which is defined as

AIC(K) = N ln ρk + 2k. · · · (54)

The model order selected is the one that minimises the AIC. The performance of the Akaike
information criterion and final prediction error is similar. For short data records, the use of the
Akaike information criterion is recommended. For larger data records, both estimators will
yield identical model order estimates, since they are functionally related to each other. An
important feature of the auto-regressive model description is that it has predictor capability in
that, given the history of the system, the subsequent output can be predicted.
There is a direct correspondence between the auto-regressive parameters and the covariance
function of the response, and this correspondence can be inverted to determine the auto-
regressive parameters or coefficients from the autocorrelation function. This is done using the
Yule–Walker equations. Multiplying Equation (51) by x[n − k] and taking the expectation of
the product gives

p
rxx [k] = − a[l]rxx [k − l] + rxu [k − l], · · · (55)
l=1

where rxu [k] = E(u[n]x[n − k]). But rxu = 0 for k > 0 and rxu = σ 2 for k = 0. Therefore
⎧
⎪
⎪
p
⎪
⎪ − a[l]rxx [k − l] for k ≥ 1,
⎨
rxx [k] = l=1
· · · (56)
⎪
⎪ p
⎪
⎩−
⎪ a[l]rxx [−l] + σ 2 for k = 0.
l=1

Equations (56) are the Yule–Walker equations. They define a nonlinear relationship
between the parameters of an AR process and the autocorrelation function. The Yule–Walker
equations in matrix form are given by
⎡ ⎤
rxx [0] rxx [−1] ··· rxx [−(p − 1)] ⎡a[1]⎤ ⎡
rxx [1]
⎤
⎢ ⎥
⎢ rxx [1] rxx [0] ··· rxx [−(p − 2)]⎥ ⎢ a[2]⎥ ⎢r [2]⎥
⎢ ⎥⎢ ⎥ ⎢ xx ⎥
⎢ .. .. .. .. ⎥⎢ . ⎥=−⎢
⎢ ⎥
⎢ .. ⎥
⎥ · · · (57)
⎢ . . . . ⎥ ⎣ .. ⎦ ⎣ . ⎦
⎣ ⎦
rxx [p − 1] rxx [p − 2] ··· rxx [0] a[p] rxx [p]
28 THE AERONAUTICAL JOURNAL

In this study, the auto-regressive parameters are estimated using the MATLAB system
identification toolbox(64) .
The flutter prediction parameter based on Jury’s stability criteria is calculated algebraically
from the auto-regressive coefficients of the estimated auto-regressive model. The system
under consideration is stable if Jury’s stability criterion(65) is satisfied. A discrete-time sys-
tem is stable if and only if all the characteristic roots are located inside the unit circle in the
Gauss plane. Jury’s determinant method judges the stability of a discrete-time system based
on the characteristic polynomial without solving for the roots. Suppose that the characteristic
equation of the discrete system is given by

G(z) = An zn + An−1 zn−1 + . . . + A1 z + A0 , · · · (58)

where n is an even number and Ai are the coefficients of the characteristic polynomial. Then
the necessary and sufficient conditions for the system to be stable are that all of the following
parameters are positive:

G(1) = An + An−1 + · · · + A1 + A0 ,
G(−1) = An − An−1 + · · · + A1 + A0 ,
F + (k) = det(Xk + Yk ) k = 1, 3, · · · , n − 1, · · · (59)
F − (k) = det(Xk − Yk ) k = 1, 3, · · · , n − 1

where Xk and Yk are matrices whose elements consist of the coefficients of Equation (58)
⎡ ⎤ ⎡ ⎤
An · · · ··· An−k+1 Ak−1 ··· A0
Xk = ⎣ 0 ··· ··· ⎦, Yk = ⎣ · · · ··· 0 ⎦ · · · (60)
0 0 An , A0 0 0.

Among these parameters, F − (n − 1) can be expressed by means of characteristic roots as

$
F − (n − 1) = An−1
n (1 − zi zj ) = det(Xn−1 − Yn−1 ). · · · (61)
i<j

Therefore, F − (n − 1) becomes zero whenever any pair of complex-conjugate characteristic

roots reaches the unit circle, which is the stability boundary.
F − (n − 1) is called the stability parameter and is utilised for flutter prediction. Although
F − (n − 1) is a better index for flutter prediction than the modal damping, its behaviour against
dynamic pressure is inferior to that of the flutter margin. To overcome this drawback, Torii(47)
proposed a flutter prediction parameter that is linearly related to the dynamic pressure. To
predict the flutter boundary of a multi-mode system, the flutter prediction parameter FN is

F − (n − 1)
FN = , · · · (62)
(F − (n − 2))2

where F − (n − 1) and F − (n − 2) are defined by Equations (61). For a binary flutter system, a
fourth-order characteristic polynomial

G(z) = A4 z4 + A3 z3 + A2 z2 + A1 z + A0 , · · · (63)
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 29

should be estimated, where A0 , A1 , A2 , A3 and A4 are the auto-regressive coefficients. Jury’s

stability criterion simplifies to

G(1) = A4 + A3 + A2 + A1 + A0 > 0,
G(−1) = A4 − A3 + A2 − A1 + A0 > 0,
F + (1) = A4 + A0 > 0,
· · · (64)
F − (1) = A4 − A0 > 0,
F + (3) = |X3 + Y3 | > 0,
F − (3) = |X3 − Y3 | > 0,

with
⎡ ⎤ ⎡ ⎤
A4 A3 A2 A2 A1 A0
X3 = ⎣ 0 A4 A3 ⎦, Y3 = ⎣ A1 A0 0 ⎦ · · · (65)
0 0 A4 , A0 0 0.

For a binary flutter system, the number of modes nm = 2 and the flutter prediction parameter
Fz is defined by

F − (2nm − 1) F − (3) |X3 − Y3 |

Fz = = = · · · (66)
(F − (2nm − 3)) 2
(F − (1))2 (A4 − A0 )2

and

Fz = C1 q2 + C2 q + C3 · · · (67)

The flutter boundary is predicted as the dynamic pressure at which the flutter prediction
parameter FZ becomes zero.
The response at each point of data set 1 is used to find the stability parameter by modelling
the response using the auto-regressive model. Figure 8(a) shows the variation of the flutter
stability parameter with velocity. It is observed that the flutter stability parameter shows a
decreasing trend from the beginning of the curve, as in the case of the flutter margin technique.
The stability parameter is positive for all the speeds considered in this case.
The response for data set 2 is modelled using the auto-regressive model to estimate the
stability parameters. Figure 8(b) shows the variation of the flutter stability parameter with
velocity. The decreasing trend of the flutter stability parameter curve is extrapolated to obtain
the flutter speed, which is found to be 326ft/s.
Using data set 3, the stability parameters are estimated and shown in Figure 8(c). The
predicted flutter speed is 316ft/s, which is closer to theoretical flutter speed as compared with
those predicted using the previous data sets.
When using data set 4, where the velocities are closer to the flutter speed, the estimated
flutter stability parameters are plotted against velocity in Figure 8(d). The flutter speed after
extrapolation is found to be 303ft/s, which is close to the theoretical flutter speed among all
the data sets considered for this model.
30 THE AERONAUTICAL JOURNAL

(a) (b)

(c) (d)

Figure 8. Flutter prediction with the auto-regressive technique. (a) Auto-regressive method—Data set 1 (b)
Auto-regressive method—Data set 2 (c) Auto-regressive method—Data set 3 (d) Auto-regressive method—
Data set 4.

10.0 SUMMARY OF RESULTS

Table 6 summarises the flutter speeds estimated using the various prediction tech-
niques, namely the velocity damping, envelope function, Houbolt–Rainey, Zimmerman–
Weissenburger flutter margin and discrete-time auto-regressive modelling approaches, for
the four different cases (data sets 1 to 4). The velocity damping approach predicts a flutter
speed when the data set considered is near the critical speed. However, when the data set
considered is far from the critical speed, the flutter speed prediction error increases. Similarly,
the envelope function technique calculates the overall damping trend by using the shape
parameter. However, the effectiveness of this approach improves when the data set is close
to the flutter speed. It is observed that the Houbolt–Rainey method predicts flutter only
at speeds very close to flutter. So this technique seems to be ineffective for predicting the
flutter onset speed. The flutter margin method predicts the most accurate flutter speed, even
when the data set considered for flutter prediction is far from the actual flutter speed. The
decrease in the overall stability of the system as the speed increases can be monitored clearly,
providing an indication towards flutter even at speeds much lower than the actual flutter
speed. The autoregressive approach also shows a decreasing trend long before flutter starts.
Table 7 provides a comparison of the flutter prediction methods.
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 31

Table 6
Comparison of flutter prediction methods results

Data sets
Method Set 1 Set 2 Set 3 Set 4
Vfpa (ft/s) % Errorb Vfp (ft/s) % Error Vfp (ft/s) % Error Vfp (ft/s) % error
Vgc 410.20 35.97 316.00 4.75 309.56 2.61 303.20 0.51
EFd – – 303.15 0.49 300.16 0.50 300.71 0.32
HRe – – 321.00 6.40 312.50 3.59 304.89 1.07
FMf 299.50 0.72 302.50 0.27 301.40 0.09 – –
ARg – – 326.00 8.06 316.00 4.75 303.00 0.44
a Predicted flutter speed; b error = ((Vfp − Vf )/Vf ) × 100, where the theoretical flutter speed Vf = 301.68ft/s;
c velocity damping; d envelope function; e Houbolt–Rainey; f flutter margin; g auto-regressive model.

Table 7
Comparison of flutter prediction methods

Effectiveness
Data close Data far away
to critical from critical
Sl. No. Method Parameter speed speed
a
1. Vg Damping Yes No
2. EFb Shape Yes No
parameter
3. HRc 1/| amplitude| Yes No
4. FMd Flutter Yes Yes
margin
5. ARe Flutter Yes Yes
stability
parameter
a Velocity damping; b Envelope function; c Houbolt–Rainey; d Flutter margin; e Auto-regressive
model.

11.0 CONCLUSION
A three-degree-of-freedom model of a typical wing section is used to evaluate various flut-
ter prediction techniques, namely the velocity damping, envelope function, Houbolt–Rainey,
Zimmerman–Weissenburger flutter margin and discrete-time auto-regressive modelling
approaches.
From Table 7, it is clear that the flutter margin and discrete-time auto-regressive mod-
elling approaches are automatic choices for flutter prediction. However, the velocity damping
method is still followed in industrial practice. Our recommendation is to use either the flutter
margin or auto-regressive method, as they are robust and can also be used with data points
obtained far before the onset of flutter, which is an added advantage.
32 THE AERONAUTICAL JOURNAL

SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit https://doi.org/10.1017/aer.
2020.84

REFERENCES
1. NAM, C., KIM, Y. and WEISSHAAR, T.A. Computational Aids in Aerservoelastic Analysis using
MATLAB, Tata McGraw-Hill Publishing Company Limited, 2001.
2. HAYES, W.B., GOODMAN, C.E. and SISK, K.D. F/A-18 E/F Super Hornet flutter clearance pro-
gram, No. AIAA Paper 2003-1940, AIAA/ASME/ASCE/AHS Structures, Structural Dynamics and
Materials Conference, Norfolk, Virginia, USA, April 2003.
3. KEHOE, M.W. A historical overview of flight flutter testing, Tech. Rep. NASA TM 4720, NASA,
Washington, DC 20546-0001, May 8–10 1995.
4. MIL-A-8870C, Military Specification - Airplane Strength and Rigidity, Vibration, Flutter, and
Divergence, Air Force Flight Dynamics Laboratory (AFFDL), 1993.
5. DESMARAIS, R.N. and BENNETT, R.M. An automated procedure for computing flutter eigenvalues, J.
Aircr., February 1974, 11, (2), pp 75–83.
6. KOENIG, K. Flight vibration test analysis–Methods, theory and application, AIAA/AHS/
IES/SetP/SFTE/DGLR Flight testing conference, No. AIAA Paper 1983-2752, Las Vegas,
Nevada, November 16–18 1983, pp 1–20.
7. ZIMMERMAN, N.H. and WEISSENBURGER, J.T. Prediction of flutter onset speed based on flight testing
at subcritical speeds, J. Aircr., July-August 1964, 1, (4), pp 190–202.
8. COOPER, J.E. Parameter estimation methods for flight flutter testing, 80th AGARD Structures
and Materials Panel, No. AGARD CP-566, Rotterdam, The Netherlands, 8-10 May 1995, pp
10–1–10–11.
9. KAYRAN, A. Flight flutter testing and aeroelastic stability of aircraft, Int. J. Aircr. Eng. Aerosp.
Technol., 2007, 79, (5), pp 494–506.
10. DAWSON, K.S. and MAXWELL, D.L. Limit cycle oscillation flight-test results for asymmetric store
configuration, J. Aircr., November-December 2005, 42, (6), pp 1589–1596.
11. COOPER, J.E., EMMET, P.R. and WRIGHT, J.R. Envelope function – A tool for analyzing flutter data,
J. Aircr., September–October 1993, 30, (5), pp 785–790.
12. DIMITRIADIS, G. and COOPER, J.E. Flutter prediction from flight flutter test data, J. Aircr., 2001, 38,
(2), pp 355–367.
13. HOUBOLT, J.C. and RAINEY, A.G. On the prediction of critical flutter conditions from subcritical
response data and some related wind-tunnel experience, Proceedings of the 1958 Flight Flutter
Testing Symposium, No. NASA-SP-385, Washington, D.C., May 15–16 1958, pp 23–29.
14. SANDFORD, M.C., ABEL, I. and GRAY, D.L. Development and demonstration of a flutter-suppression
system using active controls, Tech. Rep. NASA TR R-450, NASA Langley Research Center,
Washington D.C., 1975.
15. FOUGHNER, J.T.J. Some experience using subcritical response methods in wind-tunnel flutter model
studies, Flutter testing techniques, No. NASA SP-415, NASA, Washington D.C, October 1976.
16. DOGGETT, R.V.J. Some observations on the houbolt-rainey and peak-hold methods of flutter onset
prediction, Tech. Rep. NASA TM-102745, NASA, Langley Research Center, Hampton, Virginia,
23665-5225, November 1990.
17. BRENNER, M.J., LIND, R.C. and VORACEK, D.F. Overview of recent flight flutter testing research at
NASA Dryden, 38th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material
Conference and Exhibit, Kissimmee, Floida, No. NASA TM-4792, NASA Dryden Flight
Research Center, Kissimmee, Floida, April 7–10 1997.
18. LIND, R. A flight test to demonstrate flutter and evaluate the flutterometer, Aeronaut. J., June 2003,
107, pp 577–588.
19. LIND, R. Flight-test evaluation of flutter prediction methods, J. Aircr., September–October 2003,
40, (5), pp 964–970.
20. LIND, R. Flight testing with the flutterometer, J. Aircr., May–June 2003, 40, (5), pp 574–579.
SUDHA ET AL ASSESSMENT OF FLUTTER PREDICTION TECHNIQUES... 33

21. LIND, R. and BRENNER, M. Robust flutter magins of an F/A -18 aircraft from aeroelastic flight data,
J. Guid. Control Dyn., May–June 1997, 20, (3), pp 597–604.
22. LIND, R. and BRENNER, M. Incorporating flight data into a robust aeroelastic model, J. Aircr.,
May–June 1998, 35, (3), pp 470–477.
23. LIND, R. and BRENNER, M. Flutterometer: An on-line tool to predict robust flutter margins, J. Aircr.,
November–December 2000, 37, (6), pp 1105–1112.
24. PRICE, S.J. and LEE, B.H.K. Evaluation and extension of the flutter-margin method for flight flutter
prediction, J. Aircr., May–June 1993, 30, (3), pp 395–402.
25. LIND, R. Flutter margins for multi-mode unstable couplings with associated flutter confidence,
J. Aircr., September–October 2009, 46, (5), pp 1563–1568.
26. TORII, H. and MATSUZAKI, Y. Flutter margin evaluation for discrete-time systems, J. Aircr.,
January–February 2001, 38, (1), pp 42–47.
27. KATZ, H., FOPPE, F.G. and GROSSMAN, D.T. F-15 flight flutter test program, Flutter Testing
Techniques, No. NASA SP-415, Dryden Flight Research Center, 1976, pp 413–431.
28. ASTROM, K.J. and EYKHOFF, P. System identification - A survey, Automatica, 1971, 7, pp 123–162.
29. LJUNG, L. System Identification - Theory for the user, Prentice Hall PTR, 1999, New Jersey.
30. SUNDARAMURTHY, H., JATEGAONKAR, R.V. and BALAKRISHNA, S. Dynamic stability measurements
from tunnel unsteadiness excited random response, J. Aircr., January 1980, 17, (1), pp. 7–12.
31. PERANGELO, H.J. and WAISANEN, P.R. Application of advanced parameter identification meth-
ods for flight flutter data analysis with comparisons to current techniques, AGARD Conference
Proceedings: Flight Test Techniques, Vol. 29, Grumman Aerospace Corporation, Calverton, New
York, July 1984, pp 1–28.
32. BATILL, S.M., CAREY, D.M. and KEHOE, M.W. Digital time series analysis for flutter test data,
AIAA Dynamics Specialist Conference, Dallas, Texas, No. AIAA Paper 92-2103-CP, April 1992,
pp 215–223.
33. PINKELMAN, J.K., BATILL, S.M. and KEHOE, M.W. An investigation of the total least
square criteria in time domain based, parameter identification for flight flutter testing,
AIAA/ASME/ASCE/AHS/ASC 36th Structures, Dynamics and Materials Conference, No. AIAA
Paper 95-1247-CP, New Orleans, Louisiana, April 1995, pp 783–793.
34. PINKELMAN, J.K., BATILL, S.M. and KEHOE, M.W. Total least square criteria in time domain based,
parameter identification for flight flutter testing, J. Aircr., July–August 1996, 33, (4), pp 784–792.
35. MCNAMARA, J.J. and FRIEDMANN, P.P. Flutter-boundary identification for time-domain computa-
tional aeroelasticity, AIAA J., 2007, 45, (7), pp 1546–1555.
36. JATEGAONKAR, R.V. and BALAKRISHNA, S. Autoregressive modeling and damping evaluation from
random excitation response of structures, Proceedings, National System Conference, India, 1977,
pp B1–1–5.
37. KAY, S.M. Modern Spectral Estimation, Prentice-Hall Inc., 1988, Englewoods Cliff, New Jersey.
38. KAY, S.M. and MARPLE, S.L. Spectrum analysis - A modern prospective, Proc. IEEE, November
1981, 69, pp 1380–1419.
39. RAOL, J.R., JATEGAONKAR, R.V. and BALAKRISHNA, S. Determination of model order for dynamical
system, IEEE Trans. Syst. Man Cybern., January/February 1982, SMC-12, (1), pp 56–62.
40. AKAIKE, H. Statistical predictor identification, Ann. Inst. Stat. Math., 1970, 22, pp 203–217.
41. AKAIKE, H. A new look at the statistical model identification, IEEE Trans. Autom. Control,
December 1974, AC-19, (6), pp 716–723.
42. PAN, W. Akaike’s information criterion in generalized estimating equations, Biometrics, March
2001, 57, (1), pp 120–125.
43. BOZDOGAN, H. Akaike’s information criterion and recent developments in information complexity,
J. Math. Psychol., 2000, 44, pp 62–91.
44. MATSUZAKI, Y. and ANDO, Y. Estimation of flutter boundary from random responses due to
turbulence at subcritical speeds, J. Aircr., 1981, 18, (10), pp 862–868.
45. KUO, B.C. Automatic Control System, Prentice-Hall Inc., Englewoods Cliff, 1987, New Jersey.
46. BAE, J.S., MATSUZAKI, Y. and INMAN, D.J. Extension of flutter prediction parameter for multimode
flutter system, J. Aircr., 2010, 42, (1), pp 285–288.
47. TORII, H. and MATSUZAKI, Y. Flutter margin evaluation for three-mode discrete-time sytems, 52th
AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Material Conference, No.
AIAA Paper 2011-2144, 4-7 April 2011, pp 1–8.
34 THE AERONAUTICAL JOURNAL

48. THEODORSEN, T. General theory of aerodynamic instability and the mechanism of flutter, Tech.
Rep. NACA Technical Report 496, NACA, 1935.
49. VEPA, R. Finite State Modeling of Aeroelastic Systems, Tech. Rep. NASA CR-2779, NASA, 1977.
50. EDWARDS, J. Applications of Laplace transform methods to airfoil motion and stability calcula-
tions, AIAA J., 1979, (79-0772), pp 465–481.
51. KARPEL, M. Design for active flutter suppression and gust alleviation using state-space aeroelastic
modeling, J. Aircr., 1982, 19, (3), pp 221–227.
52. KARPEL, M. Design for active and passive flutter suppression and gust alleviation, Tech. Rep.
NASA CR-3482, NASA, 1981.
53. KARPEL, M. and STRUL, E. Minimum-state unsteady aerodynamic approximations with flexible
constraints, J. Aircr., 1996, 33, (6), pp 1190–1196.
54. WANG, Z. and CHEN, P.C. Adapted K-method for frequency-domain ASE control stability margin
analysis, AIAA/ASMe/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference,
National Harbor, Maryland, January 2014.
55. KARPEL, M. Extensions to the minimum-state aeroelastic modeling method, AIAA J., 29, (11),
2008, pp 2007–2009.
56. BISKRI, D., BOTEZ, R.M., STATHOPOULOS, N., THERIEN, S., RATHE, A. and DICKINSON, M. New mixed
method for unsteady aerodynamic force approximations for aeroservoelasticity studies, AIAA J.,
2006, 43, (5), pp 1538–1542.
57. DINU, A., BOTEZ, R.M. and COTOI, I. Chebyshev polynomials for unsteady aerodynamic calcula-
tions in aeroservoelasticity, J. Aircr., 2006, 43, (1), pp 165–171.
58. BOTEZ, R.M., DINU, A. and COTOI, I. Method based on Chebyshev polynomials for aerservoelastic
interactions on an FA-18 aircraft, J. Aircr., 2007, 44, (1), pp 330–333.
59. ROGER, K.L. Airplane math modeling methods for active control design, Tech. Rep. AGARD-CP-
228, AGARD, 1977.
60. WRIGHT, J.R. and COOPER, J.E. Introduction to Aircraft Aeroelasticity and Loads, John Wiley and
Sons, 2007, U.K.
61. ELDRED, M.S., VENKAYYA, V.B. and ANDERSON, W.J. New mode tracking methods in aeroelastic
analysis, AIAA J., July 1995, 33, (7), pp 1292–1299.
62. FUNG, Y.C. An Introduction to the Theory of Aeroelasticity, John Wiley & Sons, Inc., 1955, New
York.
63. HANCOCK, G.J., WRIGHT, J.R. and SIMPSON, A. On the teaching of the principles of wing flexure-
torsion flutter, Aeronaut. J., October 1985, 89, (888), pp 285–305.
64. LJUNG, L. System Identification Tool Box - For use with Matlab, Version 5, The Mathworks, Inc.,
2002, 3 Apple Hill Drive, Natick, MA, USA.
65. JURY, I.E. Theory and Application of the Z-Transform Method, Wiley, 1964, New York.