# The Use of Gramian Matrices for Aeroelastic Stability Analysis

Authors: Douglas Domingues Bueno; Clayton Rodrigo Marqui; Luiz Carlos Sandoval Góes; Paulo José Paupitz Gonçalves
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 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 k-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 [7–9]. After this work the approach was extended for unsteady aerodynamic and aeroelastic systems [10–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.


## 2. Aeroelastic Model for Stability Analysis

Assuming a general aeroelastic model that can be written in the state-space form, according to(1)x˙(t)=Ax(t)+Bcfc(t),y(t)=Cx(t),
where x(t) is the state vector, A is the dynamic matrix, Bc is the input matrix, fc(t) is a vector of external forces, C is the known as the output matrix, and y(t) 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 V1,…,Vj,…,Vn containing a flutter airspeed VF, if V<VF.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 matrixA and the output matrix C. A linear system, or the pair (A,C), is observable at instant of time t0, if the state x(t0) can be determined from the output y(t), with t∈[t0,t1], where t1>t0 is a finite instant of time. If this is true for all initial time t0 and all initial states x(t0), the system is said to be completely observable [14].A linear time-invariant system withm outputs is said to be completely observable, if and only if the observability matrix with dimension 2m2(2+nlag)×2m2(2+nlag) has hank m(2+nlag). With the observability matrix given by
(2)O=[CCACA2⋮CAn-1],
where n 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 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 timet<∞ written as
(3)Wo(t)=∫0t‍[eATtCTCeAt]dt
which, according to [14], can be determinated as
(4)W˙o(t)=ATWo(t)+Wo(t)A+CTC,
where W˙o(t) 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]
(5)ATWo+WoA+CTC=0.An important property between observability and Gramian matrices is that they share the same Kernel (i.e., the set of all vectorsx for which Ax=0),
(6)Ker[Wo(t)]=Ker[O(C,A)].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 t0 caused by the system's initial state x(0) such that
(7)Energy[y(t0)]=xT(0)Wo(t0)x(0).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, Wco, introduced by Zhou et al. [23], defined for both stable and unstable systems. Then, let Xg be the solution to the Riccati equation
(8)XgFg+FgTXg-XgGgGgTXg=0.
The observability Gramian matrix Wo is a submatrix of WG which can be computed by solving the following equation
(9)(Fg+GgMg)WG+WG(Fg+GgMg)T+GgGgT=0,
where
(10)Fg=[A00AT],Gg=[BcC-T],Mg=-GgTXg,WG=[WcWcoWcoWo],
and the modified output matrix C- is used instead of C. This modified output matrix C- with dimension m×m(2+nlag) is defined such that
(11)y-(t)=C-x(t)y-(t)={0⋯0um(i)(t)0⋯0}Tthen,C-i=[0⋯0(m+(i-1))1(i)0⋯0],
where the ith row C-i satisfies the equation y-i(t)=um(i)(t), and the other (m-1) rows are filled with zeros. Thus, considering C-i and the aeroelastic matrix A(Vj) defined at the airspeed Vj, the observability Gramian matrix Wo(C-i,Vj)=Wo(i,j) is computed by solving (9) for the pair (C-i,Vj). The input matrix Bc is written by [Mam-10m(1+nlag)×m]T to represent m 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,Wo(ij) is obtained [24]
(12)[I⊗(Fg+GgMg)+(FgT+MgTGgT)⊗I]vec(WG)=vec(-GgGgT),
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].


### 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σg(i,j)∈R+ is defined and obtained by computing a matrix norm of Wo. 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σg(i,j) for each pair Vj and C-i it is possible to build a matrix Σ (14). Assuming that aeroelastic instability occurs in this speed range, the flutter speed Vj is found for the largest value of σg(i,j):
(13)max(Σ)=maxij[σg]=σgmaxsuchthat,ifΔV⟶dV⟹dσgdV=0,
where nv is the length of the airspeed vector
(14)Σ=[σg(1,1)⋯σg(1,j)⋯σg(1,nv)⋯⋯⋯⋯⋯σg(i,1)⋯σg(i,j)⋯σg(i,nv)⋯⋯⋯⋯⋯σg(m,1)⋯σg(m,j)⋯σg(m,nv)].Theith Gramian parameter in each column of the matrix Σ containing σgmax is related to a measure of the energy absorbed by the ith aeroelastic mode. The values of σg in each column can then be compared to determine the mode contribution. On the other hand, the row of Σ which contains σgmax can be plotted to determine the airspeed of flutter. (This is illustrated in Figure 1.)Figure 1
Gramian parameter used to detect the flutter phenomenon.


### 3.3. Matrix Norm

In practical implementations the Gramian parameterσg(i,j) 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,
(15)σg(F)(i,j)=(∑ir,jc‍|wo(ir,jc)|)1/2,
where wo(ir,jc) is the Gramian matrix element of the irth row and jcth column, and ir,jc=1,…,m(2+nlag).


## 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,Wo(ij) is obtained [24]
(12)[I⊗(Fg+GgMg)+(FgT+MgTGgT)⊗I]vec(WG)=vec(-GgGgT),
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].


## 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σg(i,j)∈R+ is defined and obtained by computing a matrix norm of Wo. 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σg(i,j) for each pair Vj and C-i it is possible to build a matrix Σ (14). Assuming that aeroelastic instability occurs in this speed range, the flutter speed Vj is found for the largest value of σg(i,j):
(13)max(Σ)=maxij[σg]=σgmaxsuchthat,ifΔV⟶dV⟹dσgdV=0,
where nv is the length of the airspeed vector
(14)Σ=[σg(1,1)⋯σg(1,j)⋯σg(1,nv)⋯⋯⋯⋯⋯σg(i,1)⋯σg(i,j)⋯σg(i,nv)⋯⋯⋯⋯⋯σg(m,1)⋯σg(m,j)⋯σg(m,nv)].Theith Gramian parameter in each column of the matrix Σ containing σgmax is related to a measure of the energy absorbed by the ith aeroelastic mode. The values of σg in each column can then be compared to determine the mode contribution. On the other hand, the row of Σ which contains σgmax can be plotted to determine the airspeed of flutter. (This is illustrated in Figure 1.)Figure 1
Gramian parameter used to detect the flutter phenomenon.


## 3.3. Matrix Norm

In practical implementations the Gramian parameterσg(i,j) 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,
(15)σg(F)(i,j)=(∑ir,jc‍|wo(ir,jc)|)1/2,
where wo(ir,jc) is the Gramian matrix element of the irth row and jcth column, and ir,jc=1,…,m(2+nlag).


## 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 (semichordb). 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.Table 1
Physical and geometric properties of the 2D airfoil.
Aerodynamic semichord
b
=
0.15 m
Airfoil mass
M
=
5.0 kg
Plunge frequency
ω
h
=
3.0 Hz
Pitch frequency
ω
θ
=
4.5 Hz
Control surface rotation frequency
ω
β
=
12 Hz
Air density
ρ
=
1.2895 kg/m3
Lag parameters (Roger's method)
β
1
=
0.2
β
2
=
1.2
β
3
=
1.6
β
4
=
1.8
Reduced frequencies
0.1
≤
k
≤
2.0
;
Δ
k
=
0.1
Figure2
a
=
-
0.40
Figure2
c
=
0.60
Distance between c.e. to c.g.
x
-
θ
=
0.20
Distance between c.e. to c.g. (flap)
x
-
β
=
0.0125
Radius of gyration of the flap referred toa
r
-
β
=
(
6.25
×
1
0
-
3
)
1
/
2
Radius of gyration of the airfoil referred toa
r
-
θ
=
0.25
Elastic center
c.e.
Center of gravity
c.g.Figure 2
Typical section 2D airfoil.The results of the proposed approach were compared with classical eigenvalues analysis. Figures3 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 V=12.7 m/s.Figure 3
Aeroelastic frequency as a function of airspeed (typical section airfoil).Figure 4
Real part of the eigenvalues as a function of the airspeed (typical section airfoil).The results of the proposed method are presented in Figure5. The flutter speed was correctly identified. Similar results were also obtained using the other norms and therefore they are not shown herein.Gramian parameter computed by the Frobenius norm as a function of the airspeed (typical section airfoil).
(a)
(b)
(c)


### 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 AppendixC 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/m3. The values of reduced frequencies are presented in Appendix C. The model has a reference length of 2b=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, f1=9.45 Hz, f2=39.69 Hz, f3=49.45 Hz, and f4=95.10 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β1=0.55, β2=1.40, β3=1.90, and β4=2.90. The observability Gramian matrix was computed for each pair (C-i,Vj) and σg was computed using different matrix norms. Analysis was evaluated considering 61 pairs of (Vj,ρj) with constant Mach number. Using the traditional pk-method, the flutter speed was identified at 158.02 m/s. The classical V-f and V-g diagrams are shown in Figures 6 and 7.Figure 6
Aeroelastic frequency as a function of the air density (AGARD wing).Figure 7
Aeroelastic damping as a function of the air density (AGARD wing).With the proposed methodology, the Gramian parameterσg 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.Gramian parameters computed by the Frobenius norm as a function of the air density (AGARD wing).
(a)
(b)Gramian parameters computed by the Frobenius norm as a function of the air density (AGARD wing).
(a)
(b)


### 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.Figure10 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 P(ρj,Vj) into the flight envelope.(2)
Identify the pointP(ρj,Vj)=P(j) at which Gramian parameter is maximum, according to (13).(3)
Compute cross-Gramian parameters for points around the maximum.Figure 10
Comparison of computational time.Using this procedure, it is possible to obtain some reduction in the computational time to evaluate the system stability. This reduction is represented in Figure11.Figure 11
Computational time: pk-method versus Gramian parameters (AGARD 445.6 wing).


## 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 (semichordb). 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.Table 1
Physical and geometric properties of the 2D airfoil.
Aerodynamic semichord
b
=
0.15 m
Airfoil mass
M
=
5.0 kg
Plunge frequency
ω
h
=
3.0 Hz
Pitch frequency
ω
θ
=
4.5 Hz
Control surface rotation frequency
ω
β
=
12 Hz
Air density
ρ
=
1.2895 kg/m3
Lag parameters (Roger's method)
β
1
=
0.2
β
2
=
1.2
β
3
=
1.6
β
4
=
1.8
Reduced frequencies
0.1
≤
k
≤
2.0
;
Δ
k
=
0.1
Figure2
a
=
-
0.40
Figure2
c
=
0.60
Distance between c.e. to c.g.
x
-
θ
=
0.20
Distance between c.e. to c.g. (flap)
x
-
β
=
0.0125
Radius of gyration of the flap referred toa
r
-
β
=
(
6.25
×
1
0
-
3
)
1
/
2
Radius of gyration of the airfoil referred toa
r
-
θ
=
0.25
Elastic center
c.e.
Center of gravity
c.g.Figure 2
Typical section 2D airfoil.The results of the proposed approach were compared with classical eigenvalues analysis. Figures3 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 V=12.7 m/s.Figure 3
Aeroelastic frequency as a function of airspeed (typical section airfoil).Figure 4
Real part of the eigenvalues as a function of the airspeed (typical section airfoil).The results of the proposed method are presented in Figure5. The flutter speed was correctly identified. Similar results were also obtained using the other norms and therefore they are not shown herein.Gramian parameter computed by the Frobenius norm as a function of the airspeed (typical section airfoil).
(a)
(b)
(c)


## 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 AppendixC 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/m3. The values of reduced frequencies are presented in Appendix C. The model has a reference length of 2b=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, f1=9.45 Hz, f2=39.69 Hz, f3=49.45 Hz, and f4=95.10 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β1=0.55, β2=1.40, β3=1.90, and β4=2.90. The observability Gramian matrix was computed for each pair (C-i,Vj) and σg was computed using different matrix norms. Analysis was evaluated considering 61 pairs of (Vj,ρj) with constant Mach number. Using the traditional pk-method, the flutter speed was identified at 158.02 m/s. The classical V-f and V-g diagrams are shown in Figures 6 and 7.Figure 6
Aeroelastic frequency as a function of the air density (AGARD wing).Figure 7
Aeroelastic damping as a function of the air density (AGARD wing).With the proposed methodology, the Gramian parameterσg 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.Gramian parameters computed by the Frobenius norm as a function of the air density (AGARD wing).
(a)
(b)Gramian parameters computed by the Frobenius norm as a function of the air density (AGARD wing).
(a)
(b)


## 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.Figure10 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 P(ρj,Vj) into the flight envelope.(2)
Identify the pointP(ρj,Vj)=P(j) at which Gramian parameter is maximum, according to (13).(3)
Compute cross-Gramian parameters for points around the maximum.Figure 10
Comparison of computational time.Using this procedure, it is possible to obtain some reduction in the computational time to evaluate the system stability. This reduction is represented in Figure11.Figure 11
Computational time: pk-method versus Gramian parameters (AGARD 445.6 wing).


## 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.



---

*Source: 941689-2013-04-10.xml*