Validation of a Robust Flutter Prediction by Optimization
 Author: Chung ChanHoon, SangJoon Shin
 Organization: Chung ChanHoon; SangJoon Shin
 Publish: International Journal Aeronautical and Space Sciences Volume 13, Issue1, p43~57, 30 March 2012

ABSTRACT
In a modern aircraft, there are many variations in its mass, stiffness, and aerodynamic characteristics. Recently, an analytical approach was proposed, and this approach uses the idea of uncertainty to find out the most critical flight flutter boundary due to the variations in such aerodynamic characteristics. An analytical method that has been suggested to predict robust stability is the mu method. We previously analyzed the robust flutter boundary by using the mu method, and in that study, aerodynamic variations in the Mach number, atmospheric density, and flight speed were taken into consideration. The authors’ previous attempt and the results are currently quoted as varying Mach number mu analysis. In the author’s previous method, when the initial flight conditions were located far from the nominal flutter boundary, conservative predictions were obtained. However, relationships among those aerodynamic parameters were not applied. Thus, the varying Mach number mu analysis results required validation. Using an optimization approach, the varying Mach number mu analysis was found out to be capable of capturing a reasonable robust flutter boundary, i.e., with a low percentage difference from boundaries that were obtained by optimization. Regarding the optimization approach, a discrete nominal flutter boundary is to be obtained in advance, and based on that boundary, an interpolated function was established. Thus, the optimization approach required more computational effort for a larger number of uncertainty variables. And, this produced results similar to those from the mu method which had lower computational complexity. Thus, during the estimation of robust aeroelastic stability, the mu method was regarded as more efficient than the optimization method was. The mu method predicts reasonable results when an initial condition is located near the nominal flutter boundary, but it does not consider the relationships that are among the aerodynamic parameters, and its predictions are not very accurate when the initial condition is located far from the nominal flutter boundary. In order to provide predictions that are more accurate, the relationships among the uncertainties should also be included in the mu method.

KEYWORD
Robust aeroelasticity , Worst case flutter boundary , Uncertainty , Unsteady aerodynamics , Optimization

Nomenclature
MS = Generalized mass matrix
KS = Generalized stiffness matrix
Aij = Structural coupling matrix
Bij = Bending coefficient matrix
Tii = Torsional coefficient matrix
ψi = Bending mode shape function
Θi = Torsional mode shape matrix
Ξwi = Generalized threedimensional lift
Ξθi = Generalized threedimensional pitching moment
L？ = Twodimensional lift
M？1/4 = Twodimensional pitching moment
AQ0, BQ0, CQ0, DQ0= Aerodynamic state matrix
C(k) = Theodorsen’s lift deficiency function
AM = Generalized aerodynamic mass matrix
AC = Generalized aerodynamic damping matrix
AK = Generalized aerodynamic stiffness matrix
Ax = Generalized aerodynamic lag state matrix
WU, Wρ = Weighting value for the flight speed and atmospheric density
WA, WB, WC, WD = Weighting matrix for the state matrix
P = Robust aeroelastic system matrix
Zw = Uncertainty output
Ww = Uncertainty input
r(t) = Nominal aeroelastic flutter boundary function in terms of the discrete Mach number, altitude, and flight speed
M(t) = Mach number variation function in terms of the discrete parameters
H(t) = Altitude variation function in terms of the discrete parameters
V(t) = Flight speed variation function in terms of the discrete parameters
Aij, upper, Aij, lower, Bij, upper, Bij, lower, Cij, upper, Cij, lower, Dij, upper, Dij, lower = Upper and lower bounds of the state matrix
Arobust, Brobust, Crobust, Drobust = Robust aerodynamic state matrices
unew = State input for an aerodynamic state equation
uwnew = State lift input for an aerodynamic state equation
uθnew = State pitching moment input for aerodynamic state equation
S1, S2 = Distance between each initial condition and the nominal flutter boundary
d1, d2 = Distance in altitude between each initial condition and the nominal flutter boundary
h1, h2 = Distance between each mu prediction and the optimization result regarding the robust flutter boundary
a = Structural unbalance factor
b = Semichord length
1. Introduction
Advanced maneuverability is generally required for a modern combat aircraft to complete mission segments such as severe maneuvers. The aerodynamics of the wing varies under different flight conditions such as with different flight speeds and atmospheric densities. Wing structural characteristics also vary according to their mass and stiffness. Moreover, structural characteristics may reflect on the external store requirements. For example, the mode shape and natural frequencies of a wing can vary with the location of the store, and it may also change with a change in mass such as during the launch of a missile. The shape and location of the store item will also affect the wing aerodynamics. The study of such variations is important as these changes may be unpredictable. Thus, for reliable prediction of an aeroelastic stability boundary such as flutter, the above parameters should be included. However, an accurate prediction of a combat aircraft’s flutter margin is difficult practically and this is due to their high maneuverability and rapidly changing aerodynamic environment. Due to the difficulties mentioned during aircraft design, a flutter margin of 15% is generally adopted as a solution to this difficult problem. Thus, for aeroelastic analysis, a more accurate and reliable method that considers several possible sources of variation is required.
Recently, a new approach to aeroelastic analysis based on robust control theory and the mu method was used to determine the robust aeroelastic stability boundary [1,2]. The mu method is reported to be effective for the prediction of the stability boundary of an uncertain system. In particular, it obtains a more conservative result compared to the nominal flutter boundary. Thus, the use of the mu method has been suggested as a predictor for robust aeroelastic systems [1,2]. Since then, many robust analysis results have been reported. A robust flutter boundary with varying dynamic pressures was investigated by Lind [1], who found out that the results of the robust flutter margin were about 4% more conservative than the nominal flutter boundary. A frequency domain analysis was performed by Borglund [3,4], who presented the mu value and reduced frequency results. Moreover, an uncertainty quantification method has been suggested by Kurdi and Lindsley [5]. Aerodynamic uncertainty was considered because accurate prediction capability was needed under the varying aerodynamics environment. Practically, many factors which cause the variation in aerodynamic forces exist. These factors are considered as uncertainties and they cannot be predicted in a systematic manner. Recently, researchers have suggested a few methods for uncertainty quantification. Yuting showed that aerodynamic uncertainties were predicted from the differences between the doublet lattice method (DLM) and Computational Fluid Dynamics (CFD) [6]. In fact, the boundaries of the uncertainty were defined with the maximum and minimum differences between DLM and CFD.
In order to predict flutter boundary conservatively, aerodynamic uncertainties need be quantified and included in the analysis appropriately. The pilot also has to be aware of the flutter boundary at the current flight state.
Since such uncertainties cannot be predicted in a systematic manner, approximate arbitrary margin will be applied. The mu method has been suggested as one of the methods for the compensation of such limitations. However, the mu method can be applied only for a single variable problem. Therefore, all the preceding analyses were done under the assumption of a constant Mach number. Practically an aircraft experiences complex flight conditions such as varying Mach number, atmospheric density, and flight speed. Thus, the previous analyses under a constant Mach number assumption have a limitation, and a more complete analysis that considers variations in the Mach number. This is required for a more general and accurate prediction.
The mu method, when using a constant Mach number approach, produces conservative results as shown in Fig. 1 and this is adapted from Borglund [4]. However, the constant Mach number analysis is limited as it does not consider the structural and aerodynamic force variations within an aircraft. For example, Mach number and altitude variations affect aerodynamic forces when the aircraft does a maneuver. Thus, both the Mach number and altitude aspects of the flutter boundary problem should be analyzed.
The authors’ previous attempt and results are currently quoted as varying Mach number mu analysis. In such varying Mach number mu analysis, the mu analysis was performed by varying Mach numbers, altitudes, flight speeds, and atmospheric densities, and variations in aerodynamic forces within a certain flight envelope were represented as uncertainties [7]. In that analysis, the relevant variation ranges were estimated and they were prescribed by using the doublet lattice method (DLM). By comparing the resulting aerodynamics with Theodorsen’s function, the validity of a range was proved. Subsequently, a robust aeroelastic governing equation was constructed, and the mu method was applied to find out the robust flutter boundary. In a constant Mach number analysis, the mu method predicted a 4~5% conservative result that was relative to the nominal flutter boundary and only the structural uncertainties were included[2]. In the varying Mach number mu analysis, variations of the flutter Mach number are increased by 0.8？ 2.2% over the nominal value [7]. As mentioned, variations in the Mach number, atmospheric density, and flight speed were taken into consideration in that study. During analysis, two of those three variables were deemed dependent and
one was independent. The relationship between the two dependent variables should have been considered during mu analysis. However, the relationship was not considered in the varying Mach number mu analysis. For this reason, reliable worst case flutter results were obtained only when the initial condition was located near the nominal flutter boundary. However, the method predicted a more conservative result when the initial condition was located far from the nominal flutter boundary. Thus, the result of the mu analysis in the varying Mach number mu analysis work still required validation. The previous studies of other researchers’ have attempted to validate the results of the mu analysis by evaluating the eigenvalues at the obtained flutter boundary [4]. However, the validation was done under constant Mach number flight conditions. Thus, when the Mach number which corresponds to the worst case flutter boundary is varied from a value at a given initial condition, that validation exhibits a limitation. Here, we present a new and more complete validation of the mu analysis result.
In the present paper, we focus on the validation of the positive boundary, i.e., on obtaining a conservative result. The mu method is an optimization method that considers uncertainty. Therefore, our validation of the mu analysis uses an optimization approach. First, the nominal flutter boundary results in the varying Mach number mu analysis and they were validated at specific initial conditions. The nominal flutter boundary obtained previously by the authors [7] was validated with those by Goland [8] and Brown [9]. Based on the previous result for the nominal flutter boundary, an interpolation function to represent the nominal flutter boundary was then obtained. Subsequently, an optimization approach was used to find out the nearest point to the interpolated flutter boundary which is socalled as the most critical flutter boundary. Results obtained by using this optimization approach were then compared with those obtained by the varying Mach number mu analysis.
As a result, the comparison result and differences between the nominal flutter, mu analysis, and validation result are summarized at numerical results.
The mu method results are presented to be 3.2 ~ 4.3 % more conservative than the nominal flutter boundary. The Mach number is not changed in nominal flutter result, but the Mach number is changed in varying Mach number mu analysis result. As shown in the present numerical results, the varying Mach number mu analysis result shows good agreement with the Optimization method. Despite, a new method for flutter prediction is required at varying Mach number condition, only few have attempted to apply the mu method. Moreover, there are few papers that consider the validation of the mu result. In the present paper, the variation results of varying Mach mu method result are shown. Based on the optimization method, validation is performed.
2. Formulation
2.1 Structural Modeling
In the present aeroelastic analysis, Goland’s wing is used (Table 1). The present Goland’s wing is illustrated in Figure 2.
As shown in Figure 2, the present Goland’s wing is considered as a cantilevered beam and it includes static unbalance. Even thoughnt wing has NACA 0009 airfoil, its thickness effect was ignored for simplicity.
Eqs. (1) and (2) express the mass and stiffness matrices respectively and they were used in the present analysis,[7,10].
In the above equations, [
A ] is the modal coupling matrix due to static unbalance in the wing cross section, while the matrices [B ] and [T ] contain the bending and torsion stiffness coefficients, respectively. In the present paper, the first bending and the first torsion mode were only considered. Thus, the present structural model utilized two degrees of freedom. Detailed formulation of this structural analysis model is described by Chung [7].2.2 Aerodynamic Modeling
The threedimensional aerodynamic force is derived from the twodimensional lift and pitching moment [10]. The threedimensional generalized aerodynamic force is
expressed as follows,
In the above equation,
L´ and [M´ _{1/4} + (1/2 +a )bL´ ] is the twodimensional lift and pitching moment, respectively. Theodorsen’s functionC (k ) is included in both of them. The threedimensional lift and aerodynamic pitching moment can be obtained by integration.In the varying Mach number mu analysis, a statespace model was suggested by using a rational function approximation (RFA) [5] and in that paper, the DLM and RFA results were compared with respect to the Mach number. However, this approach showed a slight difference between the two results for large values of reduced frequency. Due to this limitation, the following frequency domain approach was used to obtain predictions that are even more accurate. In the frequency domain analysis, aerodynamic variation will be formulated in terms of reduced frequency. The detailed formulation of the aerodynamic analysis and the relevant RFA can be found in Chung [7].
2.3 Aeroelastic Equation
The aeroelastic formulation can be obtained by an equation that was used to analyze the general mechanical vibration [7]. The governing equation can be written as follows.
3. Robust Aeroelastic Analysis Considering Aerodynamic Variation
The mu method which is adopted in the present paper has limitations upon an object with large number of D.O.F. Due to this, the mu method is based on LMI theory. The present structural matrices are still expressed with two degrees of freedom. The aeroelastic system matrix is a combination of the structure and aerodynamic models and it features seven degrees of freedom which includes three aerodynamic lag state variables. Thus, it will be possible to apply the mu method for the present problem.
Our varying Mach number mu analysis took aerodynamic variations into account [7]. Hence, the mu method was adopted, and the aerodynamic force was assumed to be uncertain. Modeling of the aerodynamic variations and the relevant formulation were done as follows. First, the variations that are related to the flight speed and atmospheric density were formulated, respectively, as
Where,
W_{U} andW_{ρ} are the weighted scalar values that prescribe the ranges of the variations. Due to the inclusion of variations in the flight speed and atmospheric density, the aerodynamic state matrix is time varying with a weighted factor. In our varying Mach number mu analysis, the aerodynamic state matrix was obtained by using RFA, and its variations were defined by using a weighting matrix [7]. Here, the statespace model of aerodynamic forces is given as,Where, matrices
A_{Q0} andB_{Q0} are the aerodynamic system and input matrices, respectively. The output of the aerodynamic vector,y is represented in terms ofC_{Q0} andD_{Q0} matrices. The vector,x denotes an aerodynamic lag state. In Eq. (9), the output vector is further represented by a form that is similar to the Theodorsen’s lift deficiency function as,In Eqs. (8) ？ (10),
u_{new} is represented as follows,The components
u_{wnew} andu_{θnew} were previously defined by Kurdi (2007). Its details are given in Appendix.The weighting matrices are expressed as,
A final statespace form of the aeroelastic equation for the mu analysis was obtained as follows,
In Eq. (16), the nominal aeroelastic equation and the variation denoted by
W_{w} are coupled. In a relevant analysis, the altitude and the Mach number were prescribed by Chung [7].4. Optimization
4.1 Physical Meaning of the Most Critical Flutter Boundary
In the varying Mach mu analysis, we attempted to analyze the most critical flutter boundary from a robust aeroelastic system by using the mu method [7]. In the previous study by the authors, the Mach number, atmospheric density, and flight speed were considered as uncertainties. Moreover, there were two dependent variables and one independent variable but the relationship between the dependent variables was not taken into consideration. Thus, the varying Mach number mu analysis flutter prediction became conservative. When the initial condition was located near or at the nominal flutter boundary, the varying Mach number mu analysis result was deemed to be reasonable. However when the initial condition was located at a greater distance from the nominal flutter boundary, the varying Mach number mu results were conservative. Thus, the varying Mach number mu results required validation through an alternate analysis. The optimization analysis used here is based on the nominal flutter boundary and it was predicted with respect to each discrete aerodynamic condition. From each nominal flutter boundary, an interpolation function will be described. In the present validation, a socalled most critical flutter boundary term is used and this term indicates a flutter boundary that is nearest to a given initial condition. The objective of this optimization was to find out the most critical flutter boundary and compare it with the robust flutter prediction result that was obtained by using the mu method. This comparison is possible because the present optimization was capable of predicting a robust flutter boundary and this is regardless of the location of the initial condition. Physically, this implies that the most critical boundary is at a location where the flutter occurs with the least amount of kinetic energy in the fluid particles which flow nearby. This most critical boundary corresponds to a lower flutter flight speed, a lower flutter atmospheric density (i.e., higher flutter altitude), and a lower flutter Mach number. Geometrically, the most critical boundary on the nominal flutter boundary has the shortest distance from a specific initial condition. Repetition of the process under different initial conditions produced a line that was drawn normal to the nominal flutter boundary. Such points can be determined by using a minimization process. In the varying Mach number investigation based on the mu method, the physical meaning of the robust flutter boundary was set as the most critical flutter boundary after the consideration of all the possible variations. Thus, both methods should yield the same result regarding the worst case flutter boundary.
The concept of the present optimization is illustrated as shown in Fig. 3. The solid line denotes the nominal flutter boundary with respect to both the dynamic pressure and altitude. When the nominal flutter boundary is expressed in terms of a specific Mach number, an interpolation function of the nominal flutter boundary can be obtained. When an initial condition is given, the most critical flutter boundary can be obtained by minimizing the distance between the initial condition at a given Mach and the interpolated nominal flutter boundary.
In Fig. 3, the black square denotes a prescribed initial condition at Mach 1. The optimization method finds the nearest point and it is denoted as a whitefilled square, along the nominal flutter boundary. The normal direction from an initial condition to the nominal flutter boundary is to be considered. As shown in the figure, a change in the Mach number occurs while finding out the nearest point for a specific initial condition through optimization.
4.2 Interpolation of the Nominal Flutter Boundary
In the present validation, all the relevant parameters are to be represented in the nominal flutter boundary. This requires additional computations which may not be practical or effective. Thus, interpolation was adopted in the current approach. Practically, the discrete nominal flutter boundaries will be expressed in terms of the Mach number, altitude, and flight speed. The nominal flutter boundary is represented in the threedimensional space as,
Where,
The distance
r (t ) is then constructed in terms of the Mach numberM (t ), altitudeH (t ), and flight speedS (t ). In the above equations, the parameter t, which is the number of discrete parameters in the flutter boundary is introduced. Comparisons between the previously obtained nominal flutter boundary and the interpolated function obtained values that are shown in Fig. 4(a) ？ (c). The interpolation functions are illustrated in terms of the parameter, t to describe Mach number, altitude, and flight speed. The interpolation functions for the Mach number vs. t, altitude vs. t, and flight speed vs. t were obtained based on the functions that are shown in Fig. 4. The coefficients shown in Eqs. (18) ？ (20) were obtained from MATLAB. By using these coefficients, a continuous nominal flutter boundary was obtained. Based on Eqs. (18)  (20), the nominal flutter boundary curve expressed in Eq. (17) was obtained, and it is shown in Fig. 5 (a). In Fig. 5 (b) ？ (d), the functions that represent Mach vs. altitude, Mach vs. speed, and altitude vs. speed, are shown. And these figures are obtained by projecting Fig. 5 (a) into each twodimensional plane.4.3 Optimization
The most critical flutter boundary can be obtained by the optimization with a formulation and it is expressed as follows,
The performance index function indicates a distance from a certain initial condition and the variables with the subscript zero are to be minimized. The minimized value of the function was computed by using the “fminunc” command that was provided in the MATLAB Optimization Toolbox [11]. Other constraints were not considered, but an automatic constraint of t < 10 was used due to the discrete point number.
5. Numerical Result
5.1 Characteristics of the ThreeDimensional Wing
Goland’s wing was used for the present validation. Its characteristic values are summarized in Table 1.
5.2 Robust Aeroelastic Analysis Validation Result
Results of the robust aeroelastic stability analysis by using
the varying Mach number mu method and those of the present optimization method are summarized in Table 2. Three different initial conditions were used. In the first case, the initial Mach number was 0.5, and the initial altitudes were 20,000, 25,000, 30,000, and 40,000 ft. In the second case, the altitudes were 45,000, 50,000, and 60,000 ft with an initial Mach number of 0.6. In the third case, the altitudes were 55,000 ft, 60,000 ft, and 65,000 ft with an initial Mach number of 0.7.
In the present robust analysis results, the values of the weighting matrix were selected and they were according to the aerodynamic regime. In the present analysis, the values for the weighting matrices were chosen from the difference between the upper and lower bounds in the aerodynamic state matrices. The upper and lower bounds were selected by an approximation on the Theodorsen’s function, by applying RFA when the upper and lower Mach number was prescribed.
In the present paper, the results of the authors’ varying Mach number mu results were revised with an updated aerodynamic regime, and the updated results were validated by using the optimization method. The presently revised and validated results are summarized in Table 2.
In Table 2 the updated aerodynamic regime was defined with the following assumptions: (1) the initial Mach number must be located at the middle point of the Mach number region. (2) The initial altitude was given at the middle point of the defined altitude regime. The present revised mu results are shown in Table 2. These upgrade results are obtained by using the new aerodynamic regime. The new aerodynamic regimes that were used in the present research were defined with new Mach number and new altitude regimes. In the authors’ varying Mach number mu analysis, Mach number regime was assumed to be 0.5 ~ 0.7 and this was regardless of the initial Mach number. On the other hand, the initial Mach numbers for the present revised results were located at the middle point in the new Mach number regime. In the above conditions, the regime of altitude was allowed between 10% ~ +10% for each initial altitude.
In Table 2, the robust flutter boundary was predicted to be more conservative in each given initial condition. In the first case, the initial Mach number was 0.5, and the initial altitudes were 20,000, 25,000, 30,000, and 40,000 ft. In the second case, the altitudes were 45,000, 50,000, and 60,000 ft with an initial Mach number of 0.6. In the third case, the altitudes were 55,000 ft, 60,000 ft, and 65,000 ft with an initial Mach number of 0.7. In the first case, the new Mach number regime was defined to be 0.5 ~ 0.76. Moreover, the weighting value of the Mach number considered to be twice of differences between Mach number 0.5 ~ 0.76. The flutter speeds at the altitude are obtained by using mu analysis and they ranged from 518.5 to 529.14 ft/sec (Table 2a). Compared to the results that are obtained by using the nominal flutter prediction, the mu flutter speeds results have differences of about 1.3 to 2.8%. By using the mu method, the flutter altitude ranged from 19,662 to 38,467 ft, and the corresponding dynamic pressure range was from 179.43 to 86.51 slugs/ft3 (Table 2a). The mu results were validated with the optimization method. By using the optimization, the flutter altitudes were 2.2 to 4.84% lower compared to those from the mu method. Furthermore, the Mach numbers from the mu method ranged from 0.503 to 0.506 while the range that uses the optimization was from 0.50 to 0.570. Therefore, compared to the results that were obtained by optimization, the mu results were in the range of 12.8 to +1.0% .
In the second case, the Mach number regime considered 0.5 ~ 0.76, and the weight factor of the Mach number was defined with twice the differences between Mach number 0.6 and 0.76. The flutter speeds at altitude obtained by the mu analysis were from 604 to 610.7 ft/sec. It could be compared to the flutter speeds by optimization, 607 to 703 ft/sec. (Table 2b). The mu results were +0.45 to ？16.3% compared to those obtained by using optimization under the same conditions. By using the mu method, the flutter altitudes ranged from 44,946 to 58,123 ft, and the corresponding dynamic pressures ranged from 91.78 to 49.49 slugs/ft3. By using optimization, the flutter altitudes were from 0.48 to 61.33% compared to those [obtained?] from the mu method. Moreover, Mach numbers from the mu method ranged from 0.601？0.603, while those from the optimization ranged from 0.593 to 0.690. Relative to the results obtained by optimization, the mu results were 14.8 to +1.67%. In the third case, the Mach number range was the same, but the weighting factor was considered as twice the difference between the Mach number 0.5 ~ 0.7. The flutter speeds obtained by the mu analysis ranged from 700.63 to 703.36ft/sec (Table 2c). The mu methodderived flutter speeds were from 0.43 to +1% compared to those which were obtained by the optimization method under the same conditions. By using the mu method, the flutter altitudes ranged from 54,940 to 63,030 ft, and the corresponding dynamic pressure range was from 78.11 to 52.03 slugs/ft3. The flutter altitudes from the mu method were from 52.9 to 1.28% compared to those that were derived by optimization. Moreover, by using the mu method, the Mach number ranges from 0.7 to 0.704 while the Mach number calculated by using optimization was consistently 0.69. The mu Mach numbers were from 0 to 2% compared to those which were obtained by optimization.
Overall, the mu and the optimization results in Table 2 show good agreement over most of the flight conditions that
were analyzed. When the initial Mach number was relatively low, the agreement between the two predictions was close but even when the initial Mach number was large, the results from both the methods were still in good agreement. It showed a better agreement when the initial condition was located near the nominal flutter boundary. However, the most critical flutter boundary predicted by the mu method becomes conservative when the initial condition is located far from the nominal flutter boundary. Under these conditions, the two results show a significant discrepancy.
The present paper focuses on the varying Mach number analysis and the corresponding result of the constant Mach number analysis that were obtained by Lind (1999) are presented in Table 2. The robust aeroelastic formulation used by Lind (1999) was constructed by taking into consideration the structural uncertainties, and the worst case flutter
boundary that was predicted by using the mu method. When the initial Mach numbers are given as 0.5, 0.6, 0.7 and the differences in the flutter speed prediction between the constant and varying Mach number analysis were from 0.9 ~ 0.97%, 0.24 ~ 0.2 %, and 0.43 ~ 0.0 %respectively. Lind considered only the structural uncertainties in his constant Mach number analysis while the present analysis considered aerodynamic uncertainties under varying Mach numbers. Therefore, even though it may be meaningful the present comparison between the two results has limitations as both the analyses were based on the same mu method.
Based on the above comparison, the previously reported mu analysis method was validated by using this optimization approach. The comparison suggests that the mu method is capable of predicting a reasonable robust flutter boundary with less than a few percent differences from that which was derived by the optimization method. The previously described mu method does not require the additional computations that are needed in the optimization approach, such as Eqs. (17) ？ (21) for the creation of the interpolation functions. These computations are not required as the parameters in those equations were already considered as uncertainties and were used during the mu analysis. Thus, we conclude that the mu method is a more efficient analytical method for the determination of the robust aeroelastic stability of an aircraft.
The optimization method that is described here attempts to find out the shortest distance between an initial condition and a continuous nominal flutter boundary by minimizing the specific performance index function. As a result in most of the cases, the prediction that is obtained by such an optimization is located closer to the nominal flutter boundary than the result of the mu analysis. Thus, we deduce that our optimization approach predicts a more reasonable result compared to that from the varying Mach number mu method. However, the present optimization requires more computational efforts than the varying Mach number mu method. For example, optimization requires preliminary information regarding the nominal flutter boundary which
is as predicted by a separate reliable analysis, and an interpolation function for the nominal flutter boundary for each varying parameter that has to be obtained in advance. The performance index function is defined as the distance from an initial condition to the nominal flutter boundary from the interpolated function. During the optimization analysis, this performance index has to be evaluated and minimized. Due to this additional computational procedure, the optimization used here will consume more computation time compared to that of the mu method and it may even reach a numerical limitation when the number of variation parameters is increased.
Our varying Mach number mu analysis indicated that the mu method gave a conservative result when the initial conditions were located far from the nominal flutter boundary. However, the present optimization provides an accurate prediction of the most critical flutter boundary and this is regardless of the magnitude of the variations. In our above numerical treatments, when the initial condition is located far from the nominal flutter boundary, the difference between the values that were derived by the two methods becomes larger. In the present comparison, the flutter speed predictions showed up to 12 ~ 16% differences at the highest altitudes for the three Mach numbers tested while the dynamic pressure calculations were up to 13 ~ 67% different at the highest altitudes for the three Mach numbers tested. However, when the initial conditions are located near the nominal flutter boundary, the differences in the flutter speed prediction decreased to 0.8 ~ 3%, and it was also observed that the differences in the dynamic pressures was 0 ~ 1.18 %. In Fig. 6, Mach numbers 1 and 2 correspond to theers 0.5 and 0.6, respectively and they are presented in Table 2.
In the figure 6, initial condition 1 and Mach number 1 illustrate a case when theondition is located near the nominal flutter boundary at an initial Mach number of 0.5. On the other hand, theondition 2 and Mach number 2 denote a case when theonditions are located far from the nominal flutter boundary at an initial Mach number of 0.6. Distances d1 and d2 are the altitude distances between the initial altitude and the nominal flutter boundary while p1 and p2 are the dynamic pressure differences between eachondition and the nominal flutter dynamic pressure. The distances h1 and h2 are estimated to be approximately 1,828 and 120.0 ft, respectively and the corresponding dynamic pressure distances p1 and p2 are approximately 18.28 and 0.44 slug/ft2？sec, respectively. With regard to the flutter dynamic pressure in the former case, both the mu method and the optimization method yielded similar boundaries. However, the latter case in which there was a larger difference between theondition and the nominal flutter boundary, showed a marked difference. Fig. 7 provides a more detailed view of the parameters that are presented in Fig. 6 and it is based on the results which are presented in Table 2(b).
Fig. 7 shows the comparison of the two results (mu and
optimizationbased) for an initial Mach number of 0.6. When the altitude distance between the initial condition and the nominal flutter boundary was 15,000 ft, the dynamic pressure difference was 49.49 slug/ft2？sec and the distance S1 was 15,000 ft. In this example, the distance h1 between the flutter boundary predicted by the optimization method and that predicted by the mu method was approximately 3,634ft. However, when the altitude distance between the initial condition and the nominal flutter boundary was 1,984 ft, the dynamic pressure difference was only, 0.44 slug/ft2？sec and this resulted in a geometrically different distance S2, of approximately 1,984 ft. In that case, the distance h2 between the mu method and optimization method predictions decreased by about 45%.
In the present paper, the worst case flutter boundary was obtained by varying the altitude during the mu analysis. When the worst case flutter altitude was obtained, the corresponding worst case flutter speed and the atmospheric density can also be predicted. The worst case flutter Mach number can be obtained by applying the relationship between the standard atmospheric density and Mach number. However in this paper, such a relationship has not yet been considered. If the relationship among the aerodynamic factors such as Mach number, flight speed, and atmospheric density were considered in the present mu analysis then, the uncertainty effect will be further reduced. Moreover, the finally updated results will be closer to the nominal flutter boundary and the present optimization result. Thus, an improved procedure that updates the relevant functional relationships in the mu analysis has to be developed and applied. Moreover, as mentioned in the previous section, a frequencydomain analysis approach can be adopted to represent the relevant aerodynamics and obtain results that are far more accurate.
6. Conclusion
In the present paper, prediction results for the robust flutter boundary were introduced in a varying Mach number mu study by the authors is revised with new aerodynamic regime and is validated by using the optimization method. In the varying Mach number mu study, robust flutter boundary results were obtained by using the mu method. In that method, aerodynamic variation was prescribed by using a concept of uncertainty. However, the relationships among the variables were not considered during the mu analysis. In the mu analysis, when the initial condition was located far from the nominal flutter boundary, a conservative result was predicted. Thus, a thorough validation based on a separate analysis was deemed necessary by the usage of the optimization method. In order to apply the idea of optimization, the nominal flutter boundary needs to be analyzed, and an interpolation is required. By using optimization, the most critical flutter boundary which is the boundary located at the shortest distance from a given initial condition can be obtained. The robust flutter result previously obtained by using the mu method was then compared with that obtained by the optimization method. According to the present numerical results, the varying Mach number mu analysis has been validated by using the optimization analysis. First, the varying Mach number mu method is deemed capable of predicting a robust flutter boundary with low percentage differences from that obtained from the present optimization. Moreover, the varying Mach number mu method requires less computational time compared to that of the optimization method. This is due to a preliminary nominal flutter boundary prediction and its subsequent interpolation is not required during the mu analysis. Thus, we conclude that compared to the optimization method, the mu method is a more efficient analysis to find out the robust aeroelastic stability. When the number of the uncertainty parameters is increased, the present optimization will require more computational effort and time. However, even with an increase in the number of parameters, the mu method will involve less computational effort. Second, in the present paper, we have also observed that the mu method is capable of correctly predicting the nearest point only when the initial flight condition is located near the nominal flutter boundary. This is due to the interparameter variations that were not updated during the varying Mach number mu analysis. Updation of these variations should also be included in the mu analysis to compensate for this limitation. We suggest that if the relationships among the various aerodynamic parameters are taken into consideration, the mu method will be able to more reliably predict reasonable results. Thus, a more refined procedure that considers the relationships among the various flight parameters has to be investigated. It is expected that such studies will develop new suggestions for the future improvements of the mu method.
> Appendix
A. Structural Matrices
A.1 Mass Matrix
The mass matrix for the present wing is represented as,
Where, the matrix [
A ] denotes the structural mode coupling and its elements are expressed as follows.In Eq. (A.2), Θ_{i}, ψ_{i} denote the torsion and bending mode shape functions that are obtained as
A.2 Stiffness Matrix
The stiffness matrix is based on the bending and torsion stiffness.
Where, [
B ], [T ] denotes bending and torsion factor matrices and its element is obtained as follows.B. Aerodynamic Forces
B.1 Twodimensional Aerodynamic Force
The threedimensional aerodynamic force that is used in the present paper is obtained from the twodimensional lift and pitching moment. Twodimensional lift and pitching moment are expressed as follows.
B.2 Threedimensional Aerodynamic Force
The threedimensional unsteady aerodynamic forces are based on the twodimensional aerodynamic forces, and they can be obtained by integration as follows
The final expression for the aerodynamic forces is obtained as follows,
B.3 Rational Function Approximation (RFA)
For the aeroelastic equation in a statespace form, rational function approximation is adopted and its equation is given as follows.
In the present paper, the Theodorsen’s function is approximated by RFA, and it can be represented as the following state equation.
In the above equation,
u is an aerodynamic input and it can be represented as follows.C. Robust Aeroelastic Analysis
The present robust aeroelastic equation is based on the variation in aerodynamic force, and the aerodynamic forces can be expressed as follows,
The approximated Theodorsen’s function can be applied for the relevant robust aeroelastic analysis. The improved Theodorsen’s function and aerodynamic input can be derived as follows,
The corresponding state matrices are to be revised properly for the foregoing robust aeroelastic analysis
Based on the above derivation, the uncertainty matrix is obtained as,
The final form of the robust aeroelastic analysis equation that combines the aerodynamic uncertainty and nominal aeroelastic equation, is derived as,
Where,
In matrix
P, and
are the lift and pitching moment coefficients, respectively and they consider aerodynamic variation.
represents an aerodynamic lagstate variable.

[Fig. 1.] Nominal and Robust Flutter Margin at Constant Mach Number [4]

[Fig. 2.] A Threedimensional Cantilevered Wing [7]

[Fig. 3.] Nominal and Most Critical Flutter Boundaries

[Fig. 4.] Interpolation Functions and Nominal flutter boundary

[Fig. 5.] Interpolated and Nominal Flutter Boundaries in the Three and Twodimensional Space

[Table 1.] Characteristic Values of a ThreeDimensional Wing

[Table 21.] Validation of the Most Critical Flutter Analysis

[Table 22.] (b) Mach = 0.6

[Table 23.] (c) Mach = 0.7

[Fig. 6.] Change in the Flutter Dynamic Pressure with respect to an Initial Altitude

[Fig. 7.] Change in the Flutter Dynamic Pressure and Altitude with respect to an Initial Condition Location