The Use of Gramian Matrices for Aeroelastic Stability Analysis

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.Themain goal of this work is to introduce how observability gramianmatrix 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.


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 statespace models [6] and have been used for applications such as optimal placement of sensors and actuators [7][8][9].After this work the approach was extended for unsteady aerodynamic and aeroelastic systems [10][11][12][13].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.

Aeroelastic Model for Stability Analysis
Assuming a general aeroelastic model that can be written in the state-space form, according to where x() is the state vector, A is the dynamic matrix, B  is the input matrix, f  () is a vector of external forces, C is the known as the output matrix, and y() 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 A is stable for a range of velocities defined by increasing airspeed values  1 , . . .,   , . . .,   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.

Observability and Gramian Matrices
The concept of observability involves the dynamic matrix A and the output matrix C. A linear system, or the pair (A, C), is observable at instant of time  0 , if the state x( 0 ) can be determined from the output y(), with  ∈ [ 0 ,  1 ], where  1 >  0 is a finite instant of time.If this is true for all initial time  0 and all initial states x( 0 ), 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 2 2 (2 +  lag ) × 2 2 (2 +  lag ) has hank (2 +  lag ).With the observability matrix given by . . .
where  is the dimension of matrix A. 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 largeorder 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 x for which Ax = 0), 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 y at time  0 caused by the system's initial state x(0) 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, W  , introduced by Zhou et al. [23], defined for both stable and unstable systems.Then, let X  be the solution to the Riccati equation The observability Gramian matrix W  is a submatrix of W  which can be computed by solving the following equation where and the modified output matrix C is used instead of C.This modified output matrix C with dimension  × (2 +  lag ) is defined such that where the th row C  satisfies the equation   () =  () (), and the other ( − 1) rows are filled with zeros.Thus, considering C  and the aeroelastic matrix A(  ) defined at the airspeed   , the observability Gramian matrix W  (C  ,   ) = W  (, ) is computed by solving (9) for the pair (C  ,   ).The input matrix to represent  inputs.

Complex Schur Decomposition.
Using a complex Schur decomposition, the Lyapunov equation can be reshaped as a real linear system of equation and solving it, W () is obtained [24] [I ⊗ ( where I 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].

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   (, ) ∈ R + is defined and obtained by computing a matrix norm of W  .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 C  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  V is the length of the airspeed vector The th Gramian parameter in each column of the matrix Σ containing  max  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  max  can be plotted to determine the airspeed of flutter.(This is illustrated in Figure 1.)

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   ,   = 1, . . ., (2 +  lag ).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.

Numerical Simulations
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  = 12.7 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.

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  REF = 1.225 kg/m 3 .The values of reduced frequencies are presented in Appendix C. The model has a reference length of 2 = 0.5578 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,  1 = 9.45 Hz,  2 = 39.69Hz,  3 = 49.45Hz, and  4 = 95.10Hz.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  1 = 0.55,  2 = 1.40,  3 = 1.90, and  4 = 2.90.The observability Gramian matrix was computed for each pair (C  ,   ) and   was computed using different matrix norms.Analysis was evaluated considering 61 pairs of (  ,   ) with constant Mach number.Using the traditional pk-method, the flutter speed was identified at 158.02 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.

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 pkmethod.
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.
(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.

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 M, the viscous damping D, and the stiffness K.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 u  .The aerodynamic influence matrix Q depends on the parameters  (reduced frequency) and   (Mach Number) and  is the dynamic pressure caused by the flowing fluid (Q 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 Q 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 u() 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  lag is the number of lag terms and   is the th lag parameter ( = 1, . . .,  lag ).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 x() = { u  u  u  }  , where u  are states of lags required for the approximation.The dynamic matrix A is given by where the matrices M  , D  , and K  are comprised of terms containing structural and aerodynamic coefficients which are given by (A.4)

B. Norms for Computing the Gramian Parameter
The matrix norms for computing the Gramian parameters are presented as follows.

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].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).

Figure 1 :
Figure 1: Gramian parameter used to detect the flutter phenomenon.

Figure 5 :
Figure 5: Gramian parameter computed by the Frobenius norm as a function of the airspeed (typical section airfoil).

Figure 6 :
Figure 6: Aeroelastic frequency as a function of the air density (AGARD wing).

Figure 7 :
Figure 7: Aeroelastic damping as a function of the air density (AGARD wing).

6 MathematicalFigure 8 :Figure 9 :
Figure 8: Gramian parameters computed by the Frobenius norm as a function of the air density (AGARD wing).

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

Table 1 :
Physical and geometric properties of the 2D airfoil.

Table 2 :
Reduced frequencies for the case study AGARD wing 445.6.