SciELO - Scientific Electronic Library Online

vol.34 special issue 2Special Issue 2: Uncertainties 2012Effect of parametric uncertainties on the performance of a piezoelectric energy harvesting device author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Journal of the Brazilian Society of Mechanical Sciences and Engineering

Print version ISSN 1678-5878

J. Braz. Soc. Mech. Sci. & Eng. vol.34 no.spe2 Rio de Janeiro  2012 



Aeroelastic stability analysis using linear matrix inequalities



Douglas D. BuenoI; Clayton R. MarquiII; Luiz C. S. GóesIII; Paulo J. P. GonçalvezIV
IIITechnological Institute of Aeronautics - São José dos Campos - 12228-900 SP - Brazil -
IVUniversidade Estadual Paulista - Bauru - 17033-360 SP - Brazil




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




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.



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




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:


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 .




















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.



The authors gratefully acknowledge CNPq for partially funding the present research work through the INCT-EIE.



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.         [ Links ]

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.         [ Links ]

Bisplinghoff, H.A., Raymond, L., Halfman, R.L., 1996, ``Aeroelasticity,'', Dover Publications, INC.         [ Links ]

Boyd, S., Ghaoui, L.E., Balakrishnan, V., 1994, ``Linear Matrix Inequalities in System Control Theory,'' Society for Industrial and Applied Mathematics, Vol. 15.         [ Links ]

Chen, P.C., 2000, ``Damping Perturbation Method for Flutter Solution: The g-Method,'' AIAA Journal, Vol. 38, pp. 1519-1524.         [ Links ]

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.         [ Links ]

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.         [ Links ]

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.         [ Links ]

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.         [ Links ]

Lind, R., Brenner, M., 1999, ``Robust aeroservoelastic stability analysis,'' Spring-Verlag, London, April.         [ Links ]

Lyapunov, A., 1907, ``Probleme general de la stabilite du mouvement,'' Annales Fac. Sciences Toulouse, Vol. 7, pp. 203-474.         [ Links ]

Oliveira, M.C., Bernussou, J., Geromel, J.C., 1999, ``A new discrete-time robust stability condition,'' Systems & Control Letters, Vol. 37, pp. 261-265.         [ Links ]

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.         [ Links ]

Petersen, K. B. and Pedersen, M.S., 2008, ``The Matrix Cookbook,'', Version November 14.         [ Links ]

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.         [ Links ]

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.         [ Links ]

Theodorsen, T., 1935, ``General theory of aerodynamic instability and the mechanism of flutter,'' Technical report 496, NACA -- National Advisory Committee for Aeronautics.         [ Links ]

Theodorsen, T., Garrick, I.E., 1940, ``Mechanism of flutter A theoretical and experimental investigation of the flutter problem,'' NACA -- National Advisory Committee for Aeronautics.         [ Links ]

Yates, E.C., 1988, ``AGARD Standard Aeroelastic Configurations for Dynamic Response I-wing 445.6,'' Advisory Group for Aerospace Research and Development.         [ Links ]



Paper received 1 May 2012
Paper accepted 22 July 2012.




This appendix presents the submatrices introduced in Eq. (8).

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License