Abstract
The present work describes an alternative methodology for identification of aeroelastic stability in a range of varying parameters. Analysis is performed in time domain based on Lyapunov stability and solved by convex optimization algorithms. The theory is outlined and simulations are carried out on a benchmark system to illustrate the method. The classical methodology with the analysis of the system's eigenvalues is presented for comparing the results and validating the approach. The aeroelastic model is represented in state space format and the unsteady aerodynamic forces are written in time domain using rational function approximation. The problem is formulated as a polytopic differential inclusion system and the conceptual idea can be used in two different applications. In the first application the method verifies the aeroelastic stability in a range of air density (or its equivalent altitude range). In the second one, the stability is verified for a rage of velocities. These analyses are in contrast to the classical discrete analysis performed at fixed air density/velocity values. It is shown that this method is efficient to identify stability regions in the flight envelope and it offers promise for robust flutter identification.
Flutter analysis; LMI; Lyapunov function; Time domain
TECHNICAL PAPERS
IIITechnological Institute of Aeronautics - São José dos Campos - 12228-900 SP - Brazil - goes@ita.br
IVUniversidade Estadual Paulista - Bauru - 17033-360 SP - Brazil paulo.jpg@feb.unesp.br
ABSTRACT
The present work describes an alternative methodology for identification of aeroelastic stability in a range of varying parameters. Analysis is performed in time domain based on Lyapunov stability and solved by convex optimization algorithms. The theory is outlined and simulations are carried out on a benchmark system to illustrate the method. The classical methodology with the analysis of the system's eigenvalues is presented for comparing the results and validating the approach. The aeroelastic model is represented in state space format and the unsteady aerodynamic forces are written in time domain using rational function approximation. The problem is formulated as a polytopic differential inclusion system and the conceptual idea can be used in two different applications. In the first application the method verifies the aeroelastic stability in a range of air density (or its equivalent altitude range). In the second one, the stability is verified for a rage of velocities. These analyses are in contrast to the classical discrete analysis performed at fixed air density/velocity values. It is shown that this method is efficient to identify stability regions in the flight envelope and it offers promise for robust flutter identification.
Keywords: Flutter analysis, LMI, Lyapunov function, Time domain
Introduction
In problems of fluid and flexible structure interactions various physical phenomena can be observed and also, because of its catastrophic nature, flutter is an important topic involving aeroelastic stability. This phenomenon is an interaction between the structural dynamics and the aerodynamics that results in divergent and destructive oscillations of motion (Bisplinghoff et al., 1996).
Different approaches have been proposed to identify the flutter boundaries. In general, the methods are formulated in frequency domain as eigenvalue problems, for which the aeroelastic model is defined for each point in the flight envelope (a pair of altitude and velocity). Details can be found in Theodorsen and Garrick, 1940; Hassig, 1971; and Chen, 2000. This can be costly in terms of computational time and engineering analysis of data if a large number of points need to be calculated. Also, the uncertainties in the aircraft model makes the prediction of the stability boundary difficult (Chung et al., 2008).
In modern development, authors have introduced robust flutter analysis including structural and aerodynamics uncertainties. Initially, Lind and Brenner (1999) introduced structured singular value 
into the aeroelastic field and developed a -method in which the aeroelastic system is parametrized with dynamic pressure perturbations. Chung et al. (2008) included aerodynamic uncertainties by considering variation of Mach number, which is represented by the boundary of Theodorsen's lift deficiency function. Chung and Shin (2010) defined the worst case flutter boundary by considering the structural variation due to changing natural frequency and aerodynamic variation and, according to the authors, this approach is more conservative than a nominal case.
In this context, this paper proposes an alternative approach to include uncertainties in the aeroelastic stability problems. Two applications are introduced at which the problem is formulated as a polytopic differential inclusion system considering uncertainties in the air density (or its equivalent altitude) or in the airspeed. The main idea is to identify regions of stability in the flight envelope by solving a linear matrix inequality. The method is written in time domain using the state space realization and the Theodorsen's aerodynamic model using a by rational function approximation for the aerodynamic forces. This methodology offers promise for robust flutter analysis using convex optimization.
Nomenclature
s = Laplace variable
M = mass matrix
D = damping matrix
K = stiffness matrix
Q = aerodynamic matrix
u = vector of displacement
A = dynamic matrix
C = output matrix
x = vector of states
y = vector of output
b = semi chord
v = reduced frequency
Nlag = number of lags
v= number of vertices
V= airspeed
VLyap = Lyapunov energy function
q = dynamic pressure
mM= Mach number
M = airfoil mass
Greek Symbols
Subscripts
Superscripts
Aeroelastic Model
The aeroelastic model is characterized by the mass, stiffness, damping and aerodynamic matrices. These matrices are typically obtained from finite element method and the system is commonly represented in frequency domain using modal coordinates. In this case, the equation of motion is written as
where the subscript m indicates the modal coordinates. Aerodynamic matrices can be obtained using the Doublet Lattice method in a constant Mach number (mM) for each reduced frequency k (Albano and Rodden, 1969). However, in this work they are computed using the Theodorsen's aerodynamic theory. The method proposed requires a representation in time domain of the system, in Eq. (1), but because the term Qm(mM,K) has no Laplace inverse, it is necessary to write it in terms of rational function approximation, so that the aerodynamic forces can be transformed to time domain. This is done using Roger's approximation (Roger, 1977), given by Eq. (2). Its contains a polynomial part representing the forces on the section acting directly connected to the displacements u(t) and their first and second derivatives. Also, this equation has a rational part representing the influence of the wake acting on the section with a time delay.
where s is the Laplace variable, nlag is the number of lag terms and 
is the jth lag parameter (j = 1. ..., nla).
Substituting Eq. (2) into Eq. (1), it is possible to obtain, after some rearrangements, the classical aeroelastic equation of motion in a state space format as:
where 
is the state vector and 
are states of lags required for the approximation. The output matrix 
has dimension 
, where 
and 
, respectively the velocity and displacement, are output vectors. Matrix A is presented in the following form:
where
and matrices 
, for 
, are defined in the appendix.
In classical solutions the aeroelastic stability is treated as an eigenvalue problem of the dynamic matrix A. Frequency of oscillation and damping are extracted from the complex eigenvalues. This analysis is performed for fixed values of velocity and air density 
.
In Eq. (4), matrix A has a term 
that has to be inverted. For mass normalized eigenvectors, 
(identity matrix), so the inverse of the aeroelastic mass matrix can be rewritten as
For the cases where 
vanishes when j is large, the Taylor expansion can be used as an approximation of Eq. (6) and then, it can be written as
where NT is a positive natural number chosen conveniently.
For the problem considered in this paper, the assumption 
when 
has shown to be true. When this assumption is not satisfied, it is possible to use another approximation for the inverse of 
. Different methods to compute the inverse of sum of matrices can be found in Petersen and Pedersen, 2008, and Henderson and Searle, 1981.
When this assumption is satisfied, the approximated state space matrix of Eq. (3) can be written as
where the submatrices 
, are given in the appendix.
Continuous Range Stability Analysis
The objective of a continuous stability analysis is to introduce a robust procedure for checking if flutter is present in a range of a varying parameter. In the classical method, the discrete analysis is performed by reducing the intervals of a varying parameter, such as velocity or air density, which can be costly for analysis of a full flight envelope. Two procedures are presented here, the first considers a fixed velocity and analysis is carried out for a range of air density. The second procedure is formulated for fixed value of air density and varying parameter is now velocity.
Procedure 1: Analysis for air density range
Considering the aeroelastic model written in state space format to describe a range of altitude for a fixed value of velocity, as illustrated in Fig. 1 . A convex space is obtained using a polytopic differential inclusion system (PLDI) and the air density given by
where the subscript 
indicates the uncertain air density 
and 
denotes its associated uncertainty. In this case, the aeroelastic model is written as a PLDI system and the dynamic matrix is given by
and, if 
, then
where 
is the j-th uncertain parameter related to that polytopic system.
Procedure 2: Analysis for a velocity range
Considering the system of Eq. (3) to be written in state space format as a PLDI system and its convex space representing an airspeed range as shown in Fig. 2, the uncertainty range can be written as
where the subscript indicates the uncertain airspeed 
and 
denotes its associated range in terms of percentage. In this case, the aeroelastic model is written as a PLDI system and the dynamic matrix is given by
and, if 
, then
where 
and 
are understood as uncertain parameters related to that polytopic system and physically they represent the airspeed range. In this case the maximum and minimum values of 
are given by
Substituting 
by 
and 
by 
into Eq. (4) it is possible to define the dynamic matrix that describes the PLDI system and represents the velocity range 
. In this case, the convex space is geometrically represented by a plane.
Lyapunov Energy Function
Consider the energy function 
. The aeroelastic system defined in a particular pair air density and velocity is stable if
Equation (16) represents the energy dissipation over time, idea which has been proposed by Alexander Lyapunov (in 1907). The method is traditionally applied to the stability analysis in control theory because Lyapunov showed that the differential equation (Eq. (3)) is stable if and only if exists a positive-defined P such that (Boyd et al., 1994)
Equation (17) is called Lyapunov inequality on P and its first application in practical control engineering problems was introduced in the early 1960's (Boyd et al., 1994).
Quadratic stability criteria for aeroelastic analysis
Consider a PLDI system given by
where the dynamic matrix 
belongs to a convex polytopic set defined as
where 
and 
(Oliveira et al., 1999). Also, 
is the number of vertices of the polytopic system and 
denotes its jth vertex. The number of vertices is given by 
for application 1 (and 22 = 4 for application 2) and the operator Co means that matrices 
define the desired convex space 
(Boyd et al., 1994).
Substituting Eq. (19) into inequality (17), the Lyapunov inequality is rewritten as
Pre and post multiplying the last LMI by 
and after some rearrangement, the quadratic stability criteria is written as
The necessary condition to ensure the quadratic stability of the aeroelastic system into this convex space is to solve Eq. (21) at its vertices for which 
with 
and 
. In this case, all LMIs must be satisfied simultaneously (Barmish, 1985).
The vertices of the parameter box (or the convex space) are a combination of the minimum and maximum values of the parameters of the system. It is supposed that the system can assume any combination of values inside the box (Boyd et al., 1994; Silva and Lopes Jr., 2006).
Lyapunov Stability Index
In order to identify regions into the flight envelope at which system is stable, the Lyapunov stability index is defined to classify the regions of aeroelastic stabilities. Particular details are presented here and the conceptual understanding is the same for both applications.
Application 1: Robust air-density analysis
The convex space 
must be understood as a fictitious domain at which the aeroelastic stability is verified. The main reason is that although the product 
represents a physical parameter variation, the parameters 
, are introduced only for creating the convex space and they have no physical information. However, the physical domain 
and then:
• if there is X such that all inequalities in (22) are simultaneously feasible 
the system is stable in 
and stable in 
. So,
• otherwise the system is not stable in . So,
where 
. If the proposed Lyapunov stability index is 
, it is not possible to say that the system is unstable in 
because it can be stable in that physical domain, but unstable in 
. In this case a complementary discrete stability analysis must be performed in each point 
and 
.
Application 2: Robust velocity analysis
The convex space 
is a fictitious domain at which the aeroelastic stability is verified. Similarly to application 1, the parameter 
is introduced as a mathematical artifice only for creating the convex space. In this case, the physical domain 
and then:
• if there is X such that all inequalities in Eq. (22) are simultaneously feasible the system is stable in and then stable in 
. So,
• otherwise the system is not stable in . So,
where 
. If the proposed Lyapunov stability index is 
, it is not possible to say that the system is unstable in 
because it can be stable in that physical domain, but unstable in 
. In this case a complementary discrete stability analysis must be performed in each point 
and 
.
Numerical Application
In order to illustrate the approach, numerical simulations are performed on a three degree of freedom airfoil section (semichord b), for which the equations of motion describing an aeroelastic response are presented by Theodorsen (1935). An illustrative scheme is presented in Fig. 3 . The structural mass and stiffness and the aerodynamic forces matrices are not presented in this paper because this classical problem is easily found in literature. The physical displacement vector is 
, where 
and 
are the degrees of freedom of plunge, pitch and control surface rotation, respectively. Its physical and geometric properties are presented in Table 1 (note that 
).
Considering 
and 
, let 
and 
to 
be the flight envelope defined for this work. Figures and present frequency and damping plotted against airspeed for nominal system (). Additionally, flutter speeds are identified as 17.4 m/s and 17.0 m/s, respectively for and 
extracting the eigenvalues of the dynamic matrix. These reference results are used to show that this partially continuous method is able to identify stability regions in the flight envelope.
Application 1: Robust air density analysis
The stability analysis is performed according to the first approach (application 1) increasing 
in 0.1 m/s and also considering 
. In the example 
the Taylor's expansion provides the inverse aeroelastic mass matrix as shown below
where its possible to note a good approximation for the inverse matrix.
The convex space 
is presented in Fig.4 for which 
, 
and the dynamic matrices are defined by
Note that the matrix 
defines the rjth vertex based on 
and 
, such that
Lyapunov inequality is solved for the range of interest of air density and the stability index is plotted against airspeed (Fig. 7). According to these results the system is stable for V m/s. Although a classical eigenvalue analysis indicates that the system is stable if 
, only this introduced approach is able to indicate stability over all continuous interval 
for each velocity. In order to verify the system stability in this flight envelope, only a complementary discrete analysis for each air density must be performed if 
. Figure illustrates the geometric region of analysis of the flight envelope for this novel approach. This continuous analysis over to 
is an advantage mainly for flutter analysis in complex and large structures whose stability boundaries may be not easily predicted through discrete evaluation.
Application 2: Robust velocity analysis
The stability analysis is performed according to the second approach (application 2) particularly for the nominal air density and for each range 
. Different values of 
are considered in order to demonstrate that the approach is able to identify stability regions in the flight envelope. Figures 9 through 15 present the stability indexes for this partially continuous stability analysis; and Fig. 16 is obtained by performing a discrete analysis 
. If 
is large, a small number of analyses is required for evaluating the range 
. However, the aeroelastic stability region cannot be completely identified. Note that the mathematical domain can result in not identifying the stability even at low speeds. This work does not propose criteria which should be used for choosing the optimal values of .
Conclusions
This paper has presented the aeroelastic stability analysis using a polytopic differential inclusion system. Polytopic systems have been used extensively in control design; however, they are not commonly used for solving aeroelastic problems. An eigenvalue analysis is defined as discrete method over the air density and it is used to compare the results with this introduced methodology. Simulations on a simple system have shown that this approach can be used to provide robustness in the flutter identification. Also, the paper has demonstrated that Lyapunov's linear matrix inequality is a sufficient condition for verifying the aeroelastic stability.
Acknowledgements
The authors gratefully acknowledge CNPq for partially funding the present research work through the INCT-EIE.
Paper received 1 May 2012
Paper accepted 22 July 2012.
Appendix
This appendix presents the submatrices introduced in Eq. (8).
- Albano, E., Rodden, W.P., 1969, ``A Doublet-Lattice Method for Calculating Lift Distributions on Oscillationg Surfaces in Subsonic Flows,'' AIAA Journal, pp. 279-285.
- Barmish, B.R., 1985, ``Necessary and sufficient conditions for quadratic stabilizability of an uncertain system,'' Journal of Optimization Theory and Applications, Vol. 4, pp. 399-408.
- Bisplinghoff, H.A., Raymond, L., Halfman, R.L., 1996, ``Aeroelasticity,'', Dover Publications, INC.
- Boyd, S., Ghaoui, L.E., Balakrishnan, V., 1994, ``Linear Matrix Inequalities in System Control Theory,'' Society for Industrial and Applied Mathematics, Vol. 15.
- Chen, P.C., 2000, ``Damping Perturbation Method for Flutter Solution: The g-Method,'' AIAA Journal, Vol. 38, pp. 1519-1524.
- Chung, C., Shin, S., 2010, ``Worst case flutter analysis of a stored wing with structural and aerodynamic variation,'' AIAA Structures, Structural Dynamics, and Materials Conference, Orlando, Florida.
- Chung, C., Shin, S., Kim, T., 2008, ``A new robust aeroelastic analysis including aerodynamic uncertainty from varying Mach numbers,'' AIAA Structures, Structural Dynamics, and Materials Conference, Schaumburg, IL.
- Hassig, H.J., 1971, ``An Approximate True Damping Solution of the Flutter Equation by Determinant Iteration,'' Journal of Aircraft, Vol. 8(11), pp. 885-889.
- Henderson, H.V. and Searle, S.R., 1981, ``On Deriving the Inverse of a Sum of Matrices,'' Society for Industrial and Applied Mathematics, Vol. 23, No. 1, pp. 53-60.
- Lind, R., Brenner, M., 1999, ``Robust aeroservoelastic stability analysis,'' Spring-Verlag, London, April.
- Lyapunov, A., 1907, ``Probleme general de la stabilite du mouvement,'' Annales Fac. Sciences Toulouse, Vol. 7, pp. 203-474.
- Oliveira, M.C., Bernussou, J., Geromel, J.C., 1999, ``A new discrete-time robust stability condition,'' Systems & Control Letters, Vol. 37, pp. 261-265.
- Pahlavanloo, P., 2007, ``Dynamic Aeroelastic Simulation of the AGARD 445.6 Wing using Edge,'' FOI - Swedish Defence Research Agency Defence and Security, System and Technology
- Petersen, K. B. and Pedersen, M.S., 2008, ``The Matrix Cookbook,'' http://matrixcookbook.com, Version November 14.
- Roger, K., 1977, ``Airplane Math Modelling Methods for Active Control Design,'' Structural Aspects of Control, AGARD Conference Proceeding, Vol. 9, pp. 4.1-4.11.
- Silva, S., Lopes Jr., V., 2006, ``Active Flutter Suppression in a 2-D Airfoil Using Linear Matrix Inequalities Techniques,'' Journal of the Brazilian Society of Mechanical Sciences & Engineering, Vol. XXVIII, pp. 84-93.
- Theodorsen, T., 1935, ``General theory of aerodynamic instability and the mechanism of flutter,'' Technical report 496, NACA -- National Advisory Committee for Aeronautics.
- Theodorsen, T., Garrick, I.E., 1940, ``Mechanism of flutter A theoretical and experimental investigation of the flutter problem,'' NACA -- National Advisory Committee for Aeronautics.
- Yates, E.C., 1988, ``AGARD Standard Aeroelastic Configurations for Dynamic Response I-wing 445.6,'' Advisory Group for Aerospace Research and Development.
Aeroelastic stability analysis using linear matrix inequalities
Publication Dates
-
Publication in this collection
24 July 2013 -
Date of issue
2012
History
-
Received
01 May 2012 -
Accepted
22 July 2012
