Higher-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