Aeroservoelastic Modeling of Body Freedom Flutter for Control System Design

Aeroservoelastic Modeling of Body Freedom Flutter for

Control System Design

∗ Jeffrey Ouellette NASA Armstrong Flight Research Center, Edwards, California, 93523 One of the most severe forms of coupling between aeroelasticity and flight dynamics is an instability called body freedom flutter. The existing tools often assume relatively weak coupling, and are therefore unable to accurately model body freedom flutter. Because the existing tools were developed from traditional flutter analysis models, inconsistencies in the final models are not compatible with control system design tools. To resolve these issues, a number of small, but significant changes have been made to the existing approaches. A frequency domain transformation is used with the unsteady aerodynamics to ensure a more physically consistent stability axis rational function approximation of the unsteady aerodynamic model. The aerodynamic model is augmented with additional terms to account for limitations of the baseline unsteady aerodynamic model and to account for the gravity forces. An assumed modes method is used for the structural model to ensure a consistent definition of the aircraft states across the flight envelope. The X-56A stiff wing flight-test data were used to validate the current modeling approach. The flight-test data does not show body-freedom flutter, but does show coupling between the flight dynamics and the aeroelastic dynamics and the effects of the fuel weight.

Nomenclature A Steady aerodynamic coefficients A Quasi-steady aerodynamic coefficients A Added mass aerodynamic coefficients AIC Aerodynamic influence coefficients ¯ B Effective stability axis damping matrix B Structrual damping matrix b Wing span ¯ c Mean aerodynamic chord C Aerodynamic drag force coefficient at trim D C Aerodynamic lift force coefficient at trim L C Aerodynamic roll moment coefficient at trim l C Aerodynamic pitch moment coefficient at trim m C Aerodynamic yaw moment coefficient at trim n C Aerodynamic side force coefficient at trim Y C Aerodynamic modal force coefficient at trim ηh D Aerodynamic lag output matrix D Velocity nondimensionalization matrix u D Position nondimensionalization matrix x DoF Degrees of freedom E Aerodynamic lag input matrix FEM Finite element model G Vector of forces on each grid point due to gravity ∗ Aerospace Engineer, Controls and Dynamics Branch, P.O. Box 273/Mailstop 4840D, Edwards, California, 93523, AIAA Member.

1 of 19 American Institute of Aeronautics and Astronautics g Acceleration due to gravity GAF Generalized aerodynamic forces I Identity matrix Im Imaginart part i Imaginary number ¯ K Effective stability axis stiffness matrix K Stiffness matrix K Finite element stiffness matrix gg K Kernel (null space) ˜ K Assumed mode stiffness matrix hh k Reduced frequency ¯ M Effective stability axis mass matrix M Mass matrix M Finite element mass matrix gg ˜ M Assumed mode mass matrix hh MUTT Multi-utility technology testbed p Roll rate ¯ q Dynamic pressure Q Generalized aerodynamic force matrix q Generalized modal forces ˜ Q Assumed mode generalized aerodynamic force matrix hh q Pitch rate R Aerodynamic lag pole matrix Re Real part r Yaw rate RFA Rational function approximation ˆ T ∗ Transformation earth fixed axis to modal rate e h ˆ T Transformation earth fixed axis to modal displacement eh ˆ T ∗ Transformation stability axis to modal rate s h T Transformation matrix of gravity components in the stability axis grav ˜ T ( k ) Frequency domain transformation from modal coordinates to stability coordinates ˜ T Steady state transformation from modal coordinates to stability coordinates T Kinetic energy t Time u Velocity vector V Potential energy V Reference airspeed (true airspeed at trim) V Equivalent airspeed e x Displacement vector x Earth-fixed longitudinal position y Earth-fixed lateral position z Earth-fixed vertical position Subscripts aug Augmentation c Control Surface Coordinate System e Earth fixed coordinate system g Finite element coordinate system grav Contribution of gravity h Modal coordinate system s Stability (flight dynamics) coordinate system Superscripts ∗ Stability axis rational function approximation 2 of 19 American Institute of Aeronautics and Astronautics Symbols α Angle of attack β Angle of sideslip δ Control surface inputs η Modal coordinate θ Pitch angle φ Roll angle Φ Assumed mode shape ψ Yaw angle Conventions ˙ () Derivative with respect to time ˆ () Non-dimensionalized [ ] Matrix r [ · ] Diagonal matrix r { } Vector ∗ () Derivative with respect to nondimensional time I. Introduction s an aircraft airspeed is increased, the frequency of the short period mode increases which strengthens the A interactions between the rigid body short-period and the wing structural dynamics. When a vehicle has a lightly damped short-period, such as the General Dynamics (Falls Church, Virginia) RB-57H Canberra, and low frequency structural modes, the interaction can produce an instability called body freedom flutter.

The use of an active flight control system offers an appealing solution for suppressing the body freedom flutter instability. Developing a flutter suppression control system for body freedom flutter requires a model that considers flight dynamics, the structural dynamics, and the interaction. Both flight dynamics and aeroelasticity engineers have many sophisticated and standardized methods of generating models that are optimized to satisfy the needs of a given disipline. As modern aircraft become more flexible and these disciplines converge, inconsistencies between the independently developed modeling methodologies arise. The level of coupling present in body freedom flutter requires an integrated approach to model development that will avoid these conflicts.

In modern aeroelasticity the vehicle structure is discretized into many finite elements. Because the primary concern is the deformations of the structure, these finite element models (FEM)s use an inertial coordinate system. The FEM is often impractical for direct simulation due to the many thousands of degrees of freedom and small errors (e.g. local modes) due to discretization of the vehicle structure. For aeroelastic analysis, the order of the structural dynamics is reduced by considering only a subset of the vehicle structural modes.

The introduction of aerodynamics introduces coupling between these in vacuo modes. As a result, accurately capturing the low frequency dynamics still requires higher frequency modes. At these higher frequencies, the rate of response of the aerodynamics is on the same order as the structural dynamics. To correctly capture the changing of the airflow over the vehicle as the structure deforms requires an unsteady aerodynamic model. These high fidelity models include the airflow changing over the aircraft and traveling into the wake.

Generally, the aeroelastic models are designed considering a fixed airspeed. The drag forces act primarily on the velocity degree of freedom. Therefore, by removing the velocity degree of freedom a model of the drag is not necessary, but removing the velocity makes it impossible to model the phugoid. Often many of the tools used for aeroelasticity have been based on a frequency domain approach. However, the Helios (Aerovironment, Monrovia, California) mishap investigation showed the importance of a time domain form of the models for evaluating an aircraft closed loop performance.

In contrast to aeroelastic models, the rigid body flight dynamics models are of a much lower order, typically consisting of 6 degrees of freedom (DoF). Because these models do include the velocity dynamics, the correct modeling of the flight dynamics requires the consideration of non-inertial kinematics. As a result, the vehicle velocity is not equal to the derivative of the position. Unlike the structural models, the inclusion of the vehicle velocity in flight dynamics means a greater influence from drag. The drag of the vehicle can be a significant dissipation of energy. The flight dynamics also require the consideration of the gravitational forces.

3 of 19 American Institute of Aeronautics and Astronautics Correctly modeling the phugoid mode (a low frequency, lightly damped, and large magnitude oscillation in the vehicle altitude) requires both the drag and gravity effects.

To resolve the conflicts between the aeroelastic and flight dynamics modeling, the current work has returned to the fundamental principles shared by both these disciplines to develop a new modeling approach that is consistent with both applications. An effort has been made to keep these models similar to the existing models, so that superficially they appear very similar to the current state of the art. However, addressing the conflicts between the models from the independent disciplines has enabled a number of incremental improvements in the models. These include an aerodynamic model that uses a more consistent definition of the inputs and a very general formulation of gravitational forces.

The X-56A (Lockheed Martin, Bethesda Maryland) multi-utility technology testbed (MUTT) was developed specifically to study active control the sturcture of an aircraft and suppression of body freedom flutter. The wings of the aircraft can be exchanged to fly an aircraft with a relatively stiff structure or a flexible aircraft with unstable flutter modes. Unlike production aircraft today, the X-56A MUTT has structural dynamics that are acutely coupled with the flight dynamics. The stiff wings on the X-56A MUTT do not have instabilities in the flight envelope, but do still demonstrate significant interactions between the flight dynamics and structure. At low airspeeds, both wings have similar behavior. Therefore, the stiff wing flight data are used as a validation of the methods before proceeding to the flexible wing flights that do have the instabilities.

Developing models of an aircraft with body freedom flutter like the X-56A MUTT requires returning to the origins of both these modeling approaches. The resulting models are capable of capturing the requirements of both the flight controls engineers and the aeroelasticity engineer. To describe the modeling approach, first the coordinate systems used for each discipline are defined and the transformations that relate these coordinate systems. The aerodynamic modeling methodology is described that is based heavily on existing aeroelasticity modeling tools to generate a time domain model that is compatible with state space methods for control laws and is more consistent with aerodynamic models seen in flight dynamics. Augmentations to this aerodynamic model are described to address the limitations of the aeroelastic tools such as a fixed velocity and lack of gravitational forces. Finally, the approach was applied to the X-56A stiff wing vehicle and is compared against the available flight-test data.

II. Modeling Approach Most flight control development requires the models to have a smaller number of states. To achieve these requirements, a state space representation of the models was selected. These models are required to include unsteady aerodynamic effects to capture the phase difference between the vehicle dynamics and the aerodynamic forces. The importance of the unsteady aerodynamics increases with the influence of the structural dynamics.

Consideration of the low frequency dynamics is very important if the vehicle has low or negative static margin. The low static margin can cause a coupling between the statically unstable short-period and the phugoid mode. The coupling will result in an unstable phugoid-like mode, that includes more pitch response than in a traditional phugoid mode. The damping of these low frequency dynamics is dominated by the variations in aerodynamic forces with velocity. Most aeroelastic methods assume a fixed velocity for the analysis. The frequency of the low frequency phugoid mode is dominated by the gravitational forces which are also neglected in many aeroelastic methods.

A. Definition of Coordinate Systems and States The vehicle coordinate systems are separated into global and local coordinate systems. The global coordinate systems shown in figure 1, describe the total motion of the vehicle’s flight dynamics and structural dynamics as a complete system.

4 of 19 American Institute of Aeronautics and Astronautics M o d a l a x i s Z h Wind X Y s s E a r t h X h a x i s Z s S t a b i l i t y X Y /Y e h e a x i s Z e Figure 1. Global coordinate systems.

The primary global coordinate system is the stability axis system. The stability axis system is used to describe the orientation of the vehicle and the structural deformations with respect to the relative wind vector. The rigid body portion of the stability axis system is identical to traditional flight dynamics. The x-axis of the stability axis system is oriented to align with the vehicle relative wind vector. The first 12 stability axis states are defined by separate vectors of positions, Eq. (1), and rates, Eq. (2).

{ } { x } = x y z φ θ ψ (1) { } { u } = (2) V β α p q r e The angle of attack, α , and angle of sideslip, β , are angles that describe the vehicle’s velocity vector. The 3, 7 rigid body dynamics and structural dynamics are related by a mean axis system. For a linear system, orthogonality of the rigid body modes and the structural modes is sufficient to define the mean axis. The mean axis system simplifies the equations of motion by decoupling the kinetic and potential energy of the rigid body and structural dynamics. The coupling is instead captured by the aerodynamic forces. The mean axis system is more difficult to visualize because it is not fixed to a point on the aircraft. Rather it must move as the vehicle deforms. Specifically, the structural deformations are described as a linear combination of a finite set of mode shapes.

For a flat non-rotating earth, the earth-fixed coordinate system provides an inertial reference frame. The coordinates of the earth-fixed coordinate system describe the position and the orientation of the stability axis system. It is necessary to consider the orientation of the earth fixed coordinate system to capture the gravitational forces acting on the aircraft. The sign convention for the earth fixed coordinate system matches the convention used in standard rigid body flight dynamics simulations.

The final global coordinate system, the modal coordinate system, is also an inertial reference frame.

The modal coordinate system fully describes the vehicle dynamics and is traditionally used in aeroelasticity modeling. The pedigree of these existing tools is desirable as a starting place for the flight dynamics models, but is inadequate for describing velocity variations. Because the modal coordinate system is an inertial frame, a level trim condition is not an equilibrium in the modal coordinate system. If the modal coordinate system is used to define the linear equations of motion, then the linear model cannot be used to make any conclusions about the vehicle speed stability. The orientation of the modal coordinate system differs from the earth-fixed coordinate system to match the conventions used in aeroelastic modeling software.

The local coordinate systems are used to describe the dynamics of discrete elements of the vehicle. The local coordinates are described as sets of coordinate systems. The g -set of local degrees of freedom describe the grid points of the FEM. There may be many different coordinate systems used in the FEM. The definition of these coordinate systems is arbitrary and at the discretion of the engineer developing the FEM. The determination of these coordinate systems is beyond the scope of the current work, but knowledge of their definition from the FEM is required. The g -set coordinates are related to the modal coordinates by the mode 5 of 19 American Institute of Aeronautics and Astronautics shapes. Most finite element software can provide these mode shapes. The other local coordinates are the control surface, c -set, deflections.

The standard orthonormal mode shapes from a FEM have the benefit of being orthogonal with respect to both the FEM mass and stiffness matrices. Because the mode shapes change with the vehicle weight, the definition of the modal states would have to change with fuel conditions. However, the synthesis of a full envelope control design requires a consistent definition of the states. To ensure consistency, a single set of [ ] ˜ assumed mode shapes, [Φ ] , are used for all fuel cases. For these assumed modes, the mass, M , and 0 hh [ ] ˜ stiffness, K , matrices are specified to capture the kinetic and potential energy of the structure given the hh assumed mode shapes. This approach is often referred to as the assumed modes method. Using the full FEM mass, [ M ] , and stiffness matrices, [ K ] , kinetic energy in Eq. (3), T , and potential energy in Eq. (4), gg gg V , are a function of the mode shape, modal coordinate, { η } , and the FEM mass and stiffness matrices.

T T T = { ˙ η } [Φ ] [ M ] [Φ ] { ˙ η } (3) 0 gg 0 T T V = { η } [Φ ] [ K ] [Φ ] { η } (4) 0 gg 0 Similarly, in Eqs. (5) and (6) the kinetic and potential energy are expressed as a function of the modal mass and stiffness matrices.

[ ] T ˜ T = { ˙ η } M { ˙ η } (5) hh [ ] T ˜ V = { η } K { η } (6) hh Comparing the two forms of the energy, the assumed mass and stiffness matrix are a function of the mode shape and the FEM mass and stiffness matrices in Eqs. (7) and (8).

[ ] T ˜ M = [Φ ] [ M ] [Φ ] (7) hh 0 gg 0 [ ] T ˜ K = [Φ ] [ K ] [Φ ] (8) hh 0 gg 0 The mass and stiffness matrices in Eqs. (7) and (8) are no longer diagonal. However, for a sufficiently representative set of modes, the assumed mode method is able to very effectively model the structural dynamics. The determination of what is representative does reflect an increase in the complexity of generating the structural model. The application of the assumed modes method requires a sufficiently large set of assumed modes, [Φ ] , that is able to represent the deformed shape at all fuel cases. At a minimum, the assumed modes method will typically require more mode shapes than would be needed for the exact FEM modes. For more complex aircraft configurations, there can be local modes that will only appear at specific mass cases. These modes must be included in the set of assumed modes. At cases where the local modes would not normally appear, the local mode will have a very high frequency outside the bandwidth of the model.

[ ] ˜ The generalized aerodynamic forces (GAF), Q , in Eq. (9) are calculated from the aerodynamic hh influence coefficients (AIC).

[ ] T ˜ Q = [Φ ] [ AIC ( k )] [Φ ] (9) hh 0 0 Like the modal stiffness matrices, the AIC does not change with the fuel weight. Therefore the coefficients are consistent for all mass cases.

Because the unsteady aerodynamics are a function of time, generation of non-dimensional coefficients require a non-dimensional definition of time. The nondimensionalization of time, Eq. (10), is determined by the true airspeed at trim and the chord length.

2 V ˆ t = t (10) ¯ c 6 of 19 American Institute of Aeronautics and Astronautics To allow a consistent aerodynamic model across the flight envelope, it is desirable to define a set of nondimensional states for use in the aerodynamic model. The individual states are grouped together into non-dimensional position, { ˆ x } in Eq. (11), and rates, { ˆ u } in Eq. (12).

{ } { ˆ x } = (11) ˆ x ˆ y ˆ z φ θ ψ η { } ∗ ˆ { ˆ u } = (12) V β α ˆ p ˆ q ˆ r η ˆ ˆ The non-dimensional position, { x } in Eq. (13), and rates, { u } in Eq. (14), are related to the dimensional r r states by the diagonal matrices, [ D ] in Eq. (15) and [ D ] in Eq. (16).

x r u r [ ] r { ˆ x } = D { x } (13) x r [ ] r { ˆ u } = D { u } (14) u r These matrices are partitioned to distinguish the rigid and flexible transformations.

  0 0 0 0 0 0 ¯ c    0 0 0 0 0 0  ¯ c     0 0 0 0 0 0 ¯ c   [ ] r   D = 0 0 0 1 0 0 0 (15) x r     0 0 0 0 1 0 0      0 0 0 0 0 1 0  r 0 0 0 0 0 0 I r   0 0 0 0 0 0 V     0 1 0 0 0 0 0     0 0 1 0 0 0 0   [ ]   ¯ c r 0 0 0 0 0 0 D =   (16) u r 2 V   ¯ c   0 0 0 0 0 0  2 V    ¯ c 0 0 0 0 0 0   2 V ¯ c r 0 0 0 0 0 0 I r 2 V B. Transformations The structural and unsteady aerodynamic model are defined in the modal coordinate system. However, because the modal coordinate system is an inertial coordinate system, it is insufficient for describing the vehicle flight dynamics. It is also easier for a controls engineer to understand a non-inertial stability axis system, because that is a more traditional form used in flight dynamics modeling and simulation. The transformations for the structural and unsteady aerodynamic models are first derived in the time domain.

The time domain form is used directly for the structural models and to generate the final state space models.

For the aerodynamic model generation, a separate frequency domain form of the transformation is derived from the time domain form.

1. Time Domain Transformation Using the time domain allows for the kinematic relationship between the modal and stability axis to be expressed in an exact nonlinear form. For an aerodynamic model derived in the modal coordinate system, the first six modes are the dimensional rigid body translations and rotations. The displacements are nondimensionalized by the chord length. The x-axis and z-axis of the structural and stability axis coordinate systems are reversed. The modal coordinates in Eq. (17) are expressed as a function of these stability axis displacements.

{ } T ¯ c ¯ c ¯ c { η } = (17) − ˆ x ˆ y − ˆ z − φ θ − ψ 2 2 2 7 of 19 American Institute of Aeronautics and Astronautics ˆ The derivatives of these states with respect to nondimensional time, t , are given by the nondimensionalization of the Euler angle and navigation equations traditionally used in nonlinear flight dynamics in Eq. (18).

  ¯ c ¯ c ˆ ˆ   − V cos α cos β cos θ cos ψ − V sin β ( − cos φ sin ψ + sin φ sin θ cos ψ )   2 2      ¯ c  ˆ   − V sin α cos β (sin φ sin ψ + cos φ sin θ cos ψ )       ¯ c ¯ c   ˆ ˆ  V cos α cos β cos θ sin ψ + V sin β (cos φ cos ψ + sin φ sin θ sin ψ )   2 2      { } ¯ c ˆ ∗ + V sin α cos β ( − sin φ cos ψ + cos φ sin θ sin ψ ) η = (18) ¯ c ¯ c ¯ c ˆ ˆ ˆ   V cos α cos β sin θ − V sin β sin φ cos θ − V sin α cos β cos φ cos θ   2 2 2         − ˆ p − tan θ (ˆ q sin φ + ˆ r cos φ )          ˆ q cos φ − ˆ r sin φ        − ˆ q sin φ sec θ − ˆ r cos φ sec θ These equations are linearized to give a linear transformation between the modal coordinates and the stability axis states in Eq. (19). In the time domain, the linear transformation between the modal coordinates and the mean stability axis states is non-singular.

{ } [ ] { } ∗ ˆ ˆ T ∗ T ∗ η ˆ u s h e h = (19) ˆ η ˆ x 0 T eh Where the matricies are defined by Eqs. (20), (21), and (22).

  ¯ c − 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 − 0 0 0 0 [ ]   ˆ   T = (20) 0 0 0 − 1 0 0 0 eh     0 0 0 0 1 0 0       0 0 0 0 0 − 1 0 r 0 0 0 0 0 0 I r   ¯ c − 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 − 0 0 0 0 [ ]   ˆ   T ∗ = (21) 0 0 0 − 1 0 0 0   s h   0 0 0 0 1 0 0       0 0 0 0 0 − 1 0 r 0 0 0 0 0 0 I r   0 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 0 0 0 0 [ ]   ˆ   T ∗ = (22) 0 0 0 0 0 0 0   e h   0 0 0 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 [ ] ˆ Because the stability axis is a non-inertial coordinate system, the matrix, T ∗ , is nonzero. However, it is e h still singular. The time domain allows for a clear derivation of the transformation. However, the aerodynamic model from many unsteady aerodynamic methods are defined in the frequency domain. Previous efforts have transformed the frequency domain model to a time domain model before applying the change of coordinates.

Fitting the rational function approximation (RFA) to the modal GAF reduces the total error, but uses the orientation states (heading and pitch angle). Because there are no conditions placed on the coefficients, the final RFA has coefficients that are a function of these orientation states. To avoid the introduction of these 8 of 19 American Institute of Aeronautics and Astronautics erroneous coefficients, the GAF is transformed into the frequency domain before the RFA is calculated. The frequency domain transformation isolates these non-aerodynamic states, and allows constraints to be placed on the resulting coefficients.

2. Frequency Domain Transformation To determine the frequency domain transformation matrices the transformation to the non-dimensional frequency domain in Eq. (23) is applied to the time domain transformation in Eq. (19). The non-dimensional reduced frequency response of the linear equation can be analytically determined in Eqs. (24) and (25).

{ } ∗ η = ik { η } (23) [ ] [ ] ˆ ˆ ik { η } = T ∗ { ˆ u } + T ∗ { ˆ x } (24) s h e h [ ] ˆ { η } = T { ˆ x } (25) eh These equations can be solved to get the modal states, { η ( k ) } , as a function of the stability axis velocities, { ˆ u ( k ) } , in Eqs. (26) and (27).

[ ] ˜ ˆ { η ( k ) } = T ( k ) { u ( k ) } (26) ( ) [ ] [ ] [ ] − 1 [ ] [ ] − 1 r ˜ ˆ ˆ ˆ T ( k ) = ik ik I − T ∗ T T ∗ ∀ k 6 = 0 (27) r eh e h s h [ ] ˆ Because the matrix T ∗ is singular, the transformation in Eq. (27) is undefined at zero reduced frequency.

e h Therefore, the case of zero reduced frequency must be handled separately.

The zero reduced frequency transformation is derivied from the frequency domain kinematic equations in Eq. (24) at zero reduced frequency, shown in Eq. (28).

[ ] [ ] [ ] − 1 ˆ ˆ ˆ { 0 } = T ∗ { ˆ u } + T ∗ T { η } (28) eh s h e h [ ] ˆ A a nontrivial solution to Eq. (28) is defined by the nullspace, kernel, of the singular matrix, T ∗ , in e h Eqs. (29) and (30).

[ ] ˜ { η ( k ) } = T { ˆ u ( k ) } (29) [ ] [ ( ) ] [ ] ˜ ˆ ˆ T = 0 K T ∗ T ∗ (30) e h s h The definition of the kernel is not unique, but is selected in Eq. (31) to ensure consistency in the transformation.

  1 0 0 0 0 0 0   0 1 0 0 0 0 0       0 0 1 0 0 0 0 [ ]   ˜   T = (31) 0 0 0 0 1 0 0 0      0 0 0 0 0 0 0      0 0 0 0 0 0 0 r 0 0 0 0 0 0 I r The transformation indicates that the fifth and sixth modal displacements are directly correlated to the angle of attack and side slip. The remaining modal displacements are independent of the stability axis velocities.

C. Aerodynamics Traditionally, flight dynamics models have utilized quasi-steady aerodynamics. As the frequency of the dynamics increases, as is seen in aeroelastic models, it becomes necessary to include a fully unsteady 9 of 19 American Institute of Aeronautics and Astronautics aerodynamic model. Exclusion of the unsteady aerodynamics results in an error in the phase between the modes. Because the relatively low frequency flight dynamics are strongly coupled with the structural dynamics, none of the unsteady aerodynamic coupling terms can be neglected.

Frequency domain aerodynamic models have traditionally offered a very computationally efficient method for the calculation of unsteady aerodynamic models. The constitutive equations for the unsteady aerodynamics are significantly simplified by expressing the aerodynamics in the frequency domain. Many traditional aeroelasticity tools, such as doublet lattice, use the frequency domain approach. These methods give a linear relationship in Eq. (32) called the GAF, [ Q ( k )] , that relates the modal states, η ( k ) , to the generalized modal forces, { q ( k ) } , at a specified reduced frequency, k .

{ q ( k ) } = ¯ q [ Q ( k )] { η ( k ) } (32) A separate GAF is calculated at each reduced frequency. For these results to be integrated into state space models needed for the control development, the GAF needs to be converted to the time domain.

The issue with frequency domain methods is that, even for a simple 2-D airfoil, the frequency response is an irrational non-linear function of the reduced frequency. A closed form time domain representation is impossible for an irrational frequency response function. The lack of a closed form solution, requires an RFA (a transfer function) to describe these unsteady aerodynamics. For flutter suppression experiments on a B-52 11, 12 aircraft (The Boeing Company, Chicago, Illinois), Roger developed the idea of an RFA that could be used for controller design. The general form of the RFA typically used is shown in Eq. (33).

([ ] [ ] [ ] [ ] ) ( [ ] [ ]) − 1 2 r r ˆ ˆ ˆ ˆ { ˜ q ( k ) } ≈ − A + A ik + A k + [ D ] ik I − R E ik { η } (33) 0 1 2 r r r The traditional Roger’s RFA prescribes a specific form for the [ R ] and [ D ] matrices. The aerodynamic lag r r poles in [ R ] must be defined by the engineer developing the model. The traditional Roger’s RFA specified r a number of poles that are repeated for each degree of freedom in the structural model. The matrix [ D ] then has a specified structure in Eq. (34).

  1 1 · · · 1 0 0 0 0 0 · · · 0    0 0 · · · 0 1 1 · · · 1 0 0 · · · 0    [ D ] ≡ (34) . . . . . . . .

 .  . . . . . . . . .

 .  . . . . . . . .

0 0 · · · 0 0 0 · · · 0 1 1 · · · 1 Other structures for the RFA matrices exist, such as Karpel’s minimum state. These can be helpful in reducing the order of the rational function approximation. The aerodynamic lags are not associated with a specific degree of freedom. As a result, it is more difficult to ensure a consistent definition of these aerodynamic states. The Roger’s RFA remains popular because the synthesis of the RFA is much more predictable.

The issue with all of the prior traditional RFA forms in Eq. (33) is that they all use the modal coordinate system. The traditional RFA worked well when the structural dynamics could be decoupled from the flight dynamics and vice versa. To be able to model the flight dynamics, as required for body freedom flutter, it is necessary for the equations of motion to be expressed in a non-inertial reference frame. Instead of taking multiple derivatives, as is done with the modal states, it is necessary to separate the states into velocities and displacements. The positions are inertial and can be directly related to the modal reference frame. The angle of attack and angle of sideslip are treated as velocities because they describe the velocity vector of the aircraft.

In the time domain, the linear transformation between the modal coordinates and the mean stability axis states is non-singular. The non-dimensional stability axis states and positions are used to define a new stability axis RFA in Eq. (35).

[ ] ([ ] [ ] [ ]) ( [ ] [ ]) − 1 ∗ ∗ ∗ ∗ r r ∗ ˆ ˆ ˆ ˆ { ˜ q ( k ) } ≈ − A { ˆ x } − A + ik A + [ D ] ik I − R E { ˆ u } (35) r r 0 1 2 The primary difference from the traditional RFA in Eq. (33) is that the velocities are no longer a direct integral of the positions (i.e. { ˆ x } ik 6 = { ˆ u } ). As a result, the RFA transfer function now requires twice as 10 of 19 American Institute of Aeronautics and Astronautics many inputs or additional information, in the form of the kinematic equations, to describe the aerodynamic forces.

There are many existing codes for the traditional RFA in Eq. (33), and Roger’s RFA has been used on previous flight programs. The pedigree of the modal RFA has motivated most modeling efforts to use the modal form to create the equations of motion and then transform to a non-inertial frame as a final step.

Using this approach, the time domain transformation in Eq. (19) is applied to the traditional RFA, Eq. (33), to give the matrices (Eqs. (36), (37), (38), and (39)) for the stability axis RFA, Eq. (35).

[ ] [ ] [ ] [ ] [ ] ∗ ˆ ˆ ˆ ˆ ˆ ∗ A = A T + A T (36) 0 eh 1 e h [ ] [ ] [ ] [ ] [ ] [ ] [ ] − 1 ∗ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ∗ ∗ ∗ A = A T + A T T T (37) 1 2 eh s h e h s h [ ] [ ] [ ] ∗ ˆ ˆ ˆ A = A T ∗ (38) s h [ ] [ ] [ ] ∗ ˆ ˆ ˆ E = E T ∗ (39) s h r The [ D ] and [ R ] matrices are defined identically for both RFA forms.

r The physics of the unsteady aerodynamics require the steady state portion to be independent of the aerodynamic lag states. The indepencence conditions results in the constraint of Eq. (40).

[ ] ( [ ] [ ]) r r + ˆ ∗ ˆ [ D ] ik I − R [ E ] T { x } ≡ [ 0 ] ∀ k ∈ R (40) r r e h Additionally, the aerodynamic forces are independent of the inertial rigid body states. Numerical errors in the GAFs and curve fitting errors introduced by the RFA in general will cause the traditional modal RFA to violate these requirements. Attempting to correct these errors will negate many of the benefits of the least squares estimation used in calculating the RFA matrices. These adjustments will increase the total error in the approximation and will specifically add a bias to the approximation.

The use of the modal RFA clearly results in erroneous coefficients that have no physical basis. The total error of these coefficients is generally small, but they do inhibit user acceptance of the models by controls engineers and confuse any efforts for model tuning. The model reduction required for a low order controller is also complicated. Because the errors are at a low frequency, the steady state gain of the models cannot be maintained in the model reduction.

By applying the transformation to the exact frequency domain results, before these errors are introduced, the least squares estimates remain unbiased. The generalized modal forces are defined as a function of a ∗ stability axis GAF, [ Q ( k )] . The standard GAF are a frequency response using displacements as inputs. The new stability axis RFA is instead using velocities as an input. To ensure that the new GAF would reduce to the classical GAF under a unitary transformation, the integral of the velocities, Eq. (41), is used instead.

{ } ∗ { q ( k ) } = ¯ q [ Q ( k )] ˆ u ( k ) (41) ik [ ] ˜ The frequency domain transformation, T ( k ) , is applied to the modal GAF to provide a GAF based on the stability axis states in Eq. (42).

[ ] ∗ ˜ [ Q ( k )] = [ Q ( k )] T ( k ) (42) The rigid body states in the stability axis positions { ˆ x } represent the position and orientation of the vehicle, ∗ which has no effect on the aerodynamic forces. The coefficients for the positions in the [ A ] matrix are defined to be zero. The remaining states in the { ˆ x } vector are inertial states, so the derivatives are equal to the velocities. Therefore, the GAF at steady state ( k = 0 ) can be handled as a separate case. The application of a separate transformation at k = 0 makes the common practice of constraining the zero frequency a [ ] ∗ ˆ requirement. The constraint is achieved by selecting the A matrix in Eq. (43).

[ ] [ ] ∗ ˆ ˜ A = − [ Q (0)] T (43) 11 of 19 American Institute of Aeronautics and Astronautics Constraining the RFA at steady state is required to remove the earth fixed axis displacements from the RFA in Eq. (44).

[ ] ([ ] [ ] [ ] ) ( [ ] [ ]) − 1 ∗ ∗ ∗ 2 r r ∗ ˜ ˜ ˆ ˆ ˆ Q ( ik ) − Q (0) T = − A ik − A k + [ D ] ik I − R E ik (44) 0 r r 1 2 The RFA can be separated into is real part, Eq. (45), and imaginary part, Eq. (46).

[ ] [ ] ( ) [ ] − 1 [ ] [ ] ∗ ∗ 2 2 r r ∗ 2 ˜ ˜ ˆ ˆ Re Q ( ik ) − Q (0) T ≈ A k − [ D ] k I + R E k (45) 0 r r [ ] [ ] ( ) [ ] − 1 [ ] [ ] [ ] ∗ ∗ 2 r r r ∗ ˜ ˆ ˆ Im Q ( ik ) ≈ − A k + [ D ] k I + R R E k (46) r r r These equations can be rewritten into a matrix form in Eq. (47). The reduced frequencies were used to scale the problem to give the best fit.

    ( ) [ ] ∗ 1 ˆ ( ) ∗ A ˜ ˜ − 1 Re Q ( ik ) − Q (0) T 2 2   0 − k I k D k I + R   k ∗ ˆ ( ) = − ( ) (47)    A  − 1 1 2 2 ∗ ˜ I 0 − D k I + R R Im Q ( ik ) ∗ ˆ E k The least square solution of the overdetermined problem is given by the Moore–Penrose pseudoinverse.

D. Rational Function Approximation Augmentation The lifting surface methods being used for the aerodynamics models include assumptions of no viscous effects and a constant freestream velocity. These assumptions mean that the aerodynamic models cannot provide coefficients for drag or for velocity variations. These assumptions work well for modeling aeroelastic phenomenon. However, the flight dynamics modes, in particular the phugoid mode, are characterized by large variations in the vehicle velocity. The damping of the phugoid mode is also dominated by the vehicle drag, which necessitates the inclusion of the skin friction drag. The forces due to gravity also need to be added to the RFA force matrices. These augmentation matricies are added to the orginal RFA matricies, Eqs. (48), (49), (50), and (51).

( ) [ ] [ ] [ ] ∗ ∗ ∗ r ˆ ˆ [ A ] ≡ A + A D (48) x r 0 0 0 grav ( ) [ ] [ ] [ ] ∗ ∗ ∗ r ˆ ˆ [ A ] ≡ A + A D (49) u r 1 1 1 aug [ ] [ ] ¯ c ∗ ∗ r ˆ [ A ] ≡ A D (50) u r 2 2 2 V e [ ] [ ] 2 V e ∗ ∗ r ˆ [ E ] ≡ E D (51) u r ¯ c 1. Aerodynamic Force For incompressible flow the aerodynamic coefficients do not change with airspeed. However, the forces do change due to the dynamic pressure. The linearization of the changes in dynamic pressure, Eq. (52), is captured by the change in velocity.

∂ ¯ q ¯ q = 2 (52) ∂V V e e ∗ Since the ¯ q variations are approximated by the variations in velocity, these errors are isolated to the [ A ] ∗ matrix terms in Eq. (35). An additive correction to the [ A ] matrix, Eq. (53), is calculated from quasi-steady wind tunnel data. The coefficients in the RFA matrix are identical to the dimensional coefficients in the classical linearized rigid body dynamics. The only difference is the addition of the forces on the structural 12 of 19 American Institute of Aeronautics and Astronautics modes at trim, C .

ηh   − 2 C 0 C 0 0 0 0 · · · 0 D L 0 0   2 C 0 0 0 0 0 0 · · · 0 Y  0     − 2 C 0 − C 0 0 0 0 · · · 0  L D 0 0     − 2 bC 0 0 0 0 0 0 · · · 0 l [ ]     ∗ ˆ 2¯ cC 0 0 0 0 0 0 · · · 0 A = S m (53)  0    aug   − 2 bC 0 0 0 0 0 0 · · · 0 n     2 C 0 0 0 0 0 0 · · · 0 η 1     . . . . . . . .

.

. . . . . . . . .

  .

. . . . . . . .

2 C 0 0 0 0 0 0 · · · 0 ηh This correction for the AIC represents the minimum required to accurately capture the vehicle flight dynamics.

Many methods exist for further corrections that would improve the accuracy of the models, but a comprehensive exploration of these methods was beyond the scope of the present work.

2. Gravity Force Assuming the small displacements of the FEM, the gravitational force on each element is the linear function of the displacements in Eqs. (54) and (55).

  T 0 0 0 grav    0 T · · · 0  grav   { G } ≈ g [ M ] [Φ ] { x } (54) gg . . . 0 .

  . . . .

.

 . . .  0 0 · · · T grav   0 0 0 0 1 0    0 0 0 − cos θ 0 0      0 0 0 0 0 0   [ T ] = (55) grav   0 0 0 0 0 0      0 0 0 0 0 0  0 0 0 0 0 0 Because the orientation of the gravity vector is fixed in the inertial frame, the rotation matrix, [ T ] , gives grav the change in the forces in the non-inertial stability axis frame. The rotation matrix is identical to the form seen in traditional rigid body dynamics. The difference is that the matrix is repeated for each grid point in the FEM.

The gravitational force matrix can be combined with the aerodynamic forces by normalizing with respect to the dynamic pressure. Consistent with the aerodynamic forces, the dynamic pressure is used for nondimensionalization of the gravity in Eq. (56).

  T 0 0 0 grav   [ ]  0 T · · · 0  g grav T ∗ ˆ   A = [Φ ] [ M ] [Φ ] (56) 0 gg . . . 0 1 .

  . . . .

grav ¯ q .

 . . .  0 0 · · · T grav E. Dynamics and Kinematics The equations of motion in the modal frame with external forcing are identical to the classical flutter equations of motion, Eq. (57), from Bisplinghoff.

/ 2 [ M ] { ¨ η } + [ B ] [ K ] { ˙ η } + [ K ] { η } = { q ( t ) } (57) 13 of 19 American Institute of Aeronautics and Astronautics To simplify the expression of the equations of motion, a dimensional form of the transformation and RFA matrices are defined by Eqs. (58), (59), and (60).

[ ] [ ] [ ] r ˆ ∗ ∗ T ≡ T D (58) u r s h s h [ ] [ ] [ ] r ˆ ∗ ∗ T ≡ T D (59) x r e h e h [ ] [ ] r ˆ [ T ] ≡ T D (60) eh eh x r The equations of motion use the dimensionalized augmented RFA matrices from Eq. (48), (49), (50), and (51) are used to construct the state space model.

To determine the kinematic equations for the accelerations, e.g. { ¨ η } , the kinematic equations of Eq. (19) are differentiated to give Eqs. (61) and (62).

[ ] [ ] 2 V 2 V 0 0 { ¨ η } = T ∗ { ˙ u } + T ∗ { ˙ x } (61) s h e h ¯ c ¯ c { ˙ η } = [ T ] { ˙ x } (62) eh Using the relationship for { ˙ η } from Eq. (19), solve for { ˙ x } in the kinematic equation, Eq. (61), so that Eq. (63) has no { ˙ x } terms.

( ) [ ] [ ] [ ] 2 V 2 V 0 0 − 1 { ¨ η } = T ∗ { ˙ u } + T ∗ [ T ] T ∗ { u } eh ¯ c s h ¯ c e h s h ( ) [ ] [ ] 2 V 0 − 1 + T ∗ [ T ] T ∗ { x } (63) eh ¯ c e h e h For the current transformation, Eq. (63) can be simplified due to the condition in Eq. (64).

[ ] [ ] − 1 T ∗ [ T ] T ∗ = [ 0 ] (64) eh e h e h Therefore the complete set of dimensional kinematic equations are expressed by Eqs. (65), (66), and (67).

( ) [ ] [ ] [ ] 2 V 2 V 0 0 − 1 { ¨ η } = T ∗ { ˙ u } + T ∗ [ T ] T ∗ { u } (65) eh s h e h s h ¯ c ¯ c [ ] [ ] 2 V 2 V 0 0 { ˙ η } = T ∗ { u } + T ∗ { x } (66) s h e h ¯ c ¯ c { η } = [ T ] { x } (67) eh The equivalent mass (Eq. (68)), damping (Eq. (69)), and stiffness (Eq. (70)) matrices for the stability axis are derived by substituting these kinematic equations into the modal equations of motion, Eq. (57).

[ ] [ ] 2 V ∗ ¯ M ≡ [ M ] T ∗ + ¯ q [ A ] (68) s h ¯ c ( ) [ ] [ ] [ ] 2 V 2 V / 2 0 − 1 0 ¯ B ≡ [ M ] [ B ] [ K ] + T ∗ [ T ] T ∗ (69) eh e h s h ¯ c ¯ c [ ] ¯ K ≡ [ K ] [ T ] (70) eh These matrices are used to assemble the state space equations, Eq. (71).

      2 V − 1 2 V − 1 0 e   T T ∗ T T ∗ 0   ˙ x x   ¯ c eh ¯ c eh   e h s h ( ) ( )   − 1 ∗ − 1 ∗ − 1 = ¯ ¯ ¯ ¯ ¯ ˙ u   u − M K + ¯ q A − M B + ¯ q A − ¯ q M D 0 1         2 V ∗ 0 r ¯ ˙ x x 0 E − R a a r ¯ c       0 0 0 δ   ( )   − 1 ∗ − 1 ∗ − 1 ∗ ¯ ¯ ¯ ˙ + (71)  − ¯ q M A − ¯ q M A − M M + ¯ q A  δ hc hc hc hc 0 1 2     ∗ ˜ ¨ 0 E 0 δ hc 14 of 19 American Institute of Aeronautics and Astronautics III. Results The X-56A MUTT (figure 2) was specifically developed for studying the modeling and control of a vehicle with body freedom flutter. The current flight-test data with the stiff wings does not have an unstable flutter mode, but does demonstrate significant measurable coupling between the flight dynamics and the structural dynamics. Flight-test data for a wide range of fuel loads and airspeeds are available. The higher airspeeds offer a better demonstration of the aeroelastic coupling. The lower fuel cases will identify any potential inaccuracies in the assumed modes method. The higher fuel cases demonstrate the effects of a statically unstable vehicle. Given the wide variety of available flight-test data, the stiff wing flight data provide a good data set for validation of the current modeling methodology. For modeling the X-56A MUTT, the modes shapes from the FEM for the full fuel case were found to be sufficient for the assumed modes.

Figure 2. X-56A MUTT in flight.

The ZAERO ™ (ZONA Technology Incorporated, Scottsdale, Arizona), subsonic lifting surface method was selected to provide the frequency domain GAF, [ Q ( k )] in Eq. (32), for the current modeling. The lifting surface method in ZAERO ™ is helpful because it avoids many of the numerical singularities that are present in more traditional doublet lattice methods. As a result, the mesh is much less sensitive to mesh refinement.

The cost is the requirement for a finer chordwise mesh to correctly capture the lift distribution. ZAERO ™ also offers many utilities to aid in the modeling. Although ZAERO ™ has been used for the present modeling, the current proposed modeling approach (e.g. stability axis RFA) is directly applicable to any frequency domain aerodynamic model such as doublet lattice and could be extended to apply to a frequency domain model derived from higher fidelity computational fluid dynamics (CFD) models.

The pitch rate response of the new integrated model is compared to an early rigid body 6 DoF model in figure 3. The frequency responses are at a low airspeed and low fuel case. The rigid model does have corrections to account for the static deflections. The frequency and damping of the short-period and phugoid mode match very well for both models. The models begin to deviate as the frequency approaches the first wing bending mode.

15 of 19 American Institute of Aeronautics and Astronautics Rigid Flexible Magnitude, dB Phase, deg Frequency, Hz Figure 3. Comparing pitch rate due to wing flap 4 of classical rigid model and new integrated model.

The models have also been compared against the flight test data. A summary of the test cases compared against the flight data are shown in table 1. The longitudinal dynamics are identified from a multisine sweep from the pitch allocator in the control system. The first case is at low speed and low fuel weight. The model matches well to the identified frequency response of the pitch rate gyro and an accelerometer on the right wing tip in figure 4. The frequency of the first wing bending mode matches extremely close, even though the mode shapes were taken from the full fuel case. The accurate frequencies show that a sufficient number of modes are included into the assumed modes solution to correctly capturing the structural frequencies.

Table 1. Validation test cases.

Test case Fuel mass Airspeed Input 1 Low Low Pitch 2 High Low Pitch 3 Low High Pitch 4 Low High Roll 16 of 19 American Institute of Aeronautics and Astronautics Short-period First wing Flight test Flight test bending Model Model Magnitude. dB Magnitude. dB Phase, deg Phase, deg Coherence Coherence Frequency, Hz Frequency, Hz (a) Pitch rate due to pitch allocator. (b) Right outboard accelerometer due to pitch allocator.

Figure 4. Low fuel, low speed pitch allocator sweep.

The accelerometers are a higher bandwidth sensor, so they are showing better coherence across the bandwidth of the frequency sweep. There is a sharp drop in the coherence after the highest frequency used in the multisine because there is no power being input at those frequencies.

The same type of sweep at a high mass case still matches well with the models in figure 5. The effect of the turbulence is reduced slightly at higher fuel weights, so the coherence is improved slightly for the high fuel weight case. This effect is most notable at higher frequencies, where the magnitude of the response is lower and thus more susceptible to noise from turbulence. The high mass cases shows the phase in the accelerometer output at the second mode more clearly.

Flight test Flight test Model Model Magnitude. dB Magnitude. dB Phase, deg Phase, deg Coherence Coherence Frequency, Hz Frequency, Hz (a) Pitch rate due to pitch allocator. (b) Right outboard accelerometer due to pitch allocator.

Figure 5. High fuel, low speed pitch allocator sweep.

The test case shown in figure 6 is another low fuel pitch allocator sweep, but at a high airspeed. At these airspeeds the aerodynamic lags are more important for correctly capturing the aeroelastic coupling. Issues with RFA would begin to be more apparent. At these higher airspeeds and low fuel cases where the structural coupling is most prevalent, the matching between the models and the flight-test is very good.

17 of 19 American Institute of Aeronautics and Astronautics Flight test Flight test Model Model Magnitude. dB Magnitude. dB Phase, deg Phase, deg 1 1 Coherence Coherence 0 0 Frequency, Hz Frequency, Hz (a) Pitch rate due to pitch allocator. (b) Right outboard accelerometer due to pitch allocator.

Figure 6. Low fuel, high speed pitch allocator sweep.

The final test point is a multisine sweep to the roll allocator. The frequency response of the roll rate gyro and the same wing tip accelerometer match well to the model shown in figure 7. The roll axis does not show the same dependence of fuel weight, so only one case is shown.

Flight test Model Flight test Model Magnitude. dB Magnitude. dB Phase, deg Phase, deg 1 1 Coherence Coherence 0 0 Frequency, Hz Frequency, Hz (a) Roll rate due to roll allocator. (b) Right outboard accelerometer due to roll allocator.

Figure 7. Low fuel, high speed roll allocator sweep.

The frequency response of the roll axis matches very well between the flight-test data and the model. The individual modes do not have as large of a magnitude peak in the roll response, but the frequencies of the modes can been identified in the phase crossover frequencies. The higher roll inertia, relative to pitch, reduces the influence of the turbulence. As a result, coherence for the roll response is much higher than for the pitch cases.

IV. Conclusion A method has been demonstrated for converting the unsteady aerodynamics in a modal coordinate system, generated by many legacy aeroelastic modeling tools, to a stability axis system. This approach is consistent with the methodology required for flight control design. An assumed modes method was selected for describing 18 of 19 American Institute of Aeronautics and Astronautics the structural dynamics. This method allows the use of the same mode shapes to ensure consistency in the definition of the states across the flight envelope. Consistency of the model states is required for gain scheduling controllers to maintain stability across the flight envelope. These methodologies were applied to the X-56A MUTT and were compared against flight-test data for the stiff wing configuration. The models generated by the present approach are able to accurately capture the coupled flight dynamics and structural dynamics of the X-56A MUTT.

References Love, M. H., Zink, P. S., Wieselmann, P. A., and Youngren, H., “Body Freedom Flutter of High Aspect Ratio Flying Wings,” AIAA 2005-1947, 2005.

Bisplinghoff, R. L., Ashley, H., and Halfman, R. L., Aeroelasticity , Dover Publications, Inc., Cambridge, 1996.

Milne, R. D., “Dynamics of the Deformable Aeroplane,” R.&M. No. 3345, 1964.

Noll, T. E., Brown, J. M., Perez-Davis, M. E., Ishmael, S. D., Tiffany, G. C., and Gaier, M., “Investigation of the Helios Prototype Aircraft Mishap,” NASA, 2004.

Stevens, B. L. and Lewis, F. L., Aircraft Control and Simulation , John Wiley and Sons, Inc., Hoboken, 2003.

Beranek, J., Nicolai, L., Buonanno, M., Burnett, E., Atkinson, C., Holm-Hansen, B., and Flick, P., “Conceptual Design of a Multi-Utility Aeroelastic Demonstrator,” AIAA 2010-9350, 2010.

Waszak, M. R. and Schmidt, D. K., “Flight Dynamics of Aeroelastic Vehicles,” Journal of Aircraft , Vol. 25, No. 6, June 1988, pp. 563 – 571.

MSC Nastran 2013 Dynamic Analysis User’s Guide , MSC.Software Corporation, Santa Ana, CA, 2013.

Meirovitch, L., Fundamentals of Vibrations , McGraw-Hill, Inc., Boston, 2001.

Baldelli, D. H., Chen, P. C., Panza, J., and Adams, J., “Unified Rational Function Approximation Formulation for Aeroelastic and Flight Dynamics Analyses,” AIAA 2006-2025, 2006.

Roger, K. L., Hodges, G. E., and Felt, L., “Active Flutter Suppression – A Flight Test Demonstration,” Journal of Aircraft , Vol. 12, No. 6, June 1975, pp. 551 – 556.

Roger, K. L., “Airplane Math Modeling Methods for Active Control Design,” AGARD-CP-228, North Atlantic Treaty Organization. Advisory Group for Aerospace Research and Development., 1977, pp. 4.1 – 4.11.

Tiffany, S. H. and Karpel, M., “Aeroservoelastic Modeling and Applications of Using Minumum-State Approximations of the Unsteady Aerodynamics,” NASA TM-101574, 1989.

19 of 19 American Institute of Aeronautics and Astronautics