Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2013 / Article

Research Article | Open Access

Volume 2013 |Article ID 941689 |

Douglas Domingues Bueno, Clayton Rodrigo Marqui, Luiz Carlos Sandoval Góes, Paulo José Paupitz Gonçalves, "The Use of Gramian Matrices for Aeroelastic Stability Analysis", Mathematical Problems in Engineering, vol. 2013, Article ID 941689, 9 pages, 2013.

The Use of Gramian Matrices for Aeroelastic Stability Analysis

Academic Editor: Cristian Toma
Received18 Nov 2012
Revised01 Mar 2013
Accepted01 Mar 2013
Published10 Apr 2013


Most of the established procedures for analysis of aeroelastic flutter in the development of aircraft are based on frequency domain methods. Proposing new methodologies in this field is always a challenge, because the new methods need to be validated by many experimental procedures. With the interest for new flight control systems and nonlinear behavior of aeroelastic structures, other strategies may be necessary to complete the analysis of such systems. If the aeroelastic model can be written in time domain, using state-space formulation, for instance, then many of the tools used in stability analysis of dynamic systems may be used to help providing an insight into the aeroelastic phenomenon. In this respect, this paper presents a discussion on the use of Gramian matrices to determine conditions of aeroelastic flutter. The main goal of this work is to introduce how observability gramian matrix can be used to identify the system instability. To explain the approach, the theory is outlined and simulations are carried out on two benchmark problems. Results are compared with classical methods to validate the approach and a reduction of computational time is obtained for the second example.

1. Introduction

Between various physical phenomena involving fluid-structure interaction, flutter is probably the most representative topic studied in engineering applications such as aircrafts and bridges. The flutter phenomenon is an interaction between structural dynamics and aerodynamics that results in divergent and destructive oscillations of motion [1].

In 1935, Theodorsen [2] proposed a method of flutter analysis in a discrete system by including aerodynamic forces in frequency domain and formulating the analysis as a complex eigenvalue problem. Hassig [3] proposes the pk-method where the unsteady aerodynamic matrix is represented by a function of the complex eigenvalues. Using an iterative algorithm the value of a reduced frequency converges to the imaginary part of a system eigenvalue. Chen [4] also proposes a flutter method including a first-order damping term into the equation of motion known as the g-method. According to the author, this method generalizes the -methods and pk-methods for reliable damping prediction and has proved to be efficient in obtaining unlimited number of aerodynamic lag roots.

These methodologies, which are well established in the research and engineering community, were developed decades ago and have been used in the development of almost all flying commercial and military aircraft.

In this context, this paper proposes an alternative approach for detecting flutter using observability Gramian matrices. The proposed methodology is developed in time domain using state-space representation of the aeroelastic system. The elements of a Gramian matrix are related to the energy of vibration modes and can be seen as an improved observability matrix, introduced by Kalman et al. [5].

Gramian matrices have been used in the field of control engineering. Their fundamental concepts were proposed by Moore after introducing the balanced reduction for state-space models [6] and have been used for applications such as optimal placement of sensors and actuators [79]. After this work the approach was extended for unsteady aerodynamic and aeroelastic systems [1013]. According to [14], despite these efforts, Gramian matrices are still neglected by the scientific community.

The principal objective of the paper is to show that observability Gramian matrices can be used to detect flutter. The rationale for using this approach is that Gramian matrices contain, in a single index, information about the energy that is transferred from the flux to the structure and then dissipated by any damping mechanism for each flight condition considered in the analysis (flight envelope).

It is shown that the Gramians are sensitive to flutter where energy transferred from the flux in a cycle is larger than the energy dissipated by the damping mechanism.

Part of the procedure to determine the Gramian matrices requires a system defined in time domain, including the aerodynamic forces. If aerodynamic forces are written in terms of reduced frequencies, this can be done using one of the methods such as least square [15], minimum state [16] and the mixed state [8].

The paper illustrate’s the process to obtain the Gramian matrix from an aeroelastic system in time domain. Two numerical examples are used to compare the proposed method with the methods in the literature. The first example is three degrees of freedom typical section airfoil and the second is the AGARD 445.6 wing model developed using finite element method [17, 18].

The paper shows that it is possible to obtain some advantage in terms of computation time using the proposed methodology.

2. Aeroelastic Model for Stability Analysis

Assuming a general aeroelastic model that can be written in the state-space form, according to where is the state vector, is the dynamic matrix, is the input matrix, is a vector of external forces, is the known as the output matrix, and is the output vector. Using this equation the aeroelastic system can be described by different aerodynamic theories and the system of equations can be written in physical or modal coordinate systems (complementary information is presented in Appendix A). The time-invariant aeroelastic system defined by matrix is stable for a range of velocities defined by increasing airspeed values containing a flutter airspeed , if .

That stability could be verified by solving an eigenvalue problem for each discrete airspeed point in the flight envelope and checking if the real part of system eigenvalues is negative values. This can be time consuming, specifically for large dimension systems. To overcome this process, the observability Gramian method presented in this paper is based on the solution of a set of linear equations. The next section presents the bases to write the problem in appropriate format.

3. Observability and Gramian Matrices

The concept of observability involves the dynamic matrix and the output matrix . A linear system, or the pair , is observable at instant of time , if the state can be determined from the output , with , where is a finite instant of time. If this is true for all initial time and all initial states , the system is said to be completely observable [14].

A linear time-invariant system with outputs is said to be completely observable, if and only if the observability matrix with dimension has hank . With the observability matrix given by where is the dimension of matrix . This concept can be easily implemented to verify system observability but may lead to numerical overflow for systems represented by large dimension matrices [14]. Also only qualitative information about the system is provided (see [14, 19, 20]).

An alternative approach that can be applied for large-order problems is to use the observability Gramian matrix. The observability Gramian matrix is defined to express quantitative properties of the system, considering it at time written as which, according to [14], can be determinated as where indicates a time-variant property. If a linear time-invariant and stable system is considered, then the observability Gramian matrix can be computed using the Lyapunov equation [14]

An important property between observability and Gramian matrices is that they share the same Kernel (i.e., the set of all vectors for which ),

According to [21, 22], one consequence of this property is that the energy detected by an output state can be computed through the observability Gramian matrix. This is done by writing an expression for the energy detected (or observed) by the output at time caused by the system's initial state such that

Equation (5) is only defined for a stable system [14]. To include the flutter speed, it is necessary to modify the equation using the generalized ordinary cross-Gramian matrix, , introduced by Zhou et al. [23], defined for both stable and unstable systems. Then, let be the solution to the Riccati equation The observability Gramian matrix is a submatrix of which can be computed by solving the following equation where and the modified output matrix is used instead of . This modified output matrix with dimension is defined such that where the th row satisfies the equation , and the other rows are filled with zeros. Thus, considering and the aeroelastic matrix defined at the airspeed , the observability Gramian matrix is computed by solving (9) for the pair . The input matrix is written by to represent inputs.

3.1. Complex Schur Decomposition

Using a complex Schur decomposition, the Lyapunov equation can be reshaped as a real linear system of equation and solving it, is obtained [24] where is an identity matrix with appropriate dimension, vec makes a column vector out of a matrix by stacking its columns, and indicates a Kronecker product [25].

3.2. Gramian Paramater to Detect Flutter

The hypothesis introduced in this paper is that observability Gramian matrix contains information which can be used to indicate the amount of energy transferred from the air flow to the structure. In this case, this amount of energy is maximum at the airspeed at which the system becomes unstable.

In order to prove this hypothesis, a Gramian parameter is defined and obtained by computing a matrix norm of . This parameter indicates two main features: the contribution of each aeroelastic mode absorbing energy from the air flow and the airspeed where flutter occurs.

By computing for each pair and it is possible to build a matrix (14). Assuming that aeroelastic instability occurs in this speed range, the flutter speed is found for the largest value of : where is the length of the airspeed vector

The th Gramian parameter in each column of the matrix containing is related to a measure of the energy absorbed by the th aeroelastic mode. The values of in each column can then be compared to determine the mode contribution. On the other hand, the row of which contains can be plotted to determine the airspeed of flutter. (This is illustrated in Figure 1.)

3.3. Matrix Norm

In practical implementations the Gramian parameter can be obtained by different matrix norms. Matrix norms are often used to provide quantitative information. In this paper, Frobenius norm is used (15). However, similar results were obtained using other norms presented in Appendix B, where is the Gramian matrix element of the th row and th column, and .

4. Numerical Simulations

4.1. Typical Section Airfoil

The effectiveness of the approach was determined through simulations using two examples. The first was a three-degree-of-freedom typical airfoil section (semichord ). The equations of motion describing the aeroelastic response are presented in [17]. The bidimensional system was formulated in modal coordinate system considering the plunge and pitch modes and the control surface deflection to represent the third mode. Its physical and geometric properties are presented in Table 1 and an illustration of the airfoil is shown in Figure 2.

Aerodynamic semichord  m
Airfoil mass  kg
Plunge frequency  Hz
Pitch frequency  Hz
Control surface rotation frequency  Hz
Air density
Lag parameters (Roger's method)
Reduced frequencies
Figure 2
Figure 2
Distance between c.e. to c.g.
Distance between c.e. to c.g. (flap)
Radius of gyration of the flap referred to
Radius of gyration of the airfoil referred to
Elastic center c.e.
Center of gravity c.g.

The results of the proposed approach were compared with classical eigenvalues analysis. Figures 3 and 4 present aeroelastic frequency and damping, both plotted as a function of the airspeed. According to the literature, although the plunge (or bending) mode is stable, the system becomes unstable because the pitch (or torsion mode) is unstable. It is said that airfoil normally is undergoing a flutter oscillation composed of important contributions of these modes. In this example, flutter airspeed was computed equal to  m/s.

The results of the proposed method are presented in Figure 5. The flutter speed was correctly identified. Similar results were also obtained using the other norms and therefore they are not shown herein.

4.2. AGARD 445.6 Wing

The AGARD 445.6 benchmark wing was also used to demonstrate the method. The structural model for the AGARD 445.6 wing was developed using finite element method using the MSC/NASTRAN. The finite element model consisted of plate elements with single-layer orthotropic material. The model has 231 nodes and 200 elements. Rotation of nodes was neglected, allowing three degrees of freedom per node. See physical properties in Appendix C and complementary details in [18, 26].

Aerodynamic and structural matrices were obtained from MSC/NASTRAN program (Aeroelastic Solution 145) for the Mach number 0.50 and  kg/m3. The values of reduced frequencies are presented in Appendix C. The model has a reference length of  m, a sweep angle of 45 degrees at the quarter chord line, a semispan of 0.762 m and a taper ratio of 0.66. The flutter boundary is investigated only using the first four fundamental structural modes and their natural frequencies are, respectively,  Hz,  Hz,  Hz, and  Hz. Details can be found in literature (e.g., see [18, 26]).

A state-space model was obtained through a MATLAB implementation for which the parameters of lag were chosen as , , , and . The observability Gramian matrix was computed for each pair and was computed using different matrix norms. Analysis was evaluated considering pairs of with constant Mach number. Using the traditional pk-method, the flutter speed was identified at  m/s. The classical and diagrams are shown in Figures 6 and 7.

With the proposed methodology, the Gramian parameter is showed in Figures 8 and 9. Based on Figure 8, it is possible to note that the first and second aeroelastic modes contribute more to the flutter than the third and fourth modes. This information can be confirmed in Figure 7, that shows small changes of aeroelastic damping for these modes.

4.3. Comparison between pk-Method and the Proposed Approach

This section shows that it is possible to reduce the computational cost to find the flutter speed using Gramian matrices, especially for industrial applications where the number of nominal and parametric cases of analysis can be very large. Based on the second example (AGARD wing), it is demonstrated that (5) and (9) can be conveniently used to get small computational time in comparison with the pk-method.

According to previous sections, (5) is not valid for unstable systems; that is, the physics is represented if and only if the system is stable. Then, if this equation is used to compute the Gramian parameter, the solution has no physical meaning after the flutter velocity (including that point). However, the computational time to solve the linear system of equations (see Section 3.1) is smaller than that to solve the iterative eigenvalue algorithm (pk-method). The following results were obtained using the personal computer with operational system MS Windows 7 (64 bits), Intel(R) Core(TM) i5 CPU M460 2.53 GHz, and RAM 4.0 GHz.

Figure 10 shows that Gramian parameters computed from the cross-Gramian matrices have computational cost larger than classical pk-method. However, if they are computed from Gramian matrices (5), the computational time is smaller. In this context, it is possible to use the following procedure.(1)Compute Gramian parameters using observability Gramian matrices obtained from (5) for every points into the flight envelope.(2)Identify the point at which Gramian parameter is maximum, according to (13).(3)Compute cross-Gramian parameters for points around the maximum.

Using this procedure, it is possible to obtain some reduction in the computational time to evaluate the system stability. This reduction is represented in Figure 11.

5. Conclusions

Controllability and observability Gramian matrices have been used extensively in control design and for optimal placement of sensors and actuators in smart structures. They have also been used for solving aeroelastic problems mainly for writing reduced-order models.

This paper has investigated the applicability of observability Gramian matrix to detect the flutter speed in two benchmark structures. It was shown that Gramian matrices represent the transfer of energy from the flow to the structure. That transfer of energy increases when the airspeed is close to the flutter condition. A classical method was used to compare and validate the results obtained by Gramian matrices.

This approach does not compute frequency and damping for the aeroelastic modes. However it allows to determine the aeroelastic modes with largest contribution to the flutter mechanism and the airspeed where the system becomes unstable. Additionally, one advantage is related to the reduction of computational time for analysis when compared to classical pk-method.

This work is a complementary effort to apply the concepts from control theory in problems involving aeroelasticity in time domain. It can also be used in flight flutter tests using an identification method to obtain the state-space representation instead of identifying aeroelastic frequencies and damping.


A. State-Space Represenation of an Aeroelastic System

An aeroelastic system consisting of structural parameters (mass, damping, and stiffness) subject to the forces from the fluid can be modeled in the Laplace domain using a physical or modal coordinate system. In general, the modal coordinates are used to truncate the system of equations using the eigenvector extracted from structural mass and stiffness matrices. Without loss of generality, both numerical case studies were performed using modal coordinates as shown in the following equations.

The matrices representing the structural parameters are the mass , the viscous damping , and the stiffness . In most cases, these matrices are obtained from finite element method and are reduced to modal coordinates (represented by the subscript ); the system displacement vector in modal coordinates is . The aerodynamic influence matrix depends on the parameters (reduced frequency) and (Mach Number) and is the dynamic pressure caused by the flowing fluid ( can be obtained by methods such as Doublet Lattice). Note that the proposed approach is not restricted to this aerodynamic theory: where is the Laplace variable.

In this case, the problem that arises from the conversion of (A.1) to time domain is the fact that has no Laplace inverse. This is overcome by writing a rational function to represent the aerodynamic coefficients in time domain. In this work, Roger's approximation [15] is used, containing a polynomial part representing the forces acting directly related to the displacements and their first and second derivatives. Also, this equation has a rational part representing the influence of the wake acting on the structure with a time delay: where is the number of lag terms and is the th lag parameter . The parameters were chosen arbitrarily to ensure equality of (A.2) for each reduced frequency and, then, their values are different for each system presented herein. The procedure to obtain the aerodynamic matrices is discussed in [15].

Equation (A.2) is substituted into (A.1) allowing to write the aeroelastic system in time domain using a state-space representation. In this case, the state vector presented in (1) is , where are states of lags required for the approximation. The dynamic matrix is given bywhere the matrices , , and are comprised of terms containing structural and aerodynamic coefficients which are given by

B. Norms for Computing the Gramian Parameter

The matrix norms for computing the Gramian parameters are presented as follows.

B.1. -Norm

The -norm of   is the maximum of the row sums; that is,

B.2. 2-Norm

The 2-norm is computed based on the largest singular value of ,   ; that is, where denotes a Hermitian matrix.

B.3. 1-Norm

The 1-norm is the maximum of the columns sums; that is,

C. AGARD 445.6 Wind Structural Model

This appendix presents complementary information for the AGARD 445.6 wind and more details can be found in [18, 26].

Figures 12, 13, 14, and 15 show the four structural modes considered in the finite element model obtained from MSC/NASTRAN program.

The thickness distribution is governed by the airfoil shape. The material properties used are ,  GPa, ,  GPa, and kg/m3, where and are the moduli of elasticity in the longitudinal and lateral directions, is Poisson's ratio, is the shear modulus in each plane, and is the wing density.

There are 10 elements in chordwise and 20 elements in spanwise direction. There are a total of 231 nodes and 200 elements. The boundary conditions of the wing are selected in accordance with the physical model. The root is cantilevered except for the nodes at the leading and trailing edges. Additionally, the rotation around all axis degrees of freedom at all nodes is constrained to zero. In this case, there are 666 degrees of freedom in the model in physical coordinates. This work considered four modes to obtain the model in modal coordinates system.

Table 2 shows the reduced frequencies used to compute the aerodynamic matrices for the second case study (AGARD wing 445.6).

Conflict of Interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interests.


  1. H. A. Bisplingho, L. Raymond, and R. L. Halfman, Aeroelasticity, Dover, 1996.
  2. T. Theodorsen, “General theory of aerodynamic instability and the mechanism of flutter,” Tech. Rep. 496, NACA National Advisory Committee for Aeronautics, 1935. View at: Google Scholar
  3. H. J. Hassig, “An approximate true damping solution of the flutter equation by determinant iteration,” Journal of Aircraft, vol. 8, no. 11, pp. 885–889, 1971. View at: Google Scholar
  4. P. C. Chen, “Damping perturbation method for flutter solution: the g-method,” The American Institute of Aeronautics and Astronautics Journal, vol. 38, no. 9, pp. 1519–1524, 2000. View at: Google Scholar
  5. R. E. Kalman, Y. C. Ho, and K. S. Narenda, “Controllability of linear dynamics system,” Contribution to Dierential Equations, vol. 1, no. 2, pp. 189–213, 1962. View at: Google Scholar
  6. B. C. Moore, “Principal component analysis in linear systems: controllability, observability, and model reduction,” IEEE Transactions on Automatic Control, vol. 26, no. 1, pp. 17–32, 1981. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  7. A. Arbel, “Controllability measures and actuator placement in oscillatory systems,” International Journal of Control, vol. 33, no. 3, pp. 565–574, 1981. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  8. D. Biskri, R. M. Botez, N. Stathopoulos, S. Thérien, A. Rathe, and M. Dickinson, “New mixed method for unsteady aerodynamic force approximations for aeroservoelasticity studies,” Journal of Aircraft, vol. 43, no. 5, pp. 1538–1542, 2006. View at: Publisher Site | Google Scholar
  9. G. Schulz and G. Heimbold, “Dislocated actuator/sensor positioning and feedback design for flexible structures,” Journal of Guidance, Control, and Dynamics, vol. 6, no. 5, pp. 361–367, 1983. View at: Publisher Site | Google Scholar
  10. M. L. Baker, “Approximate subspace iteration for constructing internally balanced reduced order models of unsteady aerodynamic systems,” in Proceedings of the 37th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference and Exhibit, Salt Lake City, Utah, USA, April 1996, Technical Papers. Part 2 (A96-26801 06-39). View at: Google Scholar
  11. K. Willcox and J. Peraire, “Balanced model reduction via the proper orthogonal decomposition,” The American Institute of Aeronautics and Astronautics Journal, vol. 40, no. 11, pp. 2323–2330, 2002. View at: Google Scholar
  12. D. J. Lucia, P. S. Beran, and W. A. Silva, “Reduced-order modeling: new approaches for computational physics,” Progress in Aerospace Sciences, vol. 40, no. 1-2, pp. 51–117, 2004. View at: Publisher Site | Google Scholar
  13. M. J. Balas, “Reduced order modeling for aero-elastic simulation,” AFOSR Final Report AFRL-SR-AR-TR 06-0456, Air Force Office of Scientiffic Research, Arlington, Va, USA, 2006. View at: Google Scholar
  14. W. K. Gawronski, Dynamics and Control of Structures: A Modal Approach, Springer, New York, NY, USA, 1st edition, 1998. View at: MathSciNet
  15. K. Roger, “Airplane math modelling methods for active control design,” in Structural Aspectsof Control, vol. 9 of AGARD Conference Proceeding, pp. 4.1–4.11, 1977. View at: Google Scholar
  16. M. Karpel, “Design for active and passive flutter suppression and gust alleviation,” Tech. Rep. 3482, National Aeronautics and Space Administration—NASA, 1981. View at: Google Scholar
  17. T. Theodorsen and I. E. Garrick, Mechanism of Flutter a Theoretical and Experimental Investigation of the Flutter Problem, NACA National Advisory Committee for Aeronautics, 1940.
  18. E. C. Yates, AGARD Standard Aeroelastic Congurations for Dynamic Response I-Wing 445.6, Advisory Group for Aerospace Research and Development, 1988.
  19. L. A. Aguirre, “Controllability and observability of linear systems: some noninvariant aspects,” IEEE Transactions on Education, vol. 38, no. 1, pp. 33–39, 1995. View at: Publisher Site | Google Scholar
  20. L. A. Aguirre and C. Letellier, “Observability of multivariate differential embeddings,” Journal of Physics A, vol. 38, no. 28, pp. 6311–6326, 2005. View at: Publisher Site | Google Scholar | MathSciNet
  21. T. Stykel, “Gramian-based model reduction for descriptor systems,” Mathematics of Control, Signals, and Systems, vol. 16, no. 4, pp. 297–319, 2004. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  22. Y. Zhu and P. R. Pagilla, “Bounds on the solution of the time-varying linear matrix differential equation P˙(t)=AH(t)P(t)+P(t)A(t)+Q(t),” IMA Journal of Mathematical Control and Information, vol. 23, no. 3, pp. 269–277, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  23. K. Zhou, G. Salomon, and E. Wu, “Balanced realization and model reduction for unstable systems,” International Journal of Robust and Nonlinear Control, vol. 9, no. 3, pp. 183–198, 1999. View at: Google Scholar | Zentralblatt MATH | MathSciNet
  24. C. D. M. Martin and C. F. V. Loan, “Solving real linear systems with the complex Schur decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 1, pp. 177–183, 2007. View at: Publisher Site | Google Scholar
  25. R. Bartels and G. Stewart, “Solutions of the matrix equation ax+xb=c,” Communications of ACM, vol. 15, no. 9, pp. 820–826, 1972. View at: Google Scholar
  26. P. Pahlavanloo, “Dynamic aeroelastic simulation of the AGARD 445.6 wing using edge,” Tech. Rep. FOI-R-2259-SE, Defence Research Agency Defence and Security, System and Technology, Stockholm, Swedish, 2007. View at: Google Scholar

Copyright © 2013 Douglas Domingues Bueno et al. 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.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly as possible. We will be providing unlimited waivers of publication charges for accepted research articles as well as case reports and case series related to COVID-19. Review articles are excluded from this waiver policy. Sign up here as a reviewer to help fast-track new submissions.