An Overview of Flutter Prediction in Tests Based on Stability Criteria in DiscreteTime Domain
 Author: Matsuzaki Yuji
 Organization: Matsuzaki Yuji
 Publish: International Journal Aeronautical and Space Sciences Volume 12, Issue4, p305~317, 30 Dec 2011

ABSTRACT
This paper presents an overview on flutter boundary prediction in tests which is principally based on a system stability measure, named Jury’s stability criterion, defined in the discretetime domain, accompanied with the use of autoregressive movingaverage (ARMA) representation of a sampled sequence of wing responses excited by continuous air turbulences. Stability parameters applicable to two, three and multimode systems, that is, the flutter margin for discretetime systems derived from Jury’s criterion are also described. Actual applications of these measures to flutter tests performed in subsonic, transonic and supersonic wind tunnels, not only stationary flutter tests but also a nonstationary one in which the dynamic pressure increased in a fixed rate, are presented. An extension of the concept of nonstationary process approach to an analysis of flutter prediction of a morphing wing for which the instability takes place during the process of structural morphing will also be mentioned. Another extension of analytical approach to a multimode aeroelastic system is presented, too. Comparisons between the prediction based on the digital techniques mentioned above and the traditional damping method are given. A future possible application of the system stability approach to flight test will be finally discussed.

KEYWORD
Flutter test , Flutter boundary prediction , Stability criterion , Discretetime domain

1. Introduction
Flutter boundary prediction during wind tunnel and flight tests is one of the most important, but very difficult tasks imposed in the process of development of a new airplane (Bisplinghoff, et al., 1955; Fung, 1955). Flutter boundary prediction methods in which only measured wing response signals are analyzed can be classified into two categories:
A. Damping approach based on measurement of the damping coefficient of flutter mode,
B. System stability approach based on the estimation of stability of an aeroelastic system.
The category B is further classified into two subcategories:
B1) Flutter margin for continuoustime systems (FMCS) (Poirel et al., 2005; Price and Lee, 1993; Zimmerman and Weissenburger, 1964) based on Routh’s stability test,
B2) Jury’s stability criterion defined in the discretetime domain (Matsuzaki and Ando, 1981; Matsuzaki and Torii, 1990), and the flutter margin for discretetime systems (FMDS) (Bae et al., 2005; Torii and Matsuzaki, 1999, 2001) founded on the Jury criterion.
The first category A is very traditional, and it is still very popular in practical applications in spite of its intrinsic difficulties. As clearly seen in Fig. 1 of an experimental study (Irwin and Guyett, 1968) published late in the 1960’s, the damping coefficient of the flutter mode often reveals two problems in a critical speedrange very close to the flutter boundary:
a) Rapid change in its value, that is, suddenly decreasing and becoming zero, as seen both in analyses and experiments of the socalled explosive flutter,
b) Large scattering of its measured values in experiments.
As for the second difficulty b), a morerecent survey (Koenig, 1995) on the results of 16 different flutter tests pointed out that there was a scatter of about 30% in the damping coefficient while only 3% in the corresponding frequency.
To compensate the first difficulty a), a criterion for a twomode system which was characterized by a gradual change, named flutter margin (FMCS2), was proposed (Zimmerman and Weissenburger, 1964). FMCS2 is based on Routh’s stability test in the continuoustime domain, and is expressed with measured frequencies and damping coefficients of the first two modes. Therefore, FMCS2 also spreads widely due to the second problem b) mentioned above. To account for uncertainty caused by the scattering, two stochastic approaches (Poirel et al., 2005) were introduced in estimating FMCS2. The stochastic analysis, However, the stochastic analysis cannot improve the accuracy itself in the prediction of the critical boundary. An extended criterion for a threemode flutter (FMCS3) was also proposed (Price and Lee, 1993).
Hereafter, we will focus on Subcategory B2, that is, the flutter boundary prediction based on the stability criteria defined in the discretetime domain, because digital techniques have enormously developed for the last two decades. Therefore, we may now apply easily the modern estimation and identification methods to data analysis of measured responses to estimate flutter characteristics in a subcritical range. In other words, it is possible that the flutter margin, modal frequencies and damping coefficients are estimated accurately from certain duration of the response even contaminated with inevitable noises, including mechanical and electric ones.
The present overview is a revised version of the author’s paper (Matsuzaki, 2011b), which was updated from the original (Matsuzaki, 2011a) by referring to a couple of recent publications and including the Section 5.
2. DiscreteTime Approach
In the discretetime domain, the advanced counterpart of RouthHurwitz’s stability criterion (Fung, 1955) in the continuoustime case is Jury’s stability criterion being developed from the one first established by SchurCohn (Jury, 1964). Like the RouthHurwitz condition, Jury’s criterion is defined in terms of the coefficients of the characteristic equation of the system, as will be described later. Needless to say, once characteristic coefficients are estimated from sampled measured data, not only stability but also the dynamic characteristics of the system, such as the modal frequencies and damping coefficients, can be evaluated at the same time.
Early in 1980’s, the present author and his coworker (Matsuzaki and Ando, 1981) proposed an innovative approach for predicting the flutter boundary based on Jury’s stability criterion together with an autoregressive movingaverage (ARMA) representation for sampled random responses. Applying it to a flutter test of a cantilever wing model performed in a transonic wind tunnel, they showed that the Jury criterion was much more eff ective in predicting the flutter boundary than the damping coefficient of the flutter mode. This criterion was also successively applied to the prediction of divergence boundary of a forwardswept wing tested in a supersonic wind tunnel (Matsuzaki and Ando, 1984). Two decades later, a new stability parameter for a twomode flutter, named "flutter margin for discretetime systems" (FMDS2) (Torii and Matsuzaki, 1999, 2001), was derived from the Jury criterion. The new parameter is mathematically equivalent to FMCS2 in the continuoustime case, provided that the sampling period is suffi ciently short (Torii, 2002; Torii and Matsuzaki, 2001). FMDS2 was recently extended to FMDS3 for a threemode flutter (Torii and Matsuzaki, 2011). An extension of FMDS2 to a multimode flutter (FMDSN) (Bae et al., 2005) was also proposed to apply with success to a wind tunnel test of a wing with a fl ap. However, FMDSN in the case of the twomode does not agree with the original FMDS2, as pointed out in Matsuzaki (2005). Recent numerical studies (Hallissy and Cesnik, 2011; McNamara and Friedmann, 2007) showed that FMDSN did not work well as a predictor of flutter onset for multimode wing systems.
2.1 Discretetime series representation
The ARMA process for a sampled signal of random responses, Jury’s stability criterion and three versions of FMDS, i.e., FMDS2, FMDS3, and FMDSN, will be briefly summarized.
It is assumed that the stationary response,
y(t) , of a wing flying at a fixed speed in an atmosphere with a Gaussiantype turbulence of a given variance, is sampled at a constant interval T to produce a sequence of random variables, {y(kT) , or simply written,y(k); k =1, 2, …,N }. The stationary random sequence is assumed by an ARMA model (Goodwin and Sin, 1984) aswhere
α _{0}=1 andα _{2M}≠0. The order2M , AR coefficients {α m} and MA coefficients {β m} withm =1, …, 2M, are an unknown integer and real numbers to be estimated from the sampled data, {y(k) }. A sequence, {e(k); k =1, …,N 1} represents noise inputs to the system, that is, a Gaussian independent random sequence with zero mean and an unknown variance σ^{2}. Application of the ztransform (Jury, 1964) to the lefthand side of Eq. (1) yields the characteristic equation:where
and
z _{m}^{*}(m = 1, 2, ... ,M ) stands for a complex conjugate ofz _{m}. It is obvious that M represents the number of the modes consisting of the aeroelastic system. The eigenfrequencies, ω_{i} , and the damping coefficients, η_{i} , of the ith mode are respectively given asJury’s criterion (Jury, 1964) is defined by
where
X_{i} andY_{i} arei x i determinants whose elements are one of the AR coefficients or zero. The dynamic and the static instability, namely, flutter and divergence of the system withM degrees of freedom is respectively given byThe flutter margin for a two, three, or multimode system (FMDS2 [Torii and Matsuzaki, 2001], FMDS3 [Torii and Matsuzaki, 2011], or FMDSN [Bae et al., 2005]) is respectively expressed as
From Eqs. (9) and (10), we have a unified formula (FMDSM for M = 2 and 3):
It is proved (Torii, 2002; Torii and Matsuzaki, 2001) by using the Tustin transform that FMDS2 is mathematically equivalent to FMCS2 in the continuoustime domain, provided that the sampling period T is suffi ciently short. FMDSN defined by Eq. (11) in the case of M = 2 does not agree with FMDS2 given by Eq. (9).
3. Applications to Flutter Test Data
3.1 Stationary process test
The first application of Jury’s stability criterion to a flutter test was made by using a cantilever wing model, as shown in Fig. 1a, tested at Mach number 1.17 in a transonic blowdown type wind tunnel (Matsuzaki and Ando, 1981). During the test, at each blow the dynamic pressure was kept constant for about ten seconds and to measure the response of the wing randomly excited by natural flow turbulences. The response at each blow was stationary, because its mean value and standard deviation remained essentially the same with time, more accurately its covariance was the same. The measurement was repeated by increasing the dynamic pressure stepwise. As the flutter criterion, the parameter
F ^{}(3)=0, obtained by settingM =2 in Eq. (7), was applied to measured data which had been bandpassed to contain contributions mainly from the first two modes. In Fig. 1b, evaluated values ofF ?(3) are shown by solid circles against the dynamic pressure,Q , between about 0.5 and 0.9 kg/cm^{2}. The flutter boundary was predicted to be aboutQ = 0.95 kg/ cm^{2} by fitting a straight line to the solid circles given only between about 0.5 and 0.8 kg/cm^{2}, and by extrapolating the fi tted line to cross a point on the horizontal coordinate, which was close to the experimental value (Q_{F} = 0.97 kg/cm^{2}) given by a cross on the coordinate. In Figs. 1c and d, the modal frequencies,f _{1} andf _{2}, and the damping coefficients of the first two modes, η_{1} and η_{2}, calculated respectively by using Eqs. (4) and (5) are plotted by circles and triangles atthe same pressures as for
F ^{}(3) shown in Fig. 1b. It is clearly very difficult to predict the boundary simply by using the damping coefficients, as seen in Fig. 1d.Figures 2be show numerical results obtained by applying the approach based on FMDSM to experimental data of a wing model shown in Fig. 2a, which were measured at M = 2.51 at a supersonic wind tunnel. FMDS2 and FMDS3 were, respectively, applied to sampled sequences obtained by bandpassing the measured data so as to contain mainly the responses due to the first two modes (Torii and Matsuzaki, 2001) and due to the first three (Torii and Matsuzaki, 2011). In Figs. 2b and d, estimated values of FMDS2 and FMDS3 are plotted by solid circles in the pressure range up to
q = 100 kPa, with the measured flutter boundary,q_{F} = 113.5 kPa, shown by a cross on the horizontal coordinate. A dotted and a solid straight line are fitted to the circles shown belowq = 100 kPa and extended to cross the horizontal coordinate at a point representing an estimated boundary,located very close to the cross on the coordinate. The calculated damping coefficients, η_{1} and η_{2}, are also given up to
q = 100 kPa in Fig. 2c and η_{1}, η_{2} and η_{3} in Fig. 2e. In these figures, again, it is not possible to predict the boundary shown by the cross by using the estimated damping coefficients. It is clear that FMDSM for M = 2 and 3 is effective in predicting the flutter boundaryfrom a low dynamic pressure range. The effectiveness of the flutter margins was also analytically confirmed by using the socalled twodimensional simple wing without and with a flap, that is, with two and three degrees of freedom, respectively.
As for a fourmode system (M = 4) of a cantilevered wing with a flap (Bae et al., 2005), the extended parameter
F_{N} defined by Eq. (11) was used to predict its flutter boundary in a test. Figure 3 shows a comparison between the application ofF_{N} during the test and the flutter analysis by the Ug method.F_{N} decreases linearly to find the boundary at 832 Pa. In the Ug diagram, there are two curves of g which cross the horizontal coordinate. With a small amount of structural damping, for example,g = 0.02 being taken into account in this figure, the critical point given by the fourth mode at 717 Pa would shift at about 1,000 Pa. However, the other crossing point at about 800 Pa due to the second mode would move to around 830 Pa to give the flutter boundary. Hence, agreement among the Ug approach, the prediction byF_{N} and the experimental value of the boundary was good.3.2 Nonstationary process test
To determine the flutter boundary in a single test, the dynamic pressure was often increased at a fixed rate until it arrived at its critical value (Matsuzaki and Ando, 1985; Torii and Matsuzaki, 1997, 2002). Figures 4a and b illustrate time histories of the dynamic pressure which increased from 76 to
117 kPa and the signal of a strain gage glued on the surface near the root of the same wing as used for the stationary test and shown in Fig. 2a. The standard deviation of the response signal increased gradually with time, and the wing suddenly started to flutter at
q_{F} = 113.5 kPa (Torii and Matsuzaki, 2002), as mentioned in Section 3.1. For the prediction of the flutter boundary, sequential data sampled from the strain signal shown in Fig. 4b were assumed in a nonstationary ARMA process, and its coefficients were evaluated through theuse of recursive maximum likelihood estimation technique (RMLE) (Ljung, 1987), since they were timedependent. In Fig. 4c, recursively estimated values of FMDS2 are plotted by a wavy solid curve against the dynamic pressure up to the flutter point given by a solid square on the horizontal coordinate. The fluctuation of the estimated values was sufficiently small except for a low pressure range, say, below 85 kPa, where the estimation was not well converged around the true values because RMLE had started with a set of initial values arbitrarily selected. In the nonstationary case, FMDS 2 was again almost linearly decreased over the range up to
q_{F} = 113.5 kPa except below 85 kPa. For comparison, numerical values of FMDS2 of the stationary test, like those shown in Fig. 2b, are also given by solid circles, which lie mostly on the curve. Before the wing model was broken during the nonstationary test, the stationary test was first performed to measure the stationary responses4. Numerical Simulations
To exploit and confirm analytically the feasibility of the system stability method, the author and his coresearcher performed numerical simulations by using the twodimensional simple wing with only two degrees of freedom in an incompressible flow, or by introducing a panel system supported with a distributed spring and exposed to a supersonic flow. In these studies, the generalized theory of unsteady aerodynamic forces (Edwards, 1977) was used to generate responses of the wing and the panel in an arbitrary motion caused by random air turbulences.
4.1 Nonstationary process analysis
Like in Section 3.2, nonstationary responses of the twodimensional simple wing, shown in Fig. 5a, exposed to a turbulent flow whose dynamic pressure increased at a fixed rate, were represented by the ARMA process, and its coefficients were evaluated by using RMLE. The evaluation of the response was repeated fifty times by using a different sequence of flow turbulence. Figure 5b shows the mean and standard deviation of evaluated values of FMDS2, respectively, by a solid line and a small Ishaped column against the dimensionless pressure,
q/q_{F} , from 0.8 to 1.02 (Torii and Matsuzaki, 2002), although the evaluation was made fromq/q_{F} = 0, whereq_{F} stands for the flutter boundary. A dotted line, named "True value", represents FMDS2 of a flutter analysis, being redrawn from a part of the curve shown, for instance, in Fig. 4 of Torii and Matsuzaki (2001), where FMDS2 decreased in a linear manner fromq/q_{F} = 0 to 1.0. FMDS2 may, therefore, provide a good tool to predict theflutter boundary for the nonstationary process, too, because its mean value changed nearly in a linear manner. It is noted that, because of a dynamic effect due to the nonstationary process, the flutter boundary shifted to 1.02 from 1.0 of the true value obtained in the flutter analysis, as shown in Fig. 5b.
4.2 Morphing aircraft
The recent development of morphing aircraft technology has now presented a newtype of technical problems. Among them, aeroelastic instability is one of the most challenging issues (Lee and Weisshaar, 2005; Tang and Dowell, 2008). In an ordinary airplane the flight speed is the most visible variable to the instability. In other words, the higher the flight speed is, the more critical the stability becomes. On the other hand, as for the morphing wing, we may consider that, even though the flight speed is kept fixed, the wing will onset flutter in the process of morphing in which the structural adaptation happens to lead to a "coalescence" of the frequencies of aeroelastic modes of the wing. In other words, the morphing technology presents a new type of nonstationary process. Accurate measurement of modal damping coefficients is considered to be harder in such nonstationary situations. Applicability of the nonstationary approach, mentioned in Sections 3.2 and 4.1, to a simple morphing wing was studied (Matsuzaki and Torii, 2005, 2006, 2007) by assuming that the stiffness of the wing shown in Fig. 5a was representatively given by the natural frequencies of the heaving and the pitching motion, ω
_{gh} and ω_{gα} , and that the morphing mechanism was influential only on the stiffness, namely, the two natural frequencies of the wing. In a numerical simulation of morphing, only ω_{gh} increased continuously at a constant rate with ω_{gα} being fixed.In Fig. 6a, numerical results of the flutter analysis, that is, FMDS2 and the damping coefficient of the critical mode, η
_{2} , are plotted against ω_{gh} from 50.0 to 80.0 rad/sec, respectively by a solid and a chained line, while ω_{gα} was kept at 100 rad/ sec. As ω_{gh} increased from ω_{gh} = 50 rad/sec, FMDS decreased uniformly and became zero at ω_{gh.cr} = 76.2 rad/sec. On the other hand, η_{2} remained almost unchanged up to about ω_{gh} = 69 rad/sec, and then started to decrease rapidly to zero at ω_{gh.cr} . Figure 6b illustrates numerical results of the nonstationary simulation analysis in the morphing process, in which ω_{gh} increased at a constant rate of 0.5 rad/sec per second. The coefficients of the nonstationary ARMA process were again estimated by RMLE. FMDS2 decreased almost linearly in a range from above ω_{gh} = 55 rad/sec and became zero at ω_{gh} = 76.5 rad/sec, which was higher than ω_{gh.cr} due to the effect of the dynamic process. It was confirmed (Matsuzaki and Torii, 2012) that FMDS2 was useful in predicting the flutter boundary, ω_{gh.cr} , of the morphing wing from a sufficiently lower range of ω_{gh} .4.3 Multimode flutter system
Airplane wings with auxiliary airfoils or additional attachments encounter often aeroelastic instabilities involved with a higherorder mode. A very flexible wing of a highaspect ratio became critical also due to an involvement of a highorder mode (Hallissy and Cesnik, 2011). In order to develop systematically the flutter boundary prediction applicable to a multimode flutter, numerical simulation analyses (Matsuzaki, 2009, 2010) were performed by utilizing responses of a panel which was supported with a distributed spring and exposed to a supersonic flow, as shown in Fig. 7a. As is wellknown, the onset of panel flutter does not usually lead to an immediate destruction of the panel structure, so that the prediction of the flutter boundary is less serious than
in wing flutter. The reason for using the panel system is that we may easily create an aeroelastic situation in which the two modes intentionally designated are coupled to induce flutter, and that the twodimensional quasisteady supersonic flow approximation is applicable as a generalized unsteady theory (Edwards, 1977) to generate the panel response caused by air turbulence. Furthermore, the responses of the panel obtained in a subcritical speed resembled very much those of a wing model often observed during a wind tunnel test, as will be seen later. For this multimode system, it is assumed that the value of an effective spring constant of the ith mode, k
_{i} , is selectable arbitrarily, so that proper adjustment of the spring constants makes it possible to cause flutter due to an aeroelastic coupling between the two designated modes. Asthe first step, the fourmode system was analyzed, in which four different couplings between (Case I) the first and the second mode, (II) the second and the third, (III) the third and the fourth and (IV) the first and the fourth were examined. Results of flutter analysis on Case IV, that is,
F ^{？}(7)?X_{7} ?Y_{7} , the modal frequencies and the damping coefficients of the four modes are shown against the dynamic pressure, Q, in Figs. 7bd, respectively. The coupling between the first and the fourth mode is clearly seen in Fig. 7c. Jury’s parameter F^{}(7) decreased more or less in a uniform manner to become zero at the flutter boundary,Q_{F} = 4.9× 10^{?4} kg/mm^{2}. On the other hand, as shown in Fig. 7d, the damping coefficient of the first mode shown by a solid line kept its value very close to 0.01 up to aboutQ = 4.8× 10^{?4}kg/mm^{2} and with a slight increase inQ it quickly became zero at
Q_{F} .Based on this preliminary analysis, numerical simulations on flutter boundary prediction of the aeroelastic system were performed by assuming that the dynamic pressure increased in a stepwise manner. Figures 8 show time histories of panel response (Case IV) of one second at four dynamic pressures, that is, a)
Q = 1.0× 10^{?4}, b) Q = 3.0× 10^{?4}, c) Q = 4.0× 10^{?4}, d) Q = 5.0× 10^{?4} (kg/mm^{2}). The responses shown in these figures are quite similar to those of wing models which we often observe during wind tunnel tests. As the critical pressure wasQ_{F} = 4.9× 10^{?4} kg/mm^{2}, the response diverged in an oscillatory manner in Fig. 8d. At ten dynamic pressures including these three Cases a) to c),F ^{?}(7) was evaluated. In Fig. 9, the mean and the standard deviation ofF ^{?}(7) based on ten estimations at each dynamic pressure are given by a square and a pair of triangles. It is obvious that a simple extrapolation of a solid line connecting the squares will cross the horizontal coordinate at a point close to the boundary,Q_{F} , indicatedby an arrow. As for the remaining three cases, similar results were obtained.
5. Discussion
Before discussing the matters related to the system stability approach, we note that the first application of the timeseries representation of wing response was given by using the AR process in 1978 (Onoda, 1978), although the system stability was not taken into account in the analysis. A number of researchers used also the AR or the ARMA representation of wing responses in the 1980s and 90s, for example, Batill et al. (1992), Pak and Friedmann (1992), and Roy et al. (1986).
As discussed above, the flutter boundary prediction based on the system stability approach has great advantages due to advanced digital techniques, although the number of its applications is very small so far. For the present, it may serve as a complement to the damping method. This is because the latter approach has a very long history of practical applications and experiences, whereas the former needs to experience much more actual applications. Modern estimation methods supposedly allow us to eliminate much of noise effects from measured noisecontaminated data, and estimate sufficiently accurate system parameters. For this end, however, it is necessary to appropriately select initial input parameters for a computing program, such as the sampling period, T, the number of data, N, the upper and the lower limit of a band pass filter, f_{U} an f_{L}, etc., as well as specific parameters associated with the programs used. Some comparisons among estimated results of FMDS, modal frequencies, etc., obtained by changing input parameter values, were given, for instance, in Torii and Matsuzaki (1992, 1995, 2001). For the parameter calibration, the modal frequencies and damping coefficients directly measured in a test will also serve as useful references. The evaluation of Jury’s criterion, FMDSM and FMDSN is rather simple and easy, once the order 2M and the AR coefficients of the system have been estimated. In an actual application, it is recommended to evaluate all these criteria as well as the modal frequencies and damping coefficients given respectively by Eqs. (4) and (5).
Although Jury’s criterion, more specifically,
F ^{?}(2M ?1) given by Eq. (7), is applicable to an aeroelastic system with any degrees of freedom, an effort for extension of the effective FMDS3 to a flutter margin for highermode systems is important, because multiple frequencies of uncoupled modes may exist very close to or between the frequencies of the coupled modes. In such a case, it is impossible to filter the sequence of measured data into one which contains only three modes including the coupled two.As is wellknown, the damping approach is not so effective, particularly in the case of the socalled explosive flutter for which the critical damping coefficient starts to decrease quickly in a narrow range close to the boundary and becomes zero to induce a violent selfexcited oscillation. An expert of an American aircraftmanufacturing company says that a half of dozen wing models are constructed for a wind tunnel test because they are often destructed by the onset of explosive flutter. Needless to say, direct measurement of the damping coefficients from a oncejerked and quicklydamped oscillation is easily influenced by air turbulence because turbulent flows change their characteristics from one instant to another. On the other hand, by using modern digital techniques which developed remarkably over the last couple of decades, we may analyze stationary and nonstationary random responses measured for certain period, and extract more reliable values of the system parameters from measured data.
As for a wind tunnel test, from our past experiences although very limited, it seems that an ordinary wind tunnel provides sufficiently effective turbulences for the system stability approach. If the flow is calm, then one may install a vortex generator in an upstream of the test section to create required flow turbulences. On the other hand, as for a flight test in an open atmosphere, no one can expect any turbulence which will continue with the same strength for a sufficiently long range and a long period. However, recent development of unmanned aircraft technology is so remarkable that we may think of introducing a turbulencegenerating system composing of a unmanned aircraft(s) equipped with a vortex generator(s), and also controlling the system to maneuver at a fixed distance ahead of a tested airplane and to produce air turbulence satisfying requirements in a region where the tested airplane flies. Figure 10 illustrates a schematic picture of flight test in which a vortexgenerating system and a tested airplane fly together at a supersonic speed.
6. Concluding Remarks
This paper presented an overview on the flutter boundary predictions based on the system stability approach for digitalized systems, that is, Jury’s stability criterion and flutter margin for discretetime systems (FMDSM and FMDSN). The applications of the prediction method to a flutter system with multimodes and to the nonstationary flutter test in which the dynamic pressure uniformly increased were mentioned. The extension of the nonstationary approach used for ordinarytype wings to the prediction of flutter of a morphing wing was reviewed, too.
Difficulties encountered by the traditional damping method and the system stability approach were discussed. Regarding a future application of the system stability method to flight test in an open air, a conceptual idea of a turbulencegenerating system composed of an unmanned aircraft(s) which flies at a fixed distance ahead of a tested airplane was proposed.

[Fig. 1.] Flutter boundary prediction of a wing tested in a transonic windtunnel: (a) wing model, (b) Estimated F?(3) vs. QI, (c) estimated modal frequencies vs. Q, (d) estimated damping coefficients vs. Q (a cross on the coordinate: measured boundary, QF = 0.97 kg/cm2). From Matsuzaki and Ando (1981). Reprinted with permission of AIAA.

[Fig. 2.] Flutter boundary prediction of a wing tested in a supersonic windtunnel: (a) wing model, (b) estimated FMDS2 vs. q, (c) estimated two damping coefficients vs. q. From Torii and Matsuzaki (2001); (d) estimated FMDS3 vs. q, (e) estimated three damping coefficients vs. q. From Torii and Matsuzaki (2011); (a cross on the coordinate: measured flutter boundary, qF= 113.5 kPa). Reprinted with permission of AIAA.

[Fig. 3.] Analysis and experiment on flutter of a cantilever wing with a flap: comparison between the Ug method and FN. From Bae et al. (2005). Reprinted with permission of AIAA.

[Fig. 4.] Nonstationary process test and flutter boundary prediction : (a) dynamic pressure increasing at a contant rate, (b) strain response of the wing shown in Figs. 2a and (c) flutter margin for discretetime system (FMD2) evaluated by recursive maxinum likelihood estimation (a solid circle: FMDS2 evaluated for stationary date;a solid square on the coordinate: qF = 113.5 kPa). From Torii and Matsuzaki (2002)

[Fig. 5.] Flutter boundary prediction of a nonstationary process: (a) wing model, (b) mean values and standard deviations of flutter margin for discretetime system, FMDS2 vs. dimensionless dynamic pressure (a dotted line: flutter analysis). From Torii and Matsuzaki (2002).

[Fig. 6.] Flutter analysis and numerical simulation during structural morphing: (a) flutter margin for discretetime system (FMDS 2) and 2nd mode damping vs. natural frequency, ωgh, of the heaving mode. From Matsuzaki and Torii (2006), (b) estimated FMDS2 by recursive maximum likelihood estimation. From Matsuzaki and Torii (2007). Reprinted with permission of AIAA.

[Fig. 7.] Flutter boundary prediction of a fourmode system: (a) twodimensional supersonic panel flutter, (b) Jury’s criterion, F?(7)=X7?Y7, (c) modal frequencies, (d) modal damping coefficients. From Matsuzaki (2010). Reprinted with permission of AIAA.

[Fig. 8.] Wing responses and their power spectral densities for fourmode system: (a) Q = 1.0× 10?4, b) Q = 3.0× 10?4, c) Q = 4.0× 10?4, d) Q = 5.0× 10?4 kg/mm2. From Matsuzaki (2010). Reprinted with permission of AIAA.

[Fig. 9.] Mean values and standard deviations of estimated F?(7) = X7?Y7: QF = 4.9× 10?4 kg/mm2. From Matsuzaki (2010). Reprinted with permission of AIAA.

[Fig. 10.] Flight flutter test in a supersonic range, using a turbulencegenerating system: in side of a region which two Mach cones overlap, artificial turbulence satisfying requirements for flight test is developed and the tested airplane flies.