Table Of ContentHigher-Order Spectral Analysis of
F-18 Flight Flutter Data
Walter A. Silva∗
NASA Langley Research Center, Hampton, Virginia, 23681
Shane Dunn†
Defence Science and Technology Organisation, Department of Defence, Australia
Royal Australian Air Force (RAAF) F/A-18 flight flutter test data is presented and
analyzed using various techniques. The data includes high-quality measurements of forced
responses and limit cycle oscillation (LCO) phenomena. Standard correlation and power
spectraldensity(PSD)techniquesareappliedtothedataandpresented. Novelapplications
ofexperimentally-identifiedimpulseresponsesandhigher-orderspectraltechniquesarealso
appliedtothedataandpresented. Thegoalofthisresearchistodevelopmethodsthatcan
identify the onset of nonlinear aeroelastic phenomena, such as LCO, during flutter testing.
I. Introduction
In recent years, there has been a a dramatic increase in the study of nonlinear aeroelastic phenomena.
The recent survey paper by Dowell, Edwards, and Strgnac1 discusses the broad range of computational and
experimental activities aimed at understanding and addressing various types of nonlinear phenomena. This
paper, however, will focus on the analysis of data acquired during flight flutter testing of RAAF F/A-18
aircraft,2 which includes high-quality measurements of limit cycle oscillation (LCO) phenomena.
During full-scale development of the F/A-18A/B Hornet, it was discovered that when carrying heavy
stores on the outboard wing pylons and AIM-9 missiles on the wing tips, the aircraft experiences an “unac-
ceptable5.0-6.0Hzoscillationatlowaltitudeandhighspeed”.3 Thisoscillationisasingle-degree-of-freedom
LCO.
An aeroelastic LCO is the result of nonlinear dynamics in the structural system or the aerodynamic
system, or both. Quite often, as the amplitude of an initial flutter incident increases, the nonlinearities
associated with the aeroelastic system become significant. As these nonlinearities grow in magnitude, the
aeroelastic system is transformed from an initially unstable linear flutter event into a limited amplitude
nonlinear aeroelastic oscillation, or an LCO. Denegri4 discusses the different types of LCO encountered by
the F-16 aircraft during its flight flutter testing phase.
RecentapplicationsoftheVolterratheorytoexperimentalaerodynamicsandaeroelasticityareproviding
valuableknowledgeregardingnonlinearaeroelasticbehavior. Inparticular,theexperimentalidentificationof
aerodynamic impulse responses may provide insight regarding the dominant flow physics of the experiment
as well as an automatic data filtering capability.5 The application of higher-order spectra (HOS) to flutter
data6,7 and the identification of Volterra kernels from flight flutter experiments8–10 are additional examples
of this promising application for the Volterra theory.
The paper is organized as follows. A description of the flight test instrumentation and data acquisition
process is presented. This is followed by a description of the various types of data acquired and the data
that is analyzed in the paper. A description of higher-order spectral methods is then provided. Finally, the
results of the analysis of the data using various techniques is presented and discussed.
∗SeniorResearchScientist,AeroelasticityBranch,AIAAAssociateFellow.
†Head-Aeroelasticity,AirVehiclesDivision.
1of20
AmericanInstituteofAeronauticsandAstronautics
II. Flight Test Instrumentation and Data Description
Greater detail on the aircraft configuration and instrumentation for these tests is given in Arms et al11
and Keeler et al.12 A brief summary is given here. The schematic in figure 1 shows the various paths
taken by the data. For safety while the aircraft is on-condition, strip charts are monitored to ensure that
pre-determined acceleration limits are not exceeded at key locations, and that accelerations do not grow in
an unexpected manner during, or following, FECU runs. The FECU is Flutter Exciter Controller Unit and
allows the ailerons to be driven in sine dwell, sine sweep or random modes.
Figure 1. Schematic of the data acquisition process.
Thefrequenciesfortheinputsinedwells,usingtheFECUateachstabilizedflighttestpoint,werechosen
in order to excite the anti-symmetric wing first bending (AW1B) mode. Data from key accelerometers
were analyzed and the frequency and damping for the AW1B mode were estimated. In this way, damping
estimates were generated so that any dangerous trends could be assessed during the flight. If the damping
estimate during a flight dropped to less than 1.5 percent of critical, then any additional planned envelope
expansion points at similar or greater Mach no./dynamic pressure would not be attempted. If there was
sufficient fuel remaining, other pre-briefed test points at safer parts of the envelope would be flown.
TwoMATLABprogramsweredevelopedforon-conditionandpost-flightdataanalyses. Theseprograms
resulted in graphical user interfaces (GUIs) which enabled very fast analysis and visual inspection of the
quality of the measured data. The first program analyzed the data in the time-domain. An example of the
GUI for this program is shown in figure 2. In figure 2, accelerometer locations also are presented. A number
of time domain analysis approaches were available in this program, but the most reliable results were found
by examining the decay following a FECU dwell using the Eigensystem Realization Algorithm with Data
Correlation (ERA/DC) (Juang et al13). For the rare occasions where stick-raps were used, rather than the
FECU, random decrement processing (Cole14) was used, followed by ERA/DC.
For sweeps using the FECU (typically 2-11Hz over 30 sec.), a transfer function analysis was performed
in the frequency-domain using an output-error curve fit in the complex plane (Ljung15). This gave an
estimateforthemodalfrequenciesanddampingsinthefrequencyrangeofinterest. Suchatransferfunction
analysis, however, requires a reference signal that is linearly proportional to the excitation force. The best
reference signals available were aileron rotational position and aileron rotational rate. However, neither of
these signals is a true force reference. For this reason, these data were only used in a qualitative sense to
provide confirmation of trends found from the analysis of the FECU dwell data. The sweep data were also
used to ensure that no instability at an unexpected frequency outside the range of the FECU dwells was
approaching.
Between flights, stability margin estimates were generated using the two-degree-of-freedom method of
Zimmerman and Weissenburger.16 In this case, the two-degree-of-freedom system consisted of the AW1B
2of20
AmericanInstituteofAeronauticsandAstronautics
Figure 2. MATLAB Graphical User Interface (GUI) for time-domain analyses of measured data.
and anti-symmetric outboard store pitch (AOBSP) modes. Using this method, extrapolations were made to
estimate the altitude at which the system would become unstable. Such stability margin estimates are only
valid within regimes where the structural and aerodynamic forces behave linearly; therefore, analysis across
test points where there is a significant variation in compressibility and/or shock effects is not valid. For the
Mach numbers tested here, it was found that data for M=0.85 formed a self consistent set of results, M=0.9
another, and data within the range of M=0.93 to M=0.95 gave a third set of results.
Presented in figure 3 is a sample of the type of data acquired during these flight tests. Figure 3 presents
aileron sine dwell inputs and resultant accelerations at the aft launcher location during a test run that
resulted in a mild flutter condition. As can be seen, the level of positive damping (stable) decreases until
negative damping (unstable) is encountered. Quite often during flight flutter testing, LCO phenomena
are encountered. These LCO phenomena can mask strong (and dangerous) flutter phenomena which can
resurface when the nonlinearity that is creating the LCO disappears. In addition to sine dwell and sine
sweep excitations, data acquired during lateral stick raps will also be presented.
Theaircraft’sAW1Bmodeexhibitsalargedecreaseinstiffnessoveraportionoftheflightenvelope. This
nonlinear and temporary reduction in stiffness creates the illusion of a highly damped (i.e., stable) system.
But when flying into another portion of the flight envelope, the wing’s AW1B mode stiffness returns to its
normal (linear) value and the aircraft may encounter an unsafe flutter condition. A contour map of the
damping that clearly delineates these low damping (L) and high damping (H) regions of the flight envelope
is presented as figure 4.
Clearly, a strong motivation exists to develop methods that can identify regions of nonlinear behavior
prior to their effects becoming significant. With this in mind, the present research effort has focused on the
analysis of the flight flutter data using novel spectral analysis methods in addition to the more traditional
approaches.
A very large amount of data was acquired during these flight tests. Based on the signal-to-noise ratio of
the data and the clearly identified data points, the data can classified as high-quality data. As the intent of
this paper is to introduce the application of HOS to flight test data, the results presented are limited to a
small subset of the total data. This paper is not intended to be a comprehensive review of the data; some of
the references listed above already serve that purpose. In addition, due to the sensitive nature of the data,
only generic information (pre-LCO, for example) will be provided instead of the parameters of the flight
condition such as Mach number and altitude.
3of20
AmericanInstituteofAeronauticsandAstronautics
Figure 3. Aileron inputs and resultant accelerations at the aft launcher location for the test run that encoun-
tered a mild flutter incident.
Figure 4. Contours of damping of AW1B mode as a function of Mach number and altitude.
4of20
AmericanInstituteofAeronauticsandAstronautics
III. Higher-Order Spectra
The primary benefit of higher-order spectra (HOS), also known as higher-order frequency response func-
tions, is that they provide information regarding the interaction of frequencies due to a nonlinear process.
For example, bispectra have been used in the study of grid-generated turbulence to identify the nonlinear
exchange of energy from one frequency to another (related to the turbulent cascade phenomenon). Lin-
ear concepts, by definition, cannot provide this type of information. In addition, higher-order spectra are
the frequency-domain version of the Volterra series. For details regarding this relationship, the reader is
referred to the recent article by Silva.17 Some very interesting and fundamental applications using the
frequency-domainVolterratheory18,19 andexperimentalapplicationsofVolterramethods20,21 areproviding
new “windows” on the world of nonlinear aeroelasticity.
In the recent work by Hajj and Silva,6,22 the aerodynamic and structural aspects of the flutter phe-
nomenon of a wind-tunnel model are determined via a frequency domain analysis based on a hierarchy of
spectral moments. The power spectrum is used to determine the distribution of power among the frequency
components in the pressure, strain and acceleration data. The cross-power spectrum, linear coherence, and
phase relation of the same frequency components between different signals are used to characterize the
bending and torsion characteristics of the model. The nonlinear aspects of the aerodynamic loading are
determined from estimates of higher-order spectral moments, namely, the auto- and cross-bispectrum.
For a discrete, stationary, real-valued, zero-mean process, the auto-bispectrum is estimated as23
M
Bˆ [l ,l ]= 1 XX(k)[l +l ]X∗(k)[l ]X∗(k)[l ] (1)
xxx 1 2 M T 1 2 T 1 T 2
k=1
where X(k)[l] is the Discrete Fourier Transform of the kth ensemble of the time series x(t) taken over a
T
time, T, and M is the number of these ensembles. The auto-bispectrum of a signal is a two-dimensional
function of frequency and is generally complex-valued. In averaging over many ensembles, the magnitude of
the auto-bispectrum will be determined by the presence (or absence) of a phase relationship among sets of
the frequency components at l , l , and l +l . If there is a random phase relationship among these three
1 2 1 2
components,theauto-bispectrumwillaveragetoaverysmallvalue. Shouldaphaserelationshipexistamong
these frequency components, the corresponding auto-bispectrum will have a large magnitude.24 Because a
quadraticnonlinearinteractionbetweentwofrequencycomponents,l andl ,yieldsaphaserelationbetween
1 2
them and their summed component, l +l , the auto-bispectrum can be used to detect a quadratic coupling
1 2
or interaction among different frequency components of a signal. The level of such coupling in a signal
can then be associated with a normalized quantity of the auto-bispectrum, called the auto-bicoherence, and
defined as
(cid:12) (cid:12)2
(cid:12)Bˆ [l ,l ](cid:12)
(cid:12) xxx 1 2 (cid:12)
b2 [l ,l ]= (2)
xxx 1 2 M(cid:12) (cid:12)2 M(cid:12) (cid:12)2
1 P(cid:12)X(k)[l ]X(k)[l ](cid:12) 1 P(cid:12)X(k)[l +l](cid:12)
M (cid:12) T 1 T 2 (cid:12) M (cid:12) T 1 2(cid:12)
k=1 k=1
By the Schwarz inequality, the value of b2 [l ,l ] varies between zero and one. If no phase relationship
xxx 1 2
exists among the frequency components at l , l , and l +l , the value of the auto-bicoherence will be at or
1 2 1 2
near zero (due to averaging effects). If a phase relationship does exist among the frequency components at
l , l , and l +l , then the value of the auto-bicoherence will be near unity. Values of the auto-bicoherence
1 2 1 2
between zero and one indicate partial quadratic coupling.
For systems where multiple signals are considered, detection of nonlinearities can be achieved by using
the cross-spectral moments. For two signals x(t) and y(t), their cross-bispectrum is estimated as
M
Bˆ [l ,l ]= 1 XY(k)[l +l ]X∗(k)[l ]X∗(k)[l ] (3)
yxx 1 2 M T 1 2 T 1 T 2
k=1
where X(k)[l] and Y(k)[l] are the Discrete Fourier Transforms of the kth ensemble of the time series x(t) and
T T
y(t), respectively, over a time, T. The cross-bispectrum provides a measure of the nonlinear relationship
amongst the frequency components at l and l in x(t) and their summed frequency component, l +l ,
1 2 1 2
in y(t). Similar to the auto-bispectrum, the cross-bispectrum of signals x(t) and y(t) is a two-dimensional
function in frequency and is generally complex-valued. In averaging over many ensembles, the magnitude
5of20
AmericanInstituteofAeronauticsandAstronautics
of the cross-bispectrum will also be determined by the presence, or absence, of a phase relationship among
sets of the frequency components at l , l , and l +l . If there is a random phase relationship among the
1 2 1 2
threecomponents, thecross-bispectrumwillaveragetoaverysmallvalue. Shouldaphaserelationshipexist
amongst these frequency components, the corresponding cross-bispectral value will have a large magnitude.
The cross-bispectrum is then able to detect nonlinear phase coupling among different frequency components
in two signals because of its phase-preserving effect.
Similarly to defining the auto-bicoherence, one can define a normalized cross-bispectrum to quantify the
levelofquadraticcouplingintwosignals. Thisnormalizedvalueiscalledthecross-bicoherenceandisdefined
as
(cid:12) (cid:12)2
(cid:12)Bˆ [l ,l ](cid:12)
(cid:12) yxx 1 2 (cid:12)
b2 [l ,l ]= (4)
yxx 1 2 M(cid:12) (cid:12)2 M(cid:12) (cid:12)2
1 P(cid:12)X(k)[l ]X(k)[l ](cid:12) 1 P(cid:12)Y(k)[l +l](cid:12)
M (cid:12) T 1 T 2 (cid:12) M (cid:12)T 1 2(cid:12)
k=1 k=1
If no phase relationship exists amongst the frequency components at l , l in x(t) and the frequency com-
1 2
ponent at l +l in y(t), the value of the cross-bicoherence will be at or near zero. If a phase relationship
1 2
does exist amongst these frequency components, the value of the cross-bicoherence will be near unity. Val-
ues of cross-bicoherence between zero and one indicate partial quadratic coupling. A digital procedure for
computing the auto and cross-bicoherence is given by Kim and Powers23 and is summarized by Hajj et al.25
A. Example
A simple example is now presented in order to explain some of the concepts associated with HOS. The
example26 consists of a linear (y) and a nonlinear (z) time series, both with added Gaussian noise (d). The
equations for these time series are
π π
x(k)=sin(2π(f )k+ )+sin(2π(f )k+ ) (5)
1 3 2 8
y(k)=x(k)+d(k) (6)
z(k)=x(k)+0.05x2(k)+d(k) (7)
where f = 0.12 Hz and f = 0.30 Hz.
1 2
Figures 5and 6eachcontainthetimeseriesandtheassociatedmagnitudeoftheFouriertransformsfor
the linear and nonlinear time series, respectively. The dominant frequencies (in this case 0.12 and 0.30 Hz)
areclearlyvisibleinbothfigures. Thefrequencycontentforthenonlineartimeseriesindicatestheexistence
of additional frequencies. The nature of these frequencies, whether or not these frequencies are random or
the result of a nonlinear coupling process, cannot be discerned from this analysis.
The power spectrum densities (PSD) for the linear and nonlinear time series are presented as figure 7.
Here again, the existence of additional frequencies in the PSD of the nonlinear time series (as compared to
thelineartimeseries)isobvious. However,thePSDinformationcannotbeusedtodiscerniftheseadditional
frequencies are random or the result of a nonlinear coupling process.
In order to determine if these additional frequencies are random or the result of a nonlinear coupling
process, the HOS for these time series must be computed. The magnitude of the bicoherence function for
the linear time series, y, is presented as figure 8. The two frequency axes correspond to the two frequency
indices in Eq. 2. If a significant peak is observed, this implies that the sum of those two frequencies is the
result of a quadratic (nonlinear) coupling. It can be seen that the magnitude of the bicoherence function for
the linear time series is quite low in value (compared to unity) and that there exist no dominant peaks that
would be indicative of a nonlinear coupling.
The magnitude of the bicoherence function for the nonlinear time series is presented as figure 9. In
contrast to the bicoherence function for the linear time series, the bicoherence function for the nonlinear
time series has several peaks at or close to unity, indicative of nonlinear (quadratic) coupling at the x- and
y-axis frequencies indicated.
Typically,inordertoenhancevisualizationandinterpretationofbicoherencefunctions,contourplotsare
viewed. The contour plot for figure 9 is presented as figure 10. The contour plot presents the frequencies
that have coupled quadratically to generate the new frequencies that were visible in the PSD functions.
6of20
AmericanInstituteofAeronauticsandAstronautics
Figure 5. Linear time series and magnitude of frequency response.
Figure 6. Nonlinear time series and magnitude of frequency response.
7of20
AmericanInstituteofAeronauticsandAstronautics
Figure 7. Power spectrum density (PSD) functions for linear and nonlinear time series.
Figure 8. Bicoherence function for the linear time series.
8of20
AmericanInstituteofAeronauticsandAstronautics
Figure 9. Bicoherence function for the nonlinear time series.
In addition, the symmetry associated with the computation of the bicoherence function is also evident in
the contour plot. The bicoherence function basically computes the correlation between two frequencies and
the sum of those two frequencies. Therefore, for example, the bicoherence function for the first frequency
(x-axis)being2Hz andthe secondfrequency(y-axis)being 5Hzwillbe the same asthebicoherence forthe
first frequency being 5 Hz and the second frequency being 2 Hz since their sum is the same.
This symmetry is presented in figure 11. As figure 11 indicates, knowledge of Region I (two such regions
in the first quadrant) and Region II (two such regions in the fourth quadrant) is sufficient to completely
define the remainder of the quadrants due to symmetry considerations. Regions I and II indicate that the
original frequencies have been quadratically-coupled, resulting in new frequencies that are the sum (positive
and negative) of the original frequencies. The new frequencies include 0.18, 0.24, and 0.42. The limitation
of Region I to the triangle shown is due to the fact that the summation of two frequencies cannot exceed
theNyquistfrequency(0.5, inthiscase). Therefore, thecombinationofthesecondfrequency(0.30Hz)with
itself is not included since 0.60 Hz would be greater than the Nyquist frequency.
For subsequent results, although the bicoherence function is computed for all frequencies (positive and
negative), Region I in the first quadrant will be of primary interest for the sake of simplicity. This will be
sufficient to demonstrate the applicability and value of HOS to flight flutter data.
IV. Results
Inthissection,dataanalysisusingtraditionalandHOSmethodsispresented. Thedataanalyzedconsists
of data acquired using sine dwells and lateral stick raps. The sine dwell data represents data acquired
subcritically (no immediate onset of flutter or LCO). The lateral stick rap data represents data that was
acquired as an LCO event was approached. The inputs (sine dwells and lateral stick raps) are in degrees
of aileron motion and the outputs are accelerations at the forward launcher location (at the wing tip). In
addition, computed results using the lateral stick rap data, generated using a method developed by Silva,5
arepresented. Usingthismethod,animpulseresponseisidentifiedfromtheexperimentalinput/outputdata
and used to determine the level of linear/nonlinear behavior at a given condition.
9of20
AmericanInstituteofAeronauticsandAstronautics
Figure 10. Contour plot of the bicoherence function for the nonlinear time series.
Figure 11. Contour plot of the bicoherence function for the nonlinear time series with symmetry regions
identified.
10of20
AmericanInstituteofAeronauticsandAstronautics
