A flutter boundary prediction method for the three-dof wing

Yang Li1 , Li Zhou2 , Bingcai Yang3

1, 2, 3State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China

2Corresponding author

Journal of Vibroengineering, Vol. 17, Issue 8, 2015, p. 4507-4516.
Received 3 April 2015; received in revised form 15 June 2015; accepted 5 July 2015; published 30 December 2015

Copyright © 2015 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

Through curve fitting, the so-called flutter-margin method can estimate a good flutter boundary, however, is not suitable for the binary wing with a control surface. This paper develops the flutter-margin method and proposes a new stability criterion based on the dynamic equation of the three-dof wing and the Routh’s stability theory. The analysis of the trend and error-resistant ability show that: the new criterion, behaving a steady decrease to zero as instability is approached, varies in a sensibly linear manner with the velocity squared and could obtain good estimates of the flutter speed from data at a low speed, thus overcoming the defects of the velocity-damping method effectively. In addition, the criterion is more dependent upon the frequencies than on the decay rates, and is presented to be a good linear trend even when there is a 50 % damping maximum allowable error and 5 % frequency error.

Keywords: flutter margin method, two-dimensional wing, boundary estimate, Routh’s stability criterion, modal parameters.

1. Introduction

Flutter of aircraft is often a disastrous phenomenon, and thus it is essential to obtain the flutter boundary through wind tunnel and flight test [1]. But, considering the safety of life and property, it is difficult to reach the critical value of flutter. Traditionally, there is a flutter boundary prediction method so called velocity-damping method that fit and extrapolate the “trend” of the modal damping. When damping trend decreases to be zero, the airspeed is thought to be the flutter boundary [2-4]. Unfortunately, sometimes the modal aeroelastic damping of the nature flutter, especially so-called explosive flutter, will increase continually until just before instability, and then reduce suddenly. Furthermore, it is unknown which mode will go unstable, and certainly there will be some scatter in the modal damping measurements, thus lead to uncertainty about whether or not flutter will occur at the next increment of airspeed.

One method so-called “flutter-margin method” [5] is proposed to reduce this uncertainty. In this method, Zimmerman and Weissen burger make use of both the modal frequency and damping information and show that the flutter margin varies in a gradual manner with dynamic pressure. In particular when the instability is approached, the flutter-margin will decrease steadily to zero. Flutter experiments of a T-tail, a wing and a stabilizer have shown that the margin is applicable to two modes situation [5]. And more investigations of the use of the method have been done by Bennett [6] and Jennifer [7] on wind-tunnel aero elastic data and by Katz [8] and Lee [9] on flight test data.

However the margin can be used to binary instabilities only, then Price and Lee [10] developed a flutter margin method for three modes. Through the results of simulated experimental data, they proposed that the new margin will vary in a sensibly linear manner with dynamic pressure, thus could be a good tool to predict the flight flutter. But there is actually a mistake in their deducing process, and the simplified analysis of the margin is lack of rigorous.

In this paper, the flutter-margin method will be developed to the trinary flutter situation, and a new stability criterion will be proposed based on the dynamic equation of the three-dof wing and the Routh’s stability theory. The variation of this new criterion with airspeed will be described upon the simulated experimental data and the influence of modal parameter errors will be discussed latter.

2. The dynamic equation of the three-dof wing

Give a binary wing model with a control surface as shown in Fig. 1. The wing has three degrees of freedom: plunging motion h, pitch motion α and control surface yaw motion β. The coordinate origin is set to be the center of the airfoil chord and the chord length is 2b. In Fig. 1, kα, kh and kβ represent the stiffness coefficient of each dof.

Fig. 1. A binary wing model with a control surface

 A binary wing model with a control surface

We can obtain the motion differential equation of the wing as follows:

(1)
m h ¨ + m x α α ¨ + m x β β ¨ + k h h - F = 0 , m x α h ¨ + m r α 2 α ¨ + m r β 2 + m x β c - b - a - b     β ¨ + k α α - T α = 0 , m x β h ¨ + m r β 2 + m x β c - b - a - b     α ¨ + m r β 2 β ¨ + k β β - T β = 0 ,

where F is aerodynamic force, Tα is the aerodynamic moment for elastic axis, Tβ is the aerodynamic moment for hinge coordinate of the control surface.

Depend on Theodorsen’s theory [11], the quasi-steady aerodynamic force on the three-dof binary wing can be:

(2)
L = π ρ a b 2 h ¨ + V α ˙ - b a - α ¨ - V π T 4 β ˙ - b π T 1 β ¨ +   2 π ρ a V b Q a ,
(3)
T α = π ρ a b 2 b a - h ¨ - V b 1 2 - a - α ˙ - b 2 1 8 + a - 2 α ¨ - V 2 π T 4 + T 10
            + V b π - T 1 + T 8 + c - - a - T 4 - 1 2 T 11 β ˙ + b 2 π T 7 + c - - a - T 1 β ¨
            + 2 π ρ a V b 2 a - + 1 2 Q a ,
(4)
T β = π ρ a b 2 b π T 1 h ¨ + V b π 2 T 9 + T 1 - a - - 1 2 T 4 α ˙ - 2 b 2 π T 13 α ¨
          - V π 2 T 5 - T 4 T 10 β + V b 2 π 2 T 4 T 11 β ˙ + b π 2 T 3 β ¨ - ρ a V b 2 T 12 Q a ,  

where T1-T13 are configuration constants.

Substituting Eq. (2)-(4) into Eq. (1), we can obtain:

(5)
a 1 h ¨ + V a 2 h ˙ + a 3 h + b 1 α ¨ + V b 2 α ˙ + V 2 b 3 α + c 1 β ¨ + V c 2 β ˙ + V 2 c 3 β = 0 ,
(6)
a 4 h ¨ + V a 5 h ˙ + b 4 α ¨ + V b 5 α ˙ + V 2 b 6 α + b 7 α + c 4 β ¨ + V c 5 β ˙ + V 2 c 6 β = 0 ,
(7)
a 6 h ¨ + V a 7 h ˙ + b 8 α ¨ + V b 9 α ˙ + V 2 b 10 α + c 7 β ¨ + V c 8 β ˙ + V 2 c 9 β + c 10 β = 0 ,

where a1-a7, b1-b10 and c1-c10 are configuration constants in derivation process.

Writing these equations in operational form with P representing d/dt and performing some algebraic manipulations, they can be written as:

(8)
a 1 P 2 + V a 2 P + a 3 h + b 1 P 2 + V b 2 P + V 2 b 3 α + c 1 P 2 + V c 2 P + V 2 c 3 β = 0 ,
(9)
a 4 P 2 + V a 5 P h + b 4 P 2 + V b 5 P + V 2 b 6 + b 7 α + c 4 P 2 + V c 5 P + V 2 c 6 β = 0 ,
(10)
a 6 P 2 + V a 7 P h + b 8 P 2 + V b 9 P + V 2 b 10 α + c 7 P 2 + V c 8 P + V 2 c 9 + c 10 β = 0 .

These three simultaneous differential equations in the unknown h, α and β can be reduced to a single differential equation in each of the unknowns using straightforward operational techniques. If the unknown h and β are eliminated, we can get the differential equation for α:

(11)
P 6 + A 5 P 5 + A 4 P 4 + A 3 P 3 + A 2 P 2 + A 1 P + A 0 α = 0 .

The exact expressions for the coefficients A0-A5 are unimportant for our purpose, only the form or compositions are pertinent. The coefficients A0-A5 take the form:

(12)
A 5 = K 5 V , A 4 = K 41 V 2 + K 42 , A 3 = K 31 V 3 + K 32 V , A 2 = K 21 V 4 + K 22 V 2 + K 23 , A 1 = K 11 V 5 + K 12 V 3 + K 13 V , A 0 = K 01 V 6 + K 02 V 4 + K 03 V 2 + K 04 ,

where the K are configuration constants, our purpose is to gain the relationship between A and airspeed V, so it is unnecessary to know the expressions of K.

From Eq. (11)-(12), since the A depend on airspeed as well as on configuration, we also get that the resulting pitch motion will also vary with air speed. However, at any speed, Eq. (11) represents a linear differential equation whose solution is of the form:

(13)
α = j = 1 6 α 0 j e s j t ,

where sj are the six roots of the corresponding characteristic equation:

(14)
s 6 + A 5 s 5 + A 4 s 4 + A 3 s 3 + A 2 s 2 + A 1 s + A 0 = 0 .

The roots may be expressed in complex form as follows:

(15)
s 1,2 = g 1 ± i ω 1 , s 3,4 = g 2 ± i ω 2 , s 5,6 = g 3 ± i ω 3 ,

where ωi is modal frequency and gi is modal damping. Substituting Eq. (15) into Eq. (14) and we have:

(16)
A 5 = - 2 i = 1 3 g i , A 4 = i = 1 3 g i 2 + ω i 2 + 4 j = 1 : 3 j i g j , A 3 = - 2 i = 1 , 3 g i j = 1 : 3 j i g j 2 + ω j 2 + 8 k = 1 3 g k , A 2 = i = 1 3 4 g i 2 + ω i 2 j = 1 : 3 j i g j + j = 1 : 3 j i g j 2 + ω j 2 , A 1 = - 2 i = 1 3 g i j = 1 : 3 j i g j 2 + ω j 2 , A 0 = i = 1 3 g i 2 + ω i 2 .

Then the Routh’s stability criterion is used to determine the stability boundary. This is done by forming the test determinants P as follows:

(17)
P 11 = A 5 , P 21 = A 4 - A 3 A 5 ,         P 22 = A 2 - A 1 A 5 , P 31 = A 3 - A 5 P 22 P 21 ,         P 32 = A 1 - A 0 A 5 P 21 , P 41 = P 22 - P 21 P 32 P 31 ,         P 42 = A 0 , P 5 = P 32 - A 0 P 31 P 41 .

The dynamic stability boundary is given by P5= 0. And if P5 is positive, we think the system is stable. At any speed, once the modal parameters ωi and gi are measured, we can calculate the P5 by Eqs. (16), (17). In order to be convenient to polynomial fit, we construct parameters as follow:

(18)
J 1 = A 5 = C 1 V , J 2 = A 4 A 5 - A 3 = C 23 V 3 + C 21 V , J 3 = P 31 J 2 = C 36 V 6 + C 34 V 4 + C 32 V 2 , J 4 = P 41 J 3 = C 410 V 10 + C 48 V 8 + C 46 V 6 + C 44 V 4 + C 42 V 2 , J 5 = P 5 J 2 J 4 = C 518 V 18 + C 516 V 16 + C 514 V 14 + C 512 V 12 + C 510 V 10 + C 58 V 8 + C 56 V 6 + C 54 V 4 ,

where C are the function of K, so they are configuration constants as well. J1-J4 have the same function with T1-T4 that proposed by Price [10], but there is difference between J5 and T5. The T5 is:

(19)
T 5 = T 515 V 15 + T 513 V 13 + T 511 V 11 + T 59 V 9 + T 57 V 7 + T 55 V 5 + T 53 V 3 .

The reason caused the difference is that Price lost a “T51/T2” in his deducing process. We perfect the T5, and plus it with T2, i.e. J2 in Eq. (18), then we obtain J5.

If J5 is positive, we think the system is stable. And when J5 reduced to zero, the flutter is happen. Unfortunately, we can see from Fig. 2 that the variation of test function, i.e. P5, J5 and gα (all normalized with respect to their maximum value), with V2 is not particularly convenient for extrapolating from a subcritical airspeed to the flutter boundary. Therefore, we made an attempt to obtain a new test function which varies in a more convenient manner, and tried several different possibilities. The F3 shown in Fig. 2 is a “better” choice. It should be stressed that all the test functions predicted exactly the same flutter boundary.

Fig. 2. Variation of the normalized test function

Variation of the normalized test function

The test function F3 for 3-dof wing is of the form:

(20)
F 3 = J 5 J 3 J 4 .

Though F3 is not a pure linear function of V2, as shown in Fig. 2, F3 presents to be a convex function whose reduction rate decreases gradually when the airspeed approaching the flutter speed. We can use a linear fitting at low velocity to obtain a conservative prediction that is approximate to flutter. What’s more, because F3 behaves a steady decrease to zero as instability’s approach, it will be helpful to overcome the defect of the velocity-damping method.

Fig. 3. Modal parameters of the three motions at every velocity

Modal parameters of the three motions at every velocity

3. Flutter boundary prediction

Construct a 3-dof wing with a control surface which have the follow structure parameters: ρa= 1.225 kg/m3, b= 0.3 m, ωh= 50.0 rad/s, ωα= 100.0 rad/s, ωβ= 300.0 rad/s, x-α=  0.2, x-β= 0.0125, r-α2= 0.25, r-β2= 0.09, μ= 40, α-= –0.4, c-= 0.6. The flutter was calculated through the p-k method, the flutter boundary is VF= 172.2 m/s, and the graph of V-f and V-g is shown in Fig. 3.

Modal damping calculated by the p-k method can be approximated as the real decay rate of the vibration. So we can calculate F3 at every velocity by the parameters in Fig. 3 through Eq. (16)-(18) and Eq. (20).

Depend on the model we also calculate the flutter margin in the Z-W method [5], which is of the form:

(21)
F = ω 2 2 - ω 1 2 2 + g 2 2 - g 1 2 2 2 + 4 g 1 g 2 ω 2 2 + ω 1 2 2 + 2 g 2 2 + g 1 2 2 2
        - g 2 - g 1 g 2 + g 1 ω 2 2 - ω 1 2 2 + 2 g 2 2 + g 1 2 2 2 2 .

Then the F3 was compared with the F of Z-W method and the damping gα of the pitch motion, as shown in Fig. 4. When the velocity is approach to the flutter boundary, both the F3 and the F reduced gradually. However, while F3 become to be zero at the flutter point, the F is still positive. So the F is not suit for the trinary flutter situation. In addition, although the damping gα decreased to zero as well when flutter is coming, the time its downtrend happened is so closed to the flutter speed that it is not useful to predict flutter boundary earlier.

Fig. 4. Variation of the criterions with velocity squared

 Variation of the criterions with velocity squared

From Fig. 5(a) we can see that F3 isn’t a pure linearity curve. It is nearly a linear function in the front part and present to be a convex function near the flutter boundary, as shown in Fig. 5(b). In the front part, the curve can be extrapolated to be a conservative prediction through linear fitting, and we call it the “first prediction”. F3 will not decrease to zero at this point. After this velocity the curve belongs to the latter part that suit for linear fitting step by step with two or three points as the velocity increase.

Depending on linear fit, the predicted flutter boundary at every velocity is summarized in Table 1. At the front velocities, the predictions approach to 156 m/s and we define this speed to be the “first prediction”. Results of Table 1 show that: when airspeed is below than the “first prediction”, the prediction is toward to the “first prediction”; when airspeed is larger than the “first prediction”, the prediction is reaching true flutter boundary gradually. Though the “first prediction” is not equal to the true flutter boundary, in fact a little less than flutter, it is helpful for test personnel to get an approximate estimation and know how much the airspeed can be increased at a very low velocity.

Fig. 5. Variation of F3 with velocity squared

Variation of F3 with velocity squared

Table 1. Predictions of the flutter speed at all the previous data points (true flutter is 172.2 m/s)

Before first prediction (m/s)
After first prediction (m/s)
Current
Prediction
Current
Prediction
Current
Prediction
Current
Prediction
34.0
181.7
102.0
165.3
150.0
155.0
167.0
170.1
51.0
177.8
119.0
160.7
155.0
160.6
170.0
171.5
68.0
173.8
136.0
156.5
160.0
164.6
171.0
172.0
85.0
167.7
153.0
156.2
164.0
167.8
172.0
172.2

4. Damping and frequency errors for 3-dof wing situation

We have established that good predictions of the flutter speed can be obtained when accurate subcritical parameters are employed. However, because of the influence by noise, atmospheric turbulence, parameter identify method and other interferences, there are recognition errors in parameter identification progress inevitably. In order to decrease the recognition errors, Katz [8] used several sensors to measure modal parameters, and Copper [12] tried to reduce the interference of false modes. Koenig [13] gave exhaustive descriptions of the different sources of uncertainties encountered in flight vibration testing, and he noted that, even after a careful analysis of all available data, final scatter values of up to 2 % for frequencies and 35 % in damping were obtained on the good estimates. Ten years later [14], not much improvement was reported. For instance, a comparison of 16 different analysis methods is presented showing a scatter of approximately 3 % in frequency and 30 % in damping. For this reason, in order to see what is the effect on the flutter prediction of errors in the damping and frequency measurements, a maximum allowable error of 2 % and 5 % in frequency and 20 %, 30 %, 40 %, 50 % in damping were added in the previous case, the variation in F3 with V2 is shown in Fig. 6 and Fig. 7.

Though the variation in F3 presented to have a scatter when there are large errors as shown in Fig. 6(d) and Fig. 7(d), it is convenient to use linear extrapolation to estimate the flutter boundary. On the other hand, comparing Fig. 6(b) with Fig. 6(a), error in damping increased from 20 % to 30 % when error in frequency is 2 %, it is apparent that the error in damping has produced very little deterioration in the variation in F3. Likewise, comparing Fig. 7(a) with Fig. 6(a), error in frequency changed from 2 % to 5 % as error in damping is 20 % all long, it shows that the effect of the error is to produce a fairly large scatter in F3. Because the recognition accuracy of modal frequency is much higher than that of damping, thus this law of errors effect to F3 will help to predict flutter boundary accurately.

Fig. 6. The trend graph of the criterion with 2 % frequency error and the following damping errors

The trend graph of the criterion with 2 % frequency error and the following damping errors

a) 20 %

The trend graph of the criterion with 2 % frequency error and the following damping errors

b) 30 %

The trend graph of the criterion with 2 % frequency error and the following damping errors

c) 40 %

The trend graph of the criterion with 2 % frequency error and the following damping errors

d) 50 %

Fig. 7. The trend graph of the criterion with 5 % frequency error and the following damping errors

The trend graph of the criterion with 5 % frequency error and the following damping errors

a) 20 %

The trend graph of the criterion with 5 % frequency error and the following damping errors

b) 30 %

The trend graph of the criterion with 5 % frequency error and the following damping errors

c) 40 %

The trend graph of the criterion with 5 % frequency error and the following damping errors

d) 50 %

A summary of the effect of the errors is presented in Table 2. Generally, as errors of modal parameters increases, especially error in frequency, the accuracy of the estimated flutter speed decreases. When actual speed is low, like 85 m/s (50 % of the flutter speed) in the case A, we can obtain a flutter prediction. Though this prediction is not equal to the true flutter, it provides us the approximate range of the flutter. As the airspeed increased, like 100 m/s (60 % of the flutter speed) in the case B and 120 m/s (70 % of the flutter speed) in the case C, the predictions is closed to the “first prediction” (156 m/s, shown in Fig. 5) gradually. Meanwhile, when airspeed is larger than the “first prediction”, F3 will decrease to zero as the speed approach to the true flutter. Unfortunately, no matter whether the airspeed is 85 m/s, 100 m/s or 120 m/s, the damping of the three modes in this example have no symptom to reduce, thus will have no help to predict flutter boundary in the velocity-damping method.

Table 2. Effect of errors in frequency and damping measurements on the flutter predictions at the previous speed: A) 85 m/s, B) 100 m/s, C) 120 m/s – the true flutter is 172.2 m/s

Error
Prediction (m/s)
Frequency, %
Damping, %
A
B
C
2
20
154.9
167.5
162.3
30
151.8
159.9
156.5
40
174.1
165.9
162.0
50
210.4
163.9
169.5
5
20
140.9
147.9
162.5
30
171.6
167.1
152.7
40
215.9
166.9
156.3
50
193.6
169.9
163.3

5. Conclusions

A new stability criterion F3 was proposed based on the dynamic equation of the three-dof wing and the Routh’s stability theory in this paper, and its variation with airspeed was described upon simulated experimental data. Further more, the influence of modal parameter errors was discussed through different combinations of maximum allowable errors. The results show that: 1) The criterion F3 can help to predict an approximate flutter boundary with linear fit at low airspeed that will improve the security of flutter test. 2) There is a “first prediction” in the extrapolation process, and the significance of the “first prediction” is that, if airspeed is below than it, the prediction is toward to the “first prediction”, and when airspeed is larger than it, the prediction is reaching true flutter boundary gradually. 3) F3 is shown to be relatively insensitive to errors in modal damping estimations, but is very sensitive to errors in modal frequency measurements. 4) Even in the worst situation, 5 % error in frequency and 50 % error in damping, good estimates of the flutter speed can be obtained from data at speeds as low as 60 % of the flutter speed.

Acknowledgements

This research is partially supported by the National Natural Science Foundation of China (Grant No. 11172128 and No. 51475228), the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20123218110001), the Research Fund of State Key Laboratory of Mechanics and Control of Mechanical Structures (Nanjing University of Aeronautics and Astronautics) (Grant No. 0515G01), the Jiangsu Graduate Training Innovation Project (CXZZ13_0146) and the Priority Academic Program Development of Jiangsu Higher Education Institutions.

References

  1. Dowell E. H., Clark R., Cox D., et al. A Modern Course in Aeroelasticity. Chapter 5. Kluwer Academic, , , 2005. [Search CrossRef]
  2. Kordes E. E. Flutter Testing Techniques. NASA SP-415, 1975. [Search CrossRef]
  3. Walker R., Gupta N. Real-Time Flutter Analysis. NASA CR-170412, 1984. [Search CrossRef]
  4. Park C. G., Friedmann P. P. New time-domain technique for flutter boundary identification. AIAA, Dynamics Specialists Conference, 1992. [Search CrossRef]
  5. Zimmerman N. H., Weissenburger J. T. Prediction of flutter onset speed based on flight testing at subcritical speeds. Journal of Aircraft, Vol. 1, Issue 4, 1964, p. 190-202. [Search CrossRef]
  6. Bennett R. M. Application of Zimmerman Flutter-Margin Criterion to a Wind-Tunnel Model. NASA TM 84545, 1982. [Search CrossRef]
  7. Jennifer H. Stochastic Characterization of Flutter Using Historical Wind Tunnel Data. 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2007-1769, 2007. [Search CrossRef]
  8. Katz H., Foppe F. G., Grossman D. T. F-15 Flight Flutter Test Program. Flutter Testing Techniques, NASA SP-415, 1975, p. 413-431. [Search CrossRef]
  9. Lee B. H. K., Ben-Neticha Z. Analysis of flight flutter test data. Canadian Aeronautics and Space Journal, Vol. 38, Issue 4, 2002, p. 156-163. [Search CrossRef]
  10. Price S. J., Lee B. H. K. Evaluation and extension of the flutter-margin method for flight flutter prediction. Journal of Aircraft, Vol. 30, Issue 3, 1993, p. 395-401. [Search CrossRef]
  11. Theodorsen T. General Theory of Aerodynamic Instability and the Mechanism of Flutter. NACA-496 80N15047, 1979. [Search CrossRef]
  12. Copper J. E. Parameter Estimation Methods for Flight Flutter Testing. Advanced Aeroservoelastic Testing and Data Analysis, AGARD-CP-566, Vol. 10, 1995, p. 1-12. [Search CrossRef]
  13. Koenig K. Problems of System Identification in Flight Vibration Testing. AGARD-R-720, 1984. [Search CrossRef]
  14. Koenig K. Pretension and Reality of Flutter-Relevant Tests. Advanced Aeroservoelastic Testing and Data Analysis, AGARD-CP-566, Vol. 17, 1995, p. 1-12. [Search CrossRef]