Generation and Calibration of Linear Models of Aircraft with Highly Coupled Aeroelastic and Flight Dynamics

Generation and Calibration of Linear Models of

Aircraft with Highly Coupled Aeroelastic and Flight

Dynamics

∗ † Jeffrey Ouellette and Felipe D. Valdez NASA Armstrong Flight Research Center, Edwards, California, 93523 The lightweight structures and unconventional configurations being considered for the next generation of aircraft mean that any effort to predict or control the flight dynamics is impacted by the structural dynamics. 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 a relatively weak effect of structural dynamics on the flight dynamics, and are therefore incapable of modeling strong interactions like body freedom flutter. A method of combining different sources of data traditionally used for aeroelasticity and flight dynamics is described by reconciling many of the differences between these models.

By building upon past modeling efforts, a level of familiarity in the approach is achieved.

Generally the differences from the traditional approaches are subtle but significant. The traditional frequency domain flutter model in a modal coordinate system is converted to a form consistent with a time domain flight dynamics model. The time domain rational function approximation about a non-inertial coordinate system and the unique constraints for the conversion between the inertial and non-inertial coordinate systems are discussed. A consistent transformation of the states of aeroelastic models to flight dynamics models is derived, which enables the integration of data from higher fidelity computational fluid dynamics models or wind-tunnel testing. The present method of integrating multidisciplinary data was used to create models that compare well with X-56A flight-test data, including conditions past the flutter speed.

I. Nomenclature A steady aerodynamic coefficients A quasi-steady aerodynamic coefficients A added mass aerodynamic coefficients AIC aerodynamic infuence coefficient matrix AIC aerodynamic influence coefficient argmin argument of the minimum (value which minimizes the function) ASE aeroservoelastic ¯ B effective stability axis damping matrix B structural damping matrix b wing span ¯ C aerodynamic panel pressure at reference condition p ¯ c mean aerodynamic chord C rotation matrix C aerodynamic coefficients to match match C aerodynamic drag force coefficient at trim D C aerodynamic lift force coefficient at trim L ∗ Aerospace Engineer, Controls and Dynamics Branch, P.O. Box 273/Mailstop 4840D, Edwards, California, 93523, AIAA Member.

† NASA Pathways Intern, Controls and Dynamics Branch, P.O. Box 273/Mailstop 4840D, Edwards, California, 93523.

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 CFD computational fluid dynamics C aerodynamic hinge moment due to angle of attack hm α C lift curve slope L α C pitch stiffness m α C aerodynamic wing bending load due to angle of attack SW 1 B α C aerodynamic wing twist load due to angle of attack SW 1 T α C aerodynamic modal force coefficient at trim ηh D aerodynamic lag output matrix D velocity non-dimensionalization matrix u D position non-dimensionalization matrix x d distance from node k to the control point of panel j kj DLM doublet lattice method DoF degrees of freedom E aerodynamic lag input matrix ¯ F steady aerodynamic downwash matrix F ( k ) aerodynamic downwash matrix ˜ F unsteady aerodynamic downwash matrix f downwash on aerodynamic panel j j FE finite element G vector of forces on each grid point due to gravity g acceleration due to gravity GAF generalized aerodynamic forces H low order equivalent system transfer function loes H transformation between Euler angle rates and body axis angular rates I identity matrix input I coefficient input incidence matrix ha output I coefficient output incidence matrix ha I aerodynamic panel-node incidence matrix kj Im () imaginary part i imaginary number j index of polynomial term ¯ 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 L non-dimensionalizing length LOES low order equivalent system ¯ 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 n transfer function numerator coefficient i n unit vector normal to aerodynamic panel j j n number of polynomial terms NED north-east-down coordinate system alignment p roll rate ¯ q dynamic pressure ˆ Q GAF polynomial approximation ˆ Q j -th coefficient matrix of GAF polynomial approximation j Q generalized aerodynamic force matrix q generalized modal forces ˜ Q assumed mode generalized aerodynamic force matrix ˜ q approximation of the generalized modal forces q pitch rate R aerodynamic lag pole matrix r position vector Re () real part r yaw rate RFA Rational function approximation s laplace variable S aerodynamic panel area matrix S Wing reference area ˆ T ∗∗ ∗ transformation body axis rate to modal acceleration h b ˆ T ∗∗ transformation body axis to modal acceleration h b ˆ ∗ T transformation body axis to modal rate hb ˆ T ∗ transformation Earth-fixed axis to modal rate he ˆ T transformation Earth-fixed axis to modal displacement he T transformation matrix of gravity components in the stability axis grav ˜ T ( k ) frequency domain transformation from body coordinates to modal coordinates ˜ T steady-state transformation from body coordinates to modal coordinates T kinetic energy t time T nonlinear transformation modal states and body axis accelerations accel T nonlinear transformation modal states and body axis displacements disp T nonlinear transformation modal states and body axis rates rate U spline from aerodynamic panels to finite element model grid u local air velocity in the Eulerian perspective at aerodynamic panel j j,j v velocity vector V local flow velocity at aerodynamic panel j L V potential energy V reference airspeed (true airspeed at trim) V equivalent airspeed e v flexible-body deflection rates f lex v rigid-body velocities and angular rates rigid v wind-axis states describing the body axis w W AIC correction factor weights W AIC correction factor weights for modal inputs mode W AIC correction factor weights for control surface inputs surface x position of aerodynamic panel j j x earth-fixed longitudinal position x aerodynamic lag states aero x flexible-body deflections f lex x rigid-body positions and orientations rigid y earth-fixed lateral position z earth-fixed vertical position Subscripts 0 quantity at trim a analysis (a-set) degrees of freedom aug augmentation b body (mean) axis system c control surface (c-set) degrees of freedom e earth-fixed coordinate system g finite element grid (g-set) degrees of freedom grav contribution of gravity h modal (h-set) degrees of freedom j aerodynamic box grid degrees of freedom k aerodynamic box control point degrees of freedom w wind-axis system Superscripts ∗ stability axis rational function approximation Symbols α angle of attack α local angle of attack of j -th panel j β angle of sideslip γ flight path angle δ control surface inputs ¯ ε AIC correction factor bias vector ¯ ε j -th AIC correction factor bias j ε AIC correction factor vector ε j -th AIC correction factor j ζ damping ratio of i -th mode i η displacement state vector ˙ η velocity state vector θ pitch angle λ smoothing contraint parameter hp ξ modal coordinate Φ assumed mode shape φ roll angle ω angular rate vector ψ heading angle ω natural frequency of i -th mode i Conventions D total/material derivative Dt ˆ () non-dimensionalized T matrix transpose [ ] matrix r [ · ] diagonal matrix r ‖‖ vector norm { } vector C rotation from frame a to frame b b/a r vector from point a to b expressed in reference frame c.

b/a,c ∇ jacobian with respect to state vector η η ∗ () derivative with respect to non-dimensional time ˙ () derivative with respect to time ¨ () second derivative with respect to time := equal by definition ≡ identically equivalent ∀ for all II. 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. One classic example is the General Dynamics (Falls Church, Virginia) RB-57H Canberra[ 1 ]. In an effort to improve high altitude performance, the RB-57H, shown in Fig. 1, was modified with a larger wing. The wing being larger relative to the tail, the damping of the short-period mode was decreased. The larger wing was also inherently more flexible. The interaction between the reduction of short-period damping and the increase in wing flexibility caused instances of body freedom flutter.

Fig. 1 The RB-57H Canberra (General Dynamics, Falls Church, Virginia).

The use of an active flight control system offers an appealing solution for suppressing body freedom flutter instability. Developing a flutter suppression control system for body freedom flutter requires a model that considers flight dynamics, structural dynamics, and interactions. 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 discipline. 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 approach to model integration that will address these conflicts.

The X-56A (Lockheed Martin, Bethesda Maryland) Multi-Utility Technology Testbed (MUTT) shown in Fig. 2 was developed specifically to study active control of the structure of an aircraft and suppression of body freedom flutter[ 2 ]. The wings of the aircraft can be swapped to fly an aircraft with a relatively stiff structure or a flexible aircraft with unstable flutter modes. Unlike most production aircraft today, the X-56A MUTT has structural dynamics that are catastrophically coupled with the flight dynamics; with the flexible wings, the coupling is strong enough to destabilize the aircraft dynamics.

Fig. 2 The X-56A MUTT in flight.

The ability to operate a vehicle with body freedom flutter like the X-56A MUTT requires that the models be able to accurately predict the flutter and be suitable for generating control laws to stabilize the flutter mode. The field of aeroelastic analysis has focused on the prediction of instabilities in the structural dynamics, but is of limited accuracy for rigid-body dynamics. Flight dynamics models use higher fidelity tools to obtain accurate models of the rigid-body dynamics, but do not consider structural effects. The objective of the present work is to integrate the results of existing flutter models and flight dynamic models to create a model that is capable of accurately predicting the in-flight dynamics of both rigid-body and flexible-body dynamics and that is useful for synthesis of a control system.

Many of the modern techniques for control law design require a model in a time domain state-space form like the traditional rigid-body model in Eq. (1a). The Helios (Aerovironment, Monrovia, California), shown in Fig. 3, mishap investigation showed the importance of a time domain form of the models for evaluating an aircraft closed-loop performance[3].

Fig. 3 The helios (Aerovironment, Monrovia, California) in flight.

The models and control laws are typically generated at discrete points, but the aircraft will operate within a continuously varying flight envelope. As a result of the continuous flight envelope, the implementation of control laws for flight requires the ability to interpolate between the grid points of the models. Such interpolation requires that the states have a consistent meaning at each flight condition. Also, the complexity of the control laws often is dependent on the order of the model. If the control law is too complex, it will not be possible to implement on a real-time flight computer. The real-time execution of the control laws therefore require the models be of relatively low order or be able to be reduced to a low order. The classic rigid-body model, Eq. (1a), consists only of six degrees of freedom (DoF)[ 4 ], rigid-body orientation/position, x , and rigid velocities/rates, v . The six-DoF models are of a sufficiently low order for most flight control applications.

rigid The phugoid (a low frequency, lightly damped, and large magnitude oscillation in the vehicle altitude and airspeed) mode is primarily effected by the velocity variations and gravity effects. Stability being a property of a trim condition, using linear models for evaluation of dynamics that include velocity DoF require that the axis translate with the vehicle. The translating axis means that the models are linearized about the trim; however, inclusion of gravity requires a coordinate system fixed to the Earth. The presence of translating and fixed axes needed to model the phugoid mode results in non-inertial kinematics. The non-inertial kinematics mean that the vehicle velocity is not equal to the derivative of the position. The non-inertial reference frame ˙ results in more complex kinematic equations relating v and x to x in Eq. (1a).

rigid rigid rigid The traditional flutter model of Eq. (1b) adds the flexible states, x and v . The kinematic f lex f lex equations for flutter models are simply identity matrices; however, the parameters change with reduced frequency, k . The frequency dependence of the coefficients makes these flutter models unsuitable for control law development. Typically, aeroelastic models like Eq. (1b) are created for a fixed airspeed[ 5 ]; however, removing the velocity makes it impossible to model the phugoid mode.

Previous efforts for flutter suppression[ 6 , 7 ] made aeroservoelastic (ASE) models of the form of Eq. (1c) which added additional aerodynamic lag states, x , to remove the frequency dependence of the coefficients.

aero The ASE models in Eq. (1c) retain the simplistic kinematic relations of the flutter, and typically lack a velocity degree of freedom.

       ˙ x  [ · · · ] [ · · · ]  x  rigid rigid                           =   (1a) ˙ v [ · · · ] [ · · · ] v rigid rigid                                    ˙ x  0 0 I 0  x  rigid rigid                   ˙ x 0 0 0 I x  f lex     f lex    ˙ =   (1b) v [ · · · ] ( k ) [ · · · ] ( k ) [ · · · ] ( k ) [ · · · ] ( k ) v rigid rigid                 ˙ v [ · · · ] ( k ) [ · · · ] ( k ) [ · · · ] ( k ) [ · · · ] ( k ) v     f lex f lex               ˙ x 0 0 I 0 0 x     rigid rigid                   ˙ x 0 0 0 I 0 x  f lex     f lex    = (1c) ˙ v  [ · · · ] [ · · · ] [ · · · ] [ · · · ] [ · · · ]  v rigid rigid                 ˙ v [ · · · ] [ · · · ] [ · · · ] [ · · · ] [ · · · ] v     f lex f lex         ˙ x 0 0 [ · · · ] [ · · · ] [ · · · ] x aero aero       ˙ x [ · · · ] 0 [ · · · ] 0 0 x     rigid rigid                   ˙ x 0 0 0 I 0 x  f lex     f lex    = (1d) ˙ v  [ · · · ] [ · · · ] [ · · · ] [ · · · ] [ · · · ]  v rigid rigid                 ˙ v 0 [ · · · ] [ · · · ] [ · · · ] [ · · · ] v  f lex   f lex          ˙ x 0 0 [ · · · ] [ · · · ] [ · · · ] x aero aero The inconsistency in kinematic equations between Eqs. (1a) and (1c) complicates the integration of data between these models, and produces erroneous behavior. By returning to first principles a method has been developed to go directly from the flutter model of Eq. (1b) to an model of the form of Eq. (1d), which is more suitable for the integration of data used in the flight dynamics model of Eq. (1a). The resulting models are capable of capturing the requirements of both the flight controls engineers and the aeroelasticity engineers.

A linear finite element (FE) model is used for the structural dynamics. In modern aeroelasticity the vehicle structure is discretized into many finite elements. Because the primary concern is the deformations of the structure, these FE models use an inertial coordinate system. The many thousands of DoF of a FE model often make direct simulation impractical. For aeroelastic analysis, the order of the structural dynamics model is reduced by considering only a subset of the vehicle structural modes, which includes the rigid modes[ 8 ].

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.

The aerodynamics model is based primarily on a linear frequency domain model. The traditional lifting surface methods from flutter analysis were selected for evaluation, but the methods of integration being discussed are extensible to other higher fidelity time linearized tools. Due to the higher frequencies of the structural dynamics, the rate of response of the aerodynamics is on the same order as the structural dynamics. To capture the dynamics of the airflow over the vehicle as the structure deforms requires an unsteady aerodynamic model. These unsteady models include the airflow changing over the aircraft and traveling into the wake. Often many of the tools used for aeroelasticity have been based on a frequency domain approach; however, these unsteady aerodynamic tools developed for flutter are often inaccurate for rigid-body modes, and do not give any measure of drag. 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. To improve the accuracy of the rigid-body modes, data from computational fluid dynamics (CFD) or wind-tunnel data are integrated into the aerodynamic model.

To resolve the conflicts between the aeroelastic and flight dynamics modeling, the current work has returned to the fundamental principles shared by both of these disciplines in order to develop a new approach to integration of models that is consistent with both applications. An effort has been made to keep the integrated models similar to the existing models, so that superficially they appear very similar to the current state of the art. Addressing the conflicts between the models from the independent disciplines has enabled a number of incremental improvements in the models.

To describe the modeling approach, first the coordinate systems used for each discipline are defined. The frequency domain aerodynamic model and structural dynamics are translated into a form compatible with time domain states space models. The time domain state space model generated is compatible with methods for control laws and is more consistent with the aerodynamic models seen in flight dynamics. Augmentations to this aerodynamic model are described to address the limitations of the linear aerodynamic tools such as a fixed velocity and lack of gravitational forces. Two methods of integrating information from CFD or wind-tunnel data into these models to improve their accuracy are introduced.

The approach has been applied to the X-56A MUTT and is compared against the available flight-test data at conditions before and after the flutter speed.

III. Elements of the Models Most flight control development requires the plant models to have a small number of states. To achieve these requirements, a state space representation of the models was selected. The importance of higher frequency dynamics means 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, typically phugoid and spiral modes, is very important if the vehicle has low or negative static margin. Low static margin can cause a coupling between the statically unstable short-period and the phugoid mode. The coupling between rigid and flexible dynamics will result in an unstable phugoid-like mode, which 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.

The discussion of the model elements is broken down by traditional disciplines (structures and aerodynamics).

The form of the structural dynamics is discussed first. Using the knowledge of the structural dynamics, an aerodynamic model is derived that includes both these structural dynamics and rigid-body motions.

A. Coordinate Systems To describe the current modeling approach, first the coordinate systems and their relationships are defined and discussed. The mathematical translation between these coordinate systems are defined in the appendix, Eqs. (A-19a) and (A-19b). The aircraft system-level reference frames shown in Fig. 4 describe the total motion of the vehicle flight dynamics and structural dynamics as a complete system.

B o d y a x i s x b y b z b r b/e r b/h z h x h M o d a l E a r t h r y e/h a x i s h a x i s x y e e z e Fig. 4 System-level reference frames.

The definition of these coordinate systems is not new to the present work; however, a more thorough discussion of their definition and their relationship in the context of a flexible aircraft can be instructive.

For a flat non-rotating Earth, the Earth-fixed coordinate system provides an inertial reference frame[ 4 ]. 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 modal coordinate system is traditionally used in aeroelasticity modeling. In the modal coordinate system, the structural deformations are described as a linear combination of a finite set of mode shapes. The use of the modal frame in existing tools makes it desirable as a starting place for the flight dynamics models; however, the modal frame is inadequate for describing velocity variations[ 5 ]. The orientation of the modal coordinate system differs from the Earth-fixed coordinate system, and matches the conventions typically used in aeroelastic modeling. Like the axis of the Earth, the modal coordinate system is also an inertial reference frame. The aeroelastic modeling tools define the modal frame such that it translates relative to the Earth-fixed coordinates at a constant rate. The orientation of the modal frame is defined such that the orientation relative to the Earth frame does not change. The translation being constant and the orientation unchanging implies that the modal frame is an inertial frame. At trim the modal frame is aligned with the wind-axis. As the flight path angle changes, the modal frame does not remain aligned with the wind-axis.

The body axis system in Fig. 4 is generalized to a flexible vehicle by the use of a mean axis system[ 5 , 9 ].

Like the traditional body axis, the mean body axis is fixed at the center of gravity. Because the vehicle deforms, the center of gravity is not a fixed point on the aircraft. As a result, the mean axis can be more difficult to visualize because the reference point floats as the vehicle deforms. The mean axis is beneficial because it simplifies the equations of motion by decoupling the kinetic and potential energy of the rigid-body and structural dynamics.

The aerodynamic forces are most physically described in a wind-axis system that is oriented with the relative wind, shown in Fig. 5.

α
x b y α w x w Relative Wind y b z z w b Fig. 5 Wind and body (mean) reference frames.

The relationships between these system-level coordinates are summarized in Table 1.

Table 1 System-level reference frame relationships.

:= 0 := := 0 (a) Rotations for level trim. ( φ , θ α , β ) 0 0 0 0 Earth (NED) Body Wind Modal Earth (NED) - Body (mean) φ , θ , ψ - Wind α , β - Modal aligned φ , α , β - (b) Translations for level trim.

Subscript Earth (NED) Body Wind Modal Earth (NED) e - Body (mean) b x , y , z - Wind w x , y , z coincident - Modal h x , y , z ( x − x ) , ( y − y ) , ( z − z ) - h h h h h h To clarify the reference frames and the notation for the vectors, the example of the position vector r is b/e shown in Fig. 4. The notation of r means motion of the origin of frame “ b ” relative the origin of the frame b/e “ e ”. Because these are non-inertial frames, it is also necessary to describe the reference frame used to describe the components of the vector. An additional subscript is added to represent the reference frame. For the example, the vector r is described in the “ w ” reference frame.

b/e,w The vehicle’s velocity vector is always aligned with the wind x-axis. Therefore, the total velocity and the orientation of the wind-axis defined by angle of attack and angle of sideslip, α and β , respectively are sufficient for an alternative expression for the the velocity of the body axis, v , in Eq. (2c). The remaining w rigid-body states are described by the vectors in Eqs. (2a), (2b) and (2d).

{ } r = (2a) x y z b/e,b { } θ = (2b) φ θ ψ b/e,b { } v = (2c) V β α w { } ω = (2d) p q r b/e,b The rigid-body states are combined with the modal states, ξ , to complete vectors of positions, Eq. (3a), and rates, Eq. (3b)[4].

{ } T T T T { η } = (3a) r θ ξ b b/e,b b/e,b { } T T T T ˙ { ˙ η } = v ω ξ (3b) b w b/e,b For a linear system, orthogonality of the rigid-body modes and the structural modes is sufficient to define the mean axis.

The coordinate systems discussed so far are used to describe the complete vehicle system. Due to the distributed nature of the structure, is is also necessary to define local coordinate systems that describe the dynamics of discrete points on the vehicle. These local coordinate systems are grouped into sets described in Table 2. The grid points of the FE model describe the g -set of local DoF[10]. There may be many different coordinate systems used in the FE model. The definition of these coordinate systems is arbitrary and at the discretion of the engineer developing the FE model. The determination of these coordinate systems is beyond the scope of the current work, but knowledge of their definition from the FE model is required. The g -set coordinates are related to the modal coordinates, h -set, by the mode shapes. Most finite element software can provide these mode shapes. The other local coordinates are the control surface, c -set, deflections.

Table 2 Local coordinate system degree of freedom sets.

Set Name Definition a Analysis set of DoF c Control surface DoF g FE grid DoF h Modal DoF j Aerodynamic box grid DoF k Aerodynamic box control point DoF To allow a consistent aerodynamic model across the flight envelope, it is desirable to define a set of non-dimensional states for use in the aerodynamic model. The non-dimensionalization is defined in Eqs.(A-2a) ˆ and (A-2b) in the appendix. The individual states are grouped together into non-dimensional position, { η } b { } ∗ in Eq. (4a), and rates, ˆ η in Eq. (4b).

b { } T T ˆ { η } = (4a) ˆ x ˆ y ˆ z φ θ ψ ξ b { } { } T ∗ T ∗ ˆ η = ˆ (4b) b V β α ˆ p ˆ q ˆ r ξ In general, the process of adding control surface deflections, δ , can represented by adding additional modes that correspond to a unit deflection of each surface. Due to the similarity of the process of handling control surface deflections, in the present work they have generally be combined with the vector of structural modes, ξ .

B. Structural Dynamics The standard orthonormal mode shapes from an FE model have the benefit of being orthogonal with respect to both the FE mass and stiffness matrices. Ensuring the orthogonality becomes a challenge when the mass of the vehicle changes, such as during fuel burn. To maintain orthogonality, the mode shapes must change. If the models are to be used for control design, it is necessary to be able to interpolate the models between these conditions. The feasibility of the interpolation requires the modes to have a consistent definition across all flight conditions. The matching of orthogonal mode shapes thus becomes a subjective exercise that can only be partially automated. To be more robust, assumed mode shapes, [ Φ ] , are used, which do not change with fuel weight. 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. To reproduce the FE model modes, the assumed modes method will typically require more mode shapes.

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 can be challenging and is heavily dependent on attributes of the vehicle. The use of assumed modes will reflect an increase in the difficulty of generating the structural model. For more complex aircraft configurations, there can be modes that will only appear at specific mass cases, such as panel modes. These mass case specific modes must be included in the set of assumed modes. At cases where the extra modes would not normally appear, the extra mode will have a frequency much higher than any remaining modes of interest.

[ ] [ ] ˜ ˜ For these assumed modes, the mass, M , and stiffness, K , matrices are specified to capture the hh hh kinetic and potential energy of the structure given the assumed mode shapes. This approach is often referred to as the assumed modes method[ 8 ]. Using the full FE mass, [ M ] , and stiffness matrices, [ K ] , kinetic gg gg energy in Eq. (5a), T , and potential energy in Eq. (5b), V , are a function of the mode shape, modal coordinate, { η } , and the FE mass and stiffness matrices.

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

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

[ ] T ˜ M = [Φ ] [ M ] [Φ ] (7a) hh 0 gg 0 [ ] T ˜ K = [Φ ] [ K ] [Φ ] (7b) hh 0 gg 0 [ ] [ ] ˜ ˜ The mass and stiffness matrices, M and K , in Eqs. (7a) and (7b) are no longer diagonal.

hh hh C. Aerodynamics Beyond the obvious need to include additional DoF beyond typical flight dynamics models as inputs to the aerodynamic model, the modeling of aeroelasticity also requires consideration of the dynamics of the air mass, or unsteady aerodynamics. Traditionally, flight dynamics models have utilized quasi-steady aerodynamics. As the frequency of the inputs to the aerodynamic model increases, as is seen in aeroelastic models, it becomes necessary to include an unsteady 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.

The approach to generating an unsteady aerodynamic model begins with describing the classical unsteady approaches, defined for oscillations at a specific frequency, that are used in flutter analysis. The relationship and transformation between these flutter models and a classical quasi-steady flight dynamic model is demonstrated. A simple polynomial approximation to the unsteady flutter model is defined, which is useful in fitting a more sophisticated approximation. To translate these to generate motions in time, an approximation approach model is described. The method of fitting the approximation model is then discussed to accurately reproduce the underlying flutter model data.

By linearizing the aerodynamics with respect to time, it is possible to describe the aerodynamics by a frequency dependent aerodynamic influence coefficient (AIC) matrix. The columns of the AIC matrix are the linear relationship of the load on all the other surface elements induced by the harmonic oscillations of a surface element. There are specialized higher fidelity CFD codes, such as ZEUS ™ (ZONA Technology Incorporated, Scottsdale, Arizona) or the German Aerospace Center “Triangular, Adaptive and Unstructured” CFD code, that are also able to directly produce these AIC matrices. These CFD methods remain computationally expensive, so most flutter analyses simplify the geometry to a flat lifting surface and use linearized potential flow methods. The doublet lattice method (DLM)[ 12 ] and ZAERO ™ (ZONA Technology Incorporated, Scottsdale, Arizona) are two common implementations of this approach that can provide these AIC matrices in a modal reference frame. The state transformation in the frequency domain, Eq. (A-23b) derived in the appendix, is used to transform the unsteady aerodynamic model inputs to be consistent with traditional flight dynamic states. A polynomial model is presented that is useful for generating an approximation of the quasi-steady portion of the model. An approximation method is shown that allows the transformed frequency domain model to be expressed in the time domain.

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. These methods give a linear relationship in Eq. (8) called the generalized aerodynamic forces, GAF, [ Q ( k )] , that relates the modal states, η ( k ) , to the generalized modal forces, { q ( k ) } , at a specified reduced frequency, k .

h { q ( k ) } = ¯ q [ Q ( k )] { η ( k ) } (8) h The GAF shown in Eq. (9) is a function of the mode shapes, the panel areas, and the velocity of the air perpendicular to the panel (downwash). Additionally a spline matrix, [ U ] , is used to relate the mode shapes of the FE grid, to the points of the aerodynamic panels.

[ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] T T r ˜ Q ( k ) = Φ U S AIC ( k ) F ( k ) U Φ (9) 0 r 0 The assumed modes ensure that the mode shapes, [ Φ ] , are constant. Therefore, the GAF in Eq. (9) does not change with the fuel weight.

To simplify the integration with FE models, most of these unsteady aerodynamic methods generate the aerodynamics as a function of the modal axis states, η . Although the modal axis states are convenient, h they are not physically representative and are inconsistent with the rigid-body states typically seen in flight dynamics models.

[ ] ∗ To address the inconsistency, a body axis GAF, Q ( k ) , is used. The standard GAF is a frequency response using displacements as inputs. The new body axis GAF is instead using the integral of the velocities as an input. The integral of the velocities, Eq. (10), ensures that the new GAF will reduce to the classical GAF under a unitary transformation.

{ } [ ] ∗ ∗ 1 { q ( k ) } = ¯ q Q ( k ) ˆ η ( k ) (10) b ik [ ] ˜ The state transformation in the frequency domain, T ( k ) , described in Eq. (A-23b) in the appendix, is applied to the modal GAF to provide a GAF with body axis states, ˆ η , as inputs as shown in Eq. (11).

b [ ] [ ] ∗ ˜ Q ( k ) = [ Q ( k )] T ( k ) (11) { } { } The rigid-body positions and orientations, r and θ , in the body axis position vector, { ˆ η } , b/e,b b/e,b b represent the position and orientation of the vehicle, which has no effect on the aerodynamic forces. The remaining states in the { ˆ η } vector, which are the structural states, are inertial. Therefore, derivatives of the b ∗ d structural states are equal to the rates, ξ ≡ ξ . Therefore, the GAF at steady-state ( k = 0 ) can be handled ˆ d t as a separate case.

1. Polynomial Approximation For low frequencies, it is possible to fit a simple polynomial model to match the GAFs. To generate a meaningful frequency response, the real part must be an even function (symmetric when mirrored across k = 0 ) and the imaginary part must be an odd function (equal when rotated 180 deg about k = 0 ). These symmetry constraints result in the polynomial model with real coefficients shown in Eq. (12).

n [ ] [ ] [ ] ∑ ∗ j − 1 ˜ ˜ ˆ ˆ Q ( ik ) − Q (0) T ≈ Q ( k ) = ( ik ) Q (12) 0 j j =1 A polynomial model would require an excessive level of complexity to accurately capture the GAFs at moderate to high frequencies. A polynomial model is, however, an effective way of generating a model of the quasi-steady aerodynamics, Eq. (13), at very low frequencies ( k > 0 . 05 ).

[ ] [ ] [ ] d ∗ d ˆ ˆ Q ( k ) ≈ Q ( k ) = i Q (13) dk k =0 dk k =0 The quasi-steady model is helpful in fitting the more sophisticated approximation shown in section III.C.3.

The quasi-steady model is fit by the method of collocation[ 11 ]. Collocation is a method of finding the solution which is exactly equal at a fixed set of points. The polynomial model of order 2 ( n − 1) , is fit to n frequencies of the body axis GAF, Eq. (12). Because of the singularity at zero frequency in Eq. (A-23b) in the appendix, the zero frequency case is excluded from the collocation. By separating the GAF into real and imaginary parts, the solution for the polynomial coefficients are uniquely determined.

The symmetry constraints mean that the even and odd portions of the polynomial model can be solved separately, as shown in Eq. (14).

      ∗ 2( n − 1) ˆ i Im Q ( k ) ik I − ik I · · · ( ik ) I Q 2 2 2 1 ∗       2( n − 1) ˆ  i Im Q ( k )   ik I − ik I · · · ( ik ) I   Q  3 3 3 3       ≈ (14) . . . . .

.

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

.

 .   . . .   .  ∗ 2( n − 1) ˆ i Im Q ( k ) ik I − ik I · · · ( ik ) I Q n n n 2( n − 1) n [ ] ˆ To determine the Q parameter in Eq. (13), only the odd terms of Eq. (13) are needed. Therefore, the even terms are not discussed further. Because a collocation method is used, numerical complications that produce small errors in the GAF can cause significant variations in the coefficients. There are singularities in the Kernel function used by DLM or ZAERO ™ at zero reduced frequency that can cause numerical errors at very low frequencies. On the other hand, the generation of the quasi-steady coefficients, Eq. (13), entails using very low frequencies. Therefore, care must be taken to ensure that numerical errors are not caused by specifying the frequencies too close to zero.

2. Mean Axis Rational Function Approximation The difficulty of frequency domain methods is that, even for the simplest two-dimensional 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 a rational function approximation (RFA) to describe these unsteady aerodynamics. The rational function approximation is a transfer function. For flutter suppression experiments on a B-52 airplane (The Boeing Company, Chicago, Illinois), Roger[ 6 , 7 ] developed the idea of an RFA, shown in Eq. (15), that could be used for controller design.

([ ] [ ] [ ] [ ] ) ( [ ] [ ]) − 1 2 r r ˆ ˆ ˆ ˆ { ˜ q ( k ) } ≈ − A + A ik + A k + [ D ] ik I − R E ik { η } (15) 0 1 2 r r h [ ] r ˆ There are several different methods of defining the form of the [ D ] , [ R ] , and E matrices. For example, r Karpel’s minimum state[ 13 ] will typically give a low order model of the unsteady aerodynamics. The lower order model is appealing for control law development; however, the fitting of model is not convex, and does not have a unique solution. The lack of convexity makes it difficult or impossible to ensure a consistent definition of these aerodynamic states. In Karpel’s minimum state method, the aerodynamic lags are not associated with a specific degree of freedom, which further complicates the consistent definition of the states.

The Roger’s RFA is the original method[ 7 ], and remains one of the most popular. The traditional Roger’s RFA specified a number of poles that are repeated for each degree of freedom in the structural model. The r matrix [ D ] then has a specified structure, as shown in Eq. (16). The aerodynamic lag poles in [ R ] must r be defined by the engineer developing the model. A conventional method of distributing the poles over the same frequencies of the available GAF matrices has been used to select these aerodynamic lags.

  1 1 · · · 1 0 0 0 0 0 · · · 0    0 0 · · · 0 1 1 · · · 1 0 0 · · · 0    [ D ] := (16) . . . . . . . .

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

 .  . . . . . . . .

0 0 · · · 0 0 0 · · · 0 1 1 · · · 1 The challenge within all of the prior traditional RFA forms in Eq. (15) 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 integrate the RFA into a model of the form of Eq.

equation (1d) , 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.

The non-dimensional body axis states and positions are used to define a new body axis RFA as shown in Eq. (17).

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

3. Fitting Time Domain Rational Function Approximation There are many existing codes for the traditional RFA in Eq. (15), 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[ 14 ].

Using this approach, the time domain transformation, from Eq. (A-19a) in the appendix, is applied to the traditional RFA, Eq. (15), to give the matrices (Eqs. (18a), (18b), (18c), and (18d)) for the body axis RFA, Eq. (17).

[ ] [ ] [ ] [ ] [ ] ∗ ˆ ˆ ˆ ˆ ˆ A = A T + A T ∗ (18a) 0 he 1 he [ ] [ ] [ ] [ ] [ ] [ ] [ ] − 1 ∗ ˆ ˆ ˆ ˆ ˆ ˆ ˆ A = A T ∗ + A T ∗ T T ∗ (18b) 1 1 2 he hb he he [ ] [ ] [ ] ∗ ˆ ˆ ˆ A = A T ∗ (18c) hb [ ] [ ] [ ] ∗ ˆ ˆ ˆ ∗ E = E T (18d) hb [ ] [ ] ˆ ˆ The use of the modal RFA results in erroneous coefficients, A T ∗ , that have no physical basis. The he 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.

The aerodynamic forces should be independent of the inertial rigid-body positions and orientations.

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, see Eq. (18a). Attempting to correct these errors after fitting and applying the transformation will increase the total error in the approximation and will specifically add a bias to the approximation.

Instead, the fitting is applied directly to the body axis GAFs shown in Eq. (11). To fit the RFA, first several constraints are defined. These constraints give consistency between the frequency domain GAFs and the RFA. It is a common practice to constrain the RFA to exactly equal the GAF at k = 0. The application of a separate transformation at k = 0 makes the matching a necessity. The constraint is achieved by selecting [ ] ∗ ˆ the A matrix shown in Eq. (19).

[ ] ∗ [ ] ˆ ˜ A = − [ Q (0)] T (19) Applying the steady-state constraint to the RFA, Eq. (17), removes the Earth-fixed axis displacements from the RFA, resulting in Eq. (20).

[ ] ([ ] [ ] [ ] ) ∗ ∗ ∗ ( [ ] [ ]) ∗ − 1 2 r r ˜ ˜ ˆ ˆ ˆ Q ( ik ) − Q (0) T = − A ik − A k + [ D ] ik I − R E ik (20) 0 r r 1 2 ∗ ˆ In addition to the constraint on the steady flow matrix ( A ), a constraint on the quasi-steady coefficients of the RFA is defined to equal the quasi-steady polynomial model in Eq. (13). As was done with the polynomial model, the quasi-steady RFA solution is calculated from the limit of the RFA in Eq. (20). The quasi-steady RFA solution is shown in Eq. (21).

[ ] [ ] ∗ ∗ d ˜ ˜ ˆ lim Q ( ik ) − Q (0) T = − i A (21) k → 0 dk By equating the quasi-steady models of Eqs. (13) and (21), the quasi-steady constraint is defined in Eq. (22).

[ ] [ ] ∗ ˆ ˆ A := − Q (22) The final constraint comes from the separation of the deflections and velocities in the RFA of Eq. (17). For ∗ ˆ the non-inertial, mean axis rigid-body, states, the corresponding terms of A are zero, because the rigid-body position does not change the aerodynamic loads. To ensure that the steady-state rigid-body aerodynamic r r forces remain an algebraic function, the aerodynamic lag portion of the RFA, [ D ] ( ik [ I ] − [ R ]) [ E ] , r r must go to zero at steady-state. The same constraint is applied to the non-inertial states, the modal DoF, ∗ ˆ because it is convenient that the A matrix remains the complete quasi-steady coefficients. To satisfy these steady-state specifications, the final value theorem[ 15 ] applied to the RFA in Eq. (20) gives the constraint in Eq. (23).

( ) [ ] [ ] ∗ ∗ ˜ ˜ ˆ lim Q ( ik ) − Q (0) T ≡ − A (23) k → 0 ik Now that the constraints are defined, the error function to minimize in the RFA fitting is defined. To ensure that the coefficients of the RFA are real numbers, the RFA, Eq. (20), is separated into its real part, Eq. (24), and imaginary part, Eq. (25).

[ ] [ ] ( ) [ ] ∗ ∗ − 1 [ ] [ ] ∗ 2 2 r r 2 ˜ ˜ ˆ ˆ Re Q ( ik ) − Q (0) T ≈ A k − [ D ] k I + R E k (24) 0 r r [ ] [ ] ( ) [ ] ∗ ∗ − 1 [ ] [ ] [ ] ∗ 2 r r r ˜ ˆ ˆ Im Q ( ik ) ≈ − A k + [ D ] k I + R R E k (25) r r r These equations can be rewritten to define the RFA error in a matrix form as shown in Eq. (26). 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   ˆ := ( ) [  ] + ( ) (26) RF A  A  − 1 1 ∗ 2 2 ˜ I 0 − D k I + R R ∗ Im Q ( ik ) ˆ E k The parameters of the RFA are then defined by the constrained linear least squares in as shown Eq. (27).

[ ] T ∗ ∗ ∗ T ˆ ˆ ˆ = argmin [  ] [  ] A A E RF A RF A 1 2 [ ] [ ] ∗ ˆ ˆ (27) subject to A = − Q [ ] [ ] − 1 ∗ r ˆ [ D ] R E = [ 0 ] r Numerous commercial or open libraries exist to solve linear least squares with linear constraints.

IV. Calibration of the Model Forces Linear lifting surface methods such as DLM or ZAERO ™ are relatively low fidelity methods. The accuracy limitations become apparent at low frequencies and when using control surfaces. The effectiveness of control surfaces can be significantly affected by viscous effects[ 16 ] which are not modeled in potential flow.

Furthermore, these methods cannot accurately model a thick lifting body, which can be significant for a flying wing. Despite these limitations, linear methods remain the primary tool for flutter analysis because they are effective for many cases. The process of fitting the RFA in section III.C.3 ensures a good match with the GAF, but the GAF may not accurately represent the aerodynamic forces of the vehicle. Instead of using a complete high fidelity model, the low fidelity linear models are refined based on limited high fidelity data.

For example, steady CFD can be used to adjust the the unsteady linear model. The corrections are done in two stages. First, the high dimensional AIC matrices are updated directly. The second stage of updating the forces is to add a correction to the RFA to capture any missing effects.

A. Calibration of Aerodynamic Influence Coefficients Correcting the AIC matrices of Eq. (9) offers a method of extrapolating limited information to modes that may not have been examined. For example, rigid-body loads due to angle of attack can be used to improve the bending loads. Extrapolating data is problematic and attempting to overcorrect can decrease the accuracy of the models. Although there have been many efforts to apply such corrections, the extrapolating nature of the methods means that they are often problematic.

Giesing et al.[ 17 ] proposed a correction to the downwash, ε , at each panel in Eq. (28). The downwash correction is frequently called a postmultiplying correction factor.

[ ] ([ ] [ ]) [ ] r r F ( k ) = I + ε F ( k ) (28) r r Britt et al.[ 18 ] observed that the unsteady methods were quite accurate at higher frequencies, but required more corrections at lower frequencies. For the Northrop Grumman (Falls Church, Virginia) B-2 airplane, a heuristic based approach of linearly scaling the correction factors to zero at a reduced frequency of 0 . 4 was used.

To better understand the limitations of the downwash corrections and to improve the frequency dependence, the irrational correction methods of Ashley[ 19 ] are examined. Ashley’s correction for thickness and boundary-layer effects offer a more physical solution.

n j u V L j,j α j
x j
t Panel j Fig. 6 Panel downwash.

For a panel as shown in Fig. 6, Ashley[ 19 ] examined the definition of the general downwash, as shown in Eq. (29).

T := f { n } { u } (29) j j j,j V ∞ The velocity in the Eulerian perspective, u in Eq. (30), is given by the material derivative[20].

j,j D ∂ { u ( x , t ) } = { x } := { x ( x , t ) } − [ ∇ x ( x , t )] { V ( x , t ) } (30) j,j j j x j L Dt ∂t ∂ The first term of Eq. (30), { x ( x , t ) } , represents the unsteady motion of the boundary. The motion of the j ∂t boundary is well known and easily expressed from the mode shapes. Therefore, no correction of this term is considered. The second term of Eq. (30), [ ∇ x ( x , t )] { V ( x , t ) } , represents the local velocity projected onto x j L the frame of the panel. The primary interest for the correction factors is this local velocity term. Applying the local velocity of Eq. (30) to the downwash in Eq. (29), the steady downwash can be reparameterized as a local angle of attack, α , and the free stream velocity, V , as shown in Eq. (31).

j ∞ 1 ‖ V ( x , t ) ‖ T L { n } [ ∇ x ] { V } = α ( x , t ) (31) j x j L j V V ∞ ∞ Linearized potential flow solvers assume small disturbance in the flow velocity. The small disturbance assumption means that the local velocity is equal to the freestream velocity, ‖ V ( x , t ) ‖ ≈ V . Under this L ∞ [ ] [ ] ¯ ˜ assumption, the downwash can be separated into steady, F , and unsteady parts, F , in Eq. (32). Both of these matrices are independent of reduced frequency.

[ ] [ ] [ ] ¯ ˜ F ( k ) = F + ik F (32) [ ] ¯ In Eq. (33), the terms of the steady downwash matrix, F , are shown to correspond to the steady downwash from Eq. (31).

( ) [ ] [ ] [ ] ‖ V ( x , t ) ‖ T L ¯ F U Φ { η (0) } ≡ α ( x , t ) (33) 0 j h j V ∞ Because the actual airfoil has a non-zero thickness, the local velocity will differ from the freestream. The assumption of the local velocity is relaxed from equal to the freestream, to a steady velocity. The ratio of local to freestream velocity then becomes a constant correction to the downwash.

[ ] ¯ The steady downwash matrix, F , is defined by a mean surface that is defined by the outer mold line of the vehicle. The panels are also assumed to oscillate about a mean flow condition aligned with the x -axis.

The local angle of attack, α , is therefore a function only of the panel orientation. By relaxing the assumption j from no boundary layer to a steady boundary layer, the boundary layer can be modeled by altering the panel orientation. The change in orientation is well represented by a multiplicative correction applied to the local angle of attack, α .

j The boundary layer and thickness corrections are combined into a single multiplicative correction factor in Eq. (34). The correction factor is separated into an known, ¯ ε , and an unknown correction factor, ε .

Separating into known and unknown parts allows for the integration of different sources of data and to incrementally improve the AIC from higher fidelity models.

1 ‖ V ( x , t ) ‖ T L { n } [ ∇ x ] { V ( x , t ) } = (¯ ε + ε ) α ( x , t ) (34) j j L j j j V V ∞ ∞ These correction factors are integrated into the GAF by generating a corrected downwash matrix as shown in Eq. (35).

[ ] ([ ] [ ]) [ ] [ ] ∗ r r ¯ ˜ F ( k ) = ¯ ε + ε F + ik F (35) r r At steady-state ( k → 0 ), this is equivalent to the post-multiplying correction factor of Giesing et al.[ 17 ] in Eq.

(28). As the reduced frequency increases, the steady portion of the downwash becomes smaller relative to the unsteady portion. Therefore, the effect of the correction decreases as the frequency increases. The reduced effect of the correction factors gives a similar effect to the heuristic of Britt et al.[18].

From Ashley’s expansion of the downwash, one can expect that the downwash corrections would be quite effective at correcting the effects of moderate thickness and to correct for the effect of the boundary layer changing the orientation of the mean flow between the upper and lower surfaces. The method cannot correct the model to match viscous forces. Attempting to correct the model to match coefficients dominated by viscous forces will cause erroneous coefficients, and could make the model less accurate.

The use of diagonal downwash corrections means that the AIC has been corrected for a single input. In the present work, only data with respect to angle of attack are used. Moreno et al.[ 21 ] shows a method of extending this to diagonally dominant corrections, that are able to correct for additional modes. Correcting for additional inputs would require specialized data beyond what may typically be available. Therefore, the enhanced correction factor matrices were not explored.

Giesing et al.[ 17 ] proposed the corrections as the smallest magnitude of smooth corrections that adjust the AIC matrix to exactly match a select set of steady aerodynamic coefficients, C . There are many more match panels in the aerodynamic model then the number of coefficients being matched. Therefore, it is possible to match the coefficients exactly. The corrections are specified as the smallest smooth change to keep the corrections reasonable.

The first step in generating the correction factors is the determination of which coefficients to match. The determination of coefficients is heavily dependent on the vehicle and the available data. The conclusions of Giesing et al.[ 17 ] provide a very good discussion of the process of selecting coefficients, and will not be replicated here.

Once the coefficients to match are determined, a column of the incidence matrix is defined to aid in selecting the coefficients of interest from the full GAF, [ Q (0)] . The incidence matrix is defined for each input/output pair in Eqs. (36a) and (36b).

  1 derivative number a due to mode h input I = (36a) ha  0 otherwise   derivative number a acts on mode h output SL I = (36b) ha  0 otherwise The output incidence matrix, Eq. (36b), is normalized by the vehicle reference area, S , and the characteristic ¯ c length, L = . The incidence matrices are used to define the constraint for the minimization in Eq. (37) used by Giesing et al.

T { ε } = argmin { ε } [ W ] { ε } ( [ ] T [ ] [ ] [ ] [ ] T output subject to diag I Φ U S AIC (0) ha (37) [ ]) ([ ] [ ]) [ ] [ ] [ ] { } T input r r ¯ ¯ ε + ε F U Φ I = C r r 0 match ha Past experience[ 17 ] has shown it is necessary to define separate correction factors for the modes and the control surface inputs. Each set of corrections factors, ε , uses a separate definition of weighting matrix, [ W ] .

The weighting for corrections in Eq. (38a) from rigid and flexible modes, [ W ] , is the lift due to a unit mode downwash. For control surface modes, the weighting, [ W ] , shown in Eq. (38b) is determined by the surface { } ¯ force at the mean flow condition. The mean flow pressures, C , are taken from the CFD results. The p weighting matrices used are similar to Giesing et al.[ 17 ], but they have been modified to include the panel area. By adding the area, the weighting is for a given force, rather than pressure. Using the force for the weighting better handles irregular panel sizes.

∣ [ ] [ ] ∣ r ∣ ∣ [ W ] = S AIC ( k ) { 1 } (38a) mode r ∣ [ ] { }∣ r ∣ ¯ ∣ [ W ] = S C (38b) surface r p Giesing et al.[ 17 ] proposed finding smooth coefficients; however, the implementation of Eq. (37) does not actually smooth the corrections. It minimizes the total variations, but does not limit jumps between neighboring panels. For subsonic flow, one would not generally expect jumps in the local flow velocity or the boundary layer, implying that the corrections should be smooth as well. The smoothness is enforced by adding a penalty to the cost function of Eq. (37) to allow a larger total correction if it will be smoother. Adding the smoothness is effectively a multiobjective optimization by the  -constraint method. The multiobjective optimization is simultaneously minimizing the variance, as in Giesing et al, and a new smoothness parameter, λ . The smoothing parameter, λ , is selected primarily by trial and error. The parameter is increased until hp hp the correction factors are sufficiently smooth.

The formulation for the penalty for the smoothness of the corrections is defined by using a modification of the Hodrick–Prescott filter[ 22 ]. The Hodrick–Prescott filter was intended for smoothing one-dimensional data, such as time series. The smoothness is achieved by limiting the mean squared difference between adjacent terms. Tibshirani and Taylor[ 23 ] used a similar 1 -norm based smoothing and extended it to an arbitrary two-dimensional geometry by using graph theory. Using a map of the United States, the change in the parameter between adjacent states is constrained. By applying the idea of the adjacency of a graph, the Hodrick–Prescott filter can be extended to improve the smoothness of the correction factors.

The adjacency penalty of Tibshirani and Taylor[ 23 ] is implemented by construction of an incidence matrix, I , in Eq. (39). The incidence matrix describes which nodes of the aerodynamic panels lie on the edge of a kj panel, as shown in Fig. 7.

Panel j d kj Node k Fig. 7 Construction of aerodynamic panel-node incidence matrix.

The typical incidence matrix would specify a 1 when a node is on the edge of a panel. Because the panels can be of substantially different size, incidence matrix is weighted by the distance from the node to the control point.

  node k lies on edge of panel j d kj I = (39) kj  0 otherwise The weighted incidence matrix of Eq. (39), is used to augment the optimization in Eq. (37) to penalize for the smoothness of the correction factors. The augmented optimization is shown in Eq. (40).

T T T { ε } = argmin { ε } [ W ] { ε } + λ { ε } [ I ] [ I ] { ε } hp kj kj ( [ ] T [ ] [ ] [ ] [ ] T output subject to diag I Φ U S AIC (0) (40) ha [ ]) ([ ] [ ]) [ ] [ ] [ ] { } T input r r ¯ ¯ ε + ε F U Φ I = C r r 0 match ha The formulation of Eq. (40) can be equivalently expressed as a constraint on the smoothness. In the case where λ = 0 , the Eq. (40) is identical to the method of Giesing et al.[17] shown in Eq. (37).

hp The use of graph theory is a fairly rudimentary way of describing the smoothness of the parameters. Some panels may have more neighbors than others, which can cause large jumps in the corrections if the smoothing parameter is too large. For smaller values of λ , however, it works very well.

hp B. Rational Function Approximation Augmentation The second method of correcting the aerodynamics is an additive correction to the RFA matrices after the fitting of Eq. (27). The additive correction to the RFA allows for the inclusion of additional physical effects that cannot be effectively captured in the AIC, such as viscous forces. These additional physics can become critical for several of the flight dynamics modes.

The lifting surface methods being used for the aerodynamics models do not provide a value for the drag, or a change in the aerodynamic forces due to velocity changes. 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 matrices are added to the original RFA matrices, as shown in Eqs. (41a), (41b), (41c), and (41d).

( ) [ ] [ ] [ ] ∗ ∗ ∗ r ˆ ˆ [ A ] ≡ A + A D (41a) x r 0 0 0 grav ( ) [ ] [ ] [ ] ∗ ∗ ∗ r ˆ ˆ [ A ] ≡ A + A D (41b) u r 1 1 1 aug [ ] [ ] ¯ c ∗ ∗ r ˆ [ A ] ≡ A D (41c) u r 2 2 2 V e [ ] [ ] 2 V e ∗ ∗ r ˆ [ E ] ≡ E D (41d) 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. (42), is captured by the change in velocity.

∂ ¯ q ¯ q = 2 (42) ∂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. (17). An additive correction to the [ A ] matrix, Eq. (43), 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[ 4 ]. The only difference is the addition of the forces on the structural 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      − 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  (43)   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, the gravitational force on each element of the FE model is the linear function of the displacements, shown in Eqs. (44) and (45)[4].

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

.

 . . .  0 0 · · · T grav   0 0 0 0 1 0    0 0 0 − cos θ 0 0      0 0 0 0 0 0   [ T ] = (45) 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 body 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 FE model.

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 non-dimensionalization of the gravity in Eq. (46).

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

 .  . . . .

grav ¯ q .

 . . .  0 0 · · · T grav V. Integration of the Models The components of the models defined in section III are integrated by the dynamic and kinematic relations of the equations of motion. The equations of motion in the modal frame with external forcing are identical to the classical flutter equations of motion, Eq. (47), from Bisplinghoff[8].

/ 2 [ M ] { ¨ η } + [ B ] [ K ] { ˙ η } + [ K ] { η } = { q ( t ) } (47) h h h The transformation from modal to body axis, Eq. (A-19a) in the appendix, is used to generate equations of motion in the mean body axis. From Ouellette[ 24 ], the terms of Eq. (1d) can be replaced by the matrices from the previous sections to give Eq. (48). The control surface modes, δ , are separated in the equations of motion to represent the inputs to the system.

      2 V − 1 2 V − 1 0 − 1 e − 1 ˆ ˆ ˆ ˆ   D T T ∗ D D T T ∗ D 0   η x u η   x he x he   s ¯ c ¯ c b he hb d  ( ) ( )  − 1 ∗ − 1 ∗ − 1 ˙ = ¯ ¯ ¯ ¯ ¯ ˙ η   η − M K + ¯ q A − M B + ¯ q A − ¯ q M D s b 0 1     dt     2 V ∗ r 0 ¯ x x aero 0 E − R aero r ¯ c       0 0 0 δ   ( )   − 1 ∗ − 1 ∗ − 1 ∗ ¯ ¯ ¯ ˙ + (48)  − ¯ q M A − ¯ q M A − M M + ¯ q A  δ hc hc hc hc 0 1 2     ∗ ¨ ˜ 0 E 0 δ hc Equivalent mass (Eq. (49a)), damping (Eq. (49b)), and stiffness (Eq. (49c)) matrices are defined to simplify the state space of Eq. (48).

[ ] [ ] [ ] [ ] [ ] 2 V r ∗ ¯ ˆ M ≡ M T ∗ D + ¯ q A (49a) u r ¯ c hb ( ) [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] / 2 − 1 2 V 2 V 0 0 r ¯ ˆ ˆ ˆ B ≡ M B K + T ∗ T T ∗ D (49b) he u r ¯ c he hb ¯ c [ ] [ ] [ ] [ ] r ¯ ˆ K ≡ K T D (49c) he x r VI. Results The X-56A flight test data were selected as the primary example for study due to the availability of flight test data at speeds faster than flutter. The primary flutter mode of the X-56A MUTT is a body freedom flutter mode, the mechanism of which is the short-period being driven unstable by the structural dynamics. The use of a flutter suppressing control law has allowed flight operations to airspeeds past where the open-loop flutter mode would go unstable. The FE model was implemented in NASTRAN ™ (MSC.Software Corporation, Santa Ana, California). The ZAERO ™ , subsonic lifting surface method was selected to provide the frequency domain GAF, [ Q ( k )] in Eq. (8), for the current modeling.

A. Correction Factors Examination of the correction factors yields some interesting insights to the limitations of the underlying modeling data, and the validity of the corrections themselves.

For the X-56A MUTT, the AIC corrections were done in two stages. First, a set of corrections was generated from CFD data. The CFD data provided the change in lift, pitch moment, bending generalized force, and torsional generalized force with the angle of attack. After generating a set of correction factors ¯ from the CFD, the CFD based factors were then used as the bias term, ε in Eq. (35), in the generation of corrections from wind-tunnel tables. Because the solution of the correction factors, Eq. (40), is giving the smallest change from the bias, ¯ ε , the the CFD correction factors will influence the wind-tunnel correction factors. The use of CFD data provides more information on the distributed pressures. The more distributed information does a better job getting an accurate distribution of the correction factors by matching additional output. The wind-tunnel data are generally more accurate, and are used as a fine tuning of the correction factors. The primary exception to the accuracy of the wind-tunnel tables was the value of C . The initial m α flight tests showed the wind-tunnel data was underestimating C . The wind-tunnel data were then tuned m α to match the CFD results, which were more consistent with flight.

The CFD data were only available for angle of attack variations as inputs. Because the data were mapped to the AIC aerodynamic boxes, it was possible to generate the coefficients for more outputs as shown in Table 3.

Table 3 Coefficients used for AIC tuning.

CFD coefficients Wind-tunnel coefficients Lift curve slope Lift curve slope C C L L α α C Pitch stiffness C Pitch stiffness m m α α Wing bending due to angle of Hinge moment due to angle of C C SW 1 B hm α α attack attack (except body flap) Wing twist due to angle of attack C SW 1 T α Hinge moment due to angle of C hm α attack The resulting distribution of correction factors across the aerodynamic panels is shown in Fig. 8a.

1 1 0.5 0.5 0 0 −0.5 −0.5 −1 −1 .

(a) CFD based. (b) Wind-tunnel based.

Fig. 8 Distribution of AIC correction factors, ε .

The largest correction factors are associated with the forward part of the centerbody. These centerbody correction factors are associated with the thickness not modeled in ZAERO ™ . The inclusion of hinge moments improves the corrections associated with the boundary layer at the trailing edge.

The wind-tunnel coefficients are limited to the rigid-body modes; however, additional inputs are available.

The selected wind-tunnel coefficients for matching are shown in Table 3. The resulting wind-tunnel based corrections shown in Fig. 8b appear almost identical to the CFD based corrections. The magnitude of the adjustment for the centerbody thickness is slightly smaller, but very similar. Also the corrections along the hinge line of the control surfaces is slightly higher.

B. Model Comparison Two different frequency domain models are generated from the flight data for comparison to the models. The first model is a low order equivalent system (LOES) transfer function, Eq. (50). The order of the LOES models was selected to capture the short-period, symmetric wing bending, and symmetric wing torsion modes.

∑ i n s i i =0 H = (50) loes ∏ 2 2 ( s + 2 ζ ω s + ω ) i i i =1 i The parameters of the transfer function are determined by a time domain output error system identification method[ 25 ]. The transfer function of the open-loop plant, without a control system, is estimated directly from the closed-loop flight data. All 10 surfaces are excited in symmetric pairs by an orthogonal multisine input. There is some error in the transfer functions because they are fit at a slightly different flight condition.

The main source of error is that the model is of insufficient order to capture all of the dynamics.

The second flight derived model is a high order nonparametric frequency response. The frequency response is estimated from a single surface excitation, and is corrected to get the open-loop response. The frequency response are the higher fidelity of the two models, but are prone to significant variability.

1. Frequency and Damping Trends The frequency and damping trends of the short-period mode with respect to airspeed are shown in Fig. 9.

The variation of the frequency and damping with respect to fuel weight for the predictive models are shown by the contours. The frequency and damping from the fit of the model of Eq. (50) to the flight data are shown by the discrete diamonds.

Fuel, % Fuel, % 1 −1 1 0.8 0.8 −0.5 0.8 0.6 0.6 0.6 Damping 0.4 0.4 0.4 Relative Frequency 0.5 0.2 0.2 0.2 Model Model Flight Data Flight Data 0 1 0 0 Airspeed Airspeed . .

(a) Frequency relative to flutter frequency. (b) Damping.

Fig. 9 Frequency and damping of the short-period mode.

Figure 9b shows that the short-period mode is going unstable and becoming the body-freedom flutter mode. The damping shows a weak relationship to the fuel weight, increasing slightly as fuel weight increases.

Although the fuel weight has a small effect of the damping, the flutter speed does change significantly with fuel weight. The models do show a reasonably accurate prediction of the damping when compared to the flight data. The flight data does appear to show the short-period becoming unstable more quickly with airspeed than predicted by the models. The weak relationship with fuel weight makes it difficult to make conclusive observations on the accuracy of the fuel weight effects on damping in the model. The frequency of the short-period in Fig. 9a shows a strong dependence on the fuel weight. The frequency of the short-period mode is comparing well with the flight test analysis.

The frequency and damping trends of the symmetric wing torsion mode with respect to airspeed are shown in Fig. 10.

Fuel, % Fuel, % 1 1 −1 1 0.8 0.8 0.8 −0.5 0.6 0.6 0.6 0.4 Damping 0.4 0.4 Relative Frequency 0.5 0.2 0.2 0.2 Model Model Flight Data Flight Data 0 1 0 0 Airspeed Airspeed . .

(a) Frequency relative to in vacuo frequency. (b) Damping.

Fig. 10 Frequency and damping of the first symmetric torsion mode.

The frequency and damping of the torsion mode show almost no dependence on the fuel weight. The damping of the torsion mode in Fig. 10b is nearly constant with airspeed. It appears that the structural damping term is strongly dominant, so the model damping shows very little change with airspeed, unlike the other modes. The frequency of the torsion mode in Fig. 10a shows a slight decrease with airspeed. The model appears to capture the trend of the frequency well, and the error appears to be a bias. Because the model captures the trend of the frequency well, it is likely an error of the in vacuo frequency from the FE model.

The frequency and damping of the bending mode shown in Fig. 11 do show more significant differences from the model predictions.

Fuel, % Fuel, % 1 −1 1 0.8 0.8 0.8 −0.5 0.6 0.6 0.6 Damping 0.4 0.4 0.4 Relative Frequency 0.5 0.2 0.2 0.2 Model Model Flight Data Flight Data 0 1 0 0 Airspeed Airspeed . .

(a) Frequency relative to in vacuo frequency. (b) Damping.

Fig. 11 Frequency and damping of the first symmetric bending mode.

The frequency in Fig. 11a shows a moderate influence of fuel weight; however, the variability of the flight test results is larger than the influence of the fuel weight shown by the models. Therefore, it is difficult to make conclusions as to the accuracy of the models in capturing the fuel effects. There is, however, a clear error in the change in frequency with airspeed. The model shows the correct trend, but the flight test data show a much larger change in frequency. The damping in Fig. 11b matches fairly well at low speeds. At higher airspeeds, the damping from the flight data departs from the model and trends toward 1 . The models do not capture the change in the damping trend.

The only CFD used for the the calibration of the AIC matrices in section VI.A was the angle of attack.

The corrections were sufficient for giving a good representation of the short-period mode, which is dominated by changes in angle of attack. The corrections cannot, however, be expected to significantly improve the model with respect to the wing bending mode. As the bending mode effects become dominant, the limitations of the AIC correction become apparent.

2. Frequency Responses The frequency response for the pitch rate and the vertical accelerometer at the wing tip are shown in Fig. 12 for a low speed case. The input is symmetric excitation of the outboard-most flaps. These input and output pairs give a good excitation of the longitudinal dynamics.

Bending Short-Period Bending Short-Period Torsion Torsion Flight Flight Model Model Magnitude, dB Magnitude, dB LOES LOES Phase, deg Phase, deg 1 1 Coherence Coherence 0 0 Frequency, Hz Frequency, Hz (a) Pitch rate. (b) Left outboard accelerometer.

Fig. 12 Symmetric outboard surface excitation at low speed.

At low speeds the predictive model matches exceptionally well to the flight test data. The LOES model does show a difference in the phase of the accelerometer. The variations in phase are primarily a result of the signal processing used before fitting the model. The nonparametric flight data show an unrealistically high coherence. This is partially due to a high signal to noise, but also due to bias in estimating the coherence.

There were insufficient repetitions of the maneuver to generate an accurate estimate of the coherence, so it is heavily biased toward the maximum. The predictive model does show a slightly higher magnitude at the short-period mode.

The frequency response for the pitch rate and the vertical accelerometer at the wing tip are shown in Fig.

13 for a case very close to the flutter speed. At the higher flight speeds, the limitations of the predictive model is more apparent.

Body Freedom Body Freedom Flutter Torsion Flutter Flight Flight Torsion Model Model Magnitude, dB Magnitude, dB LOES LOES Phase, deg Phase, deg 1 1 Coherence Coherence 0 0 Frequency, Hz Frequency, Hz (a) Pitch rate. (b) Left outboard accelerometer.

Fig. 13 Symmetric outboard surface excitation at the flutter speed.

At these higher speeds, the damped frequency of the bending mode and the short-period mode are very close. Because the bending mode is heavily damped, it is nearly impossible to distinguish the bending mode from the frequency response. The prominent peak is the body freedom flutter (short-period) mode. The predictive models are showing a lower frequency for the peak. It is also worth noting the variability in the prediction of the peak from the flight data based models. The variation appears to come from the differences in the methods and assumptions in determining the models from the flight data, rather than the frequency resolution of the flight data. The frequency resolution of the maneuvers being used to generate these estimates is much smaller than the variation in the frequency of the peak. The coherence in 13 remains high, but does appear to be more realistic. Although errors in the predictive models are showing in the frequency response at flutter, they remain a fairly useful model for designing and evaluating a control system.

VII. 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 wind-axis system. Methods of correcting these low fidelity models based on limited high fidelity data were discussed. An assumed modes method was selected for describing 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 wind across the flight envelope. The approach results in models that are very useful for flight control design.

These methodologies were applied to the X-56A Multi-Utility Technology Testbed (MUTT) and were compared against flight-test data for the flex wing configuration for airspeeds both below and above the flutter speed. The models generated by the present approach are able to accurately capture the coupled flight dynamics and structural dynamics of the X-56A MUTT.

Appendix The details of the non-dimensionalization and the complete derivation of the coordinate transforms are provided here in the appendix. The non-dimensionalization is expanded in a matrix format for clarity with the previous equations. The state transformation matrices are derived from the fundamental kinematic equations and translated into a form more usable for the current linear models.

A. Non-dimensionalization Because the unsteady aerodynamics are a function of time, generation of non-dimensional coefficients requires a non-dimensional definition of time. The non-dimensionalization of time, Eq. (A-1), is determined by the true airspeed at trim and the chord length.

2 V ˆ t = t (A-1) ¯ c { } ∗ The non-dimensional position, { ˆ η } in Eq. (A-2a), and rates, ˆ η in Eq. (A-2b), are related to the b b r r dimensional states by the diagonal matrices, [ D ] in Eq. (A-3a) and [ D ] in Eq. (A-3b).

x r u r [ ] r { ˆ η } = D { η } (A-2a) x r b b { } ∗ [ ] r ˆ η = D { ˙ η } (A-2b) u r b b The non-dimensionalizing matrices of Eqs. (A-3a) and (A-3b) 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 = (A-3a) x r 0 0 0 1 0 0 0     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 =   (A-3b) 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. Transformation Between Flight Dynamic (Mean) and Modal Axis Systems The structural and unsteady aerodynamic model are defined in the modal coordinate system. As discussed in section III, an inertial coordinate system, such as a modal coordinate system, is insufficient for describing the vehicle flight dynamics. The non-inertial body axis system is a more traditional form used in flight dynamics modeling and simulation, which improves the familiarity of the models for a controls engineer. 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, i.e., mass and stiffness, and for kinematic relationships in the state space models. For the aerodynamic model generation, a frequency domain form of the transformation is derived from the time domain form.

The goal of the time domain transformation is to generate kinematic relationships describing the motions in the modal frame, Eqs. (A-4), (A-5), and (A-6), as a function of the known states in Eqs (4a) and (4b).

{ } T T T T { η } = r θ ξ (A-4) b b/h,h b/h,h { } T T T T ˙ ˙ { ˙ η } = (A-5) ˙ r θ ξ b b/h,h b/h,h { } T T T T ¨ ¨ ¨ { η } = ¨ r θ ξ (A-6) b b/h,h b/h,h First the full nonlinear transformations of the vectorized form in Eq. (A-7) are calculated.

        ¨ η T (¨ η , ˙ η )    accel  h b b ˙ η = T ( ˙ η , η ) (A-7) rate h b b         η T ( η ) disp h b Then the transformations are linearized for implementation of the state space model.

The transformation functions of Eq. (A-7) are calculated for each term of the state matrix. Because the modal and Earth frames are both inertial, it is easiest to begin in these frames. The vectors describing the position of the body coordinate relative to the modal coordinate, r are shown in Fig. 4. The position of b/h,h the body coordinate relative to the modal coordinate, r , in Eq. (A-8a) is defined by the position of the b/h,h body frame, r , and the modal frame, r , relative to the Earth frame. Similarly, the orientation of the b/e,e h/e,h body coordinate relative to the modal coordinate, θ , is defined in Eq. (A-8b).

b/h,h T r = C r − r (A-8a) b/h,h b/e,e h/e,h e/h T θ = C θ (A-8b) b/h,h b/e,e e/h The rotation between the modal axis and Earth axis is fixed by the definition of the coordinate systems.

Therefore, it is simple to differentiate Eq. (A-8a) to calculate the velocity as shown in Eq. (A-9).

T ˙ r = C ˙ r − ˙ r (A-9) b/h,h b/e,e h/e,h e/h As the modal frame was described in section III.A, the motion of the modal frame relative to the Earth-fixed frame is a constant known parameter depending on the trim in Eq. (A-10).

{ } T := ˙ r (A-10) h/e,h − V cos ψ V sin ψ 0 0 0 0 0 The Euler angles for the rotation matrix, C , are not in the state vectors, Eqs. (A-4), (A-5), and (A-6).

e/w Therefore, the rotation matrix is expanded in Eq. (A-11) to multiple rotation matrices using the angles in Table 1.

T T T ˙ r = C C C v − ˙ r (A-11) b/h,h b/e,w h/e,h e/h b/e w/b The complete transformations also require the accelerations. The equation of Coriolis[ 4 ] is used to differentiate Eq. (A-11) to derive the translational accelerations in Eq. (A-12).

(( ) ) T T T ¨ r = C C C ω + ω × v + ˙ v (A-12) b/h,h w/b,w b/e,w b/e,w b/e,w e/h b/e w/b As with the translation, the orientation of the modal frame is fixed relative to the Earth frame, so differentiating Eq. (A-8b) is straightforward for the relative rotational rate in Eq. (A-13).

T ˙ ˙ θ = C θ (A-13) b/h,h b/e,e e/h Because the rotations for the Euler angles are not orthogonal, the translation of the reference frame of the angular rates from Earth to the body frame require the Euler rate matrix of Eq. (A-14)[4].

  1 0 − sin θ   H ( φ, θ, ψ ) =  0 cos φ sin φ cos θ  (A-14) b/e 0 − sin φ cos φ cos θ Using Eq. (A-14), the rate of the angle between the body frame and the modal frame can be expressed as shown in Eq. (A-15) as a function of the body axis.

T − 1 T ˙ θ = C H C ω (A-15) b/h,h b/e,w e/h b/e w/b From the properties of angular rates[ 4 ], the angular accelerations in the body frame and Earth-fixed frame are equal. Therefore, the expression of the rotational portion of the acceleration in Eq. (A-16) is straightforward.

T − 1 T ¨ θ = C H C ˙ ω (A-16) b/h,h b/e,w e/h w/b b/e This completes all the terms required for the transformations of Eq. (A-7).

The next step for the generation of the linear models is to linearize the transformation from the states in the modal frame to the mean axis frame. The nonlinear transformation of Eq. (A-7) is linearized about the wing level trim condition in Eqs. (A-17a) and (A-17b).

{ } T T η = (A-17a) 0 0 0 0 α 0 ξ b 0 0 0 { } T T ˙ η = (A-17b) 1 0 α 0 0 0 0 b 0 { } T T ¨ η = (A-17c) 0 0 0 0 0 0 0 b The linearized approximation is calculated using the Jacobian of the nonlinear transformation, evaluated at the wing level trim. The non-dimensionalization of Eqs. (A-2a) and (A-2b) are applied to give a non-dimensional transformation in Eqs. (A-18a) and (A-18b).

  ∗∗ [ ]   ( ) ( ) ∗∗ ˆ η − 1 ¯ c − 1 b η ≈ (A-18a) ∇ T ¨ η , ˙ η D ∇ T ¨ η , ˙ η D ¨ η accel ˙ η accel ∗ h b b u b b u b 0 0 2 V b 0 0 e   ˆ η b { } [ ] { } ( ) ( ) ∗ ∗ − 1 − 1 ∇ T ˙ η , η D ∇ T ˙ η , η D η ˙ η rate η rate ˆ η h b b u b b x b 0 0 b 0 0 b ( ) ≈ (A-18b) − 1 0 ∇ T η D η η disp ˆ η h b x b 0 b The Jacobian is evaluated to give the complete non-dimensional linear transformation in Eqs. (A-19a) and (A-19b).

{ } [ ] { } ∗ ∗ ˆ ˆ ∗ ∗ T ( α ) T ( α ) η 0 0 ˆ η h b hb he = (A-19a) ˆ η ˆ 0 T η h he b   ∗∗   [ ] ∗∗ ˆ η ˆ ˆ b T ∗∗ ∗ ( α ) T ∗∗ ( α ) η = 0 0 (A-19b) ∗ h h b h b   ˆ η b Where the matrices are defined by Eqs. (A-20a), (A-20b), (A-20c), (A-20d), and (A-20e).

  ¯ c − 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 − 0 0 0 0 [ ]   ˆ   T = (A-20a) 0 0 0 − 1 0 0 0 he     0 0 0 0 1 0 0       0 0 0 0 0 − 1 0 0 0 0 0 0 0 I   ¯ c − 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 − 0 0 0 0 [ ]   ˆ   T ∗ ( α ) = (A-20b) 0 0 0 − sec α 0 0 0 e   hb   0 0 0 0 1 0 0       0 0 0 − tan α 0 − 1 0 0 0 0 0 0 0 I   0 0 0 0 0 0 0   ¯ c ¯ c 0 0 0 − sin α 0 0   2 2   ¯ c   0 0 0 0 0 0 [ ]   ˆ   ∗ T ( α ) = 0 0 0 0 0 0 0 (A-20c) e   he   0 0 0 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0 0 0 0   ¯ c − 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 − 0 0 0 0 [ ]   ˆ   ∗∗ ∗ T ( α ) = 0 0 0 − sec α 0 0 0 (A-20d) e   h b   0 0 0 0 1 0 0       0 0 0 − tan α 0 − 1 0 0 0 0 0 0 0 I   0 0 0 0 0 0 0   ¯ c 0 0 0 0 0 0     ¯ c   0 0 0 0 0 0 [ ]   ˆ   ∗∗ T = 0 0 0 0 0 0 0 (A-20e)   h b   0 0 0 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 [ ] ˆ The body axis being a non-inertial coordinate system, the matrix, T ∗ , is nonzero; however, it is still he singular.

The time domain transformation of Eqs. (A-19a) and (A-19b) does fully describe the the relationship [ ] ˆ between the modal and body (mean) coordinate systems. The nonzero T ∗ matrix, however, prevents the he direct application to the aerodynamic model which is defined in the frequency domain. The time domain transformation must be translated into the frequency domain to be applied to the unsteady aerodynamic model. To determine the state transformation matrices in the frequency domain, the Fourier transformation to the non-dimensional frequency domain, Eq. (A-21), is applied to the time domain transformation, Eq.

(A-19a).

{ } ∗ η = ik { η } (A-21) h h Applying the transform to Eqs. (A-19a) and (A-19b), the non-dimensional reduced frequency response of the linear equation can be analytically determined in Eqs. (A-22a) and (A-22b).

{ } [ ] [ ] ∗ ˆ ˆ ik { η } = T ∗ ˆ η + T ∗ { ˆ η } (A-22a) h b b hb he [ ] ˆ { η } = T { ˆ η } (A-22b) he h b These equations can be solved to get the modal states, { η ( k ) } , as a function of the body axis velocities, { } ∗ ˆ η ( k ) , in Eqs. (A-23a) and (A-23b).

b { } ∗ [ ] ˜ ˆ { η ( k ) } = T ( k ) η ( k ) (A-23a) h b ( ) − 1 [ ] [ ] [ ] [ ] [ ] − 1 r ˜ ˆ ˆ ˆ ∗ ∗ T ( k ) = ik ik I − T T T ∀ k 6 = 0 (A-23b) r he he hb [ ] ˆ Because the matrix T ∗ is singular, the transformation in Eq. (A-23b) is undefined at zero reduced he frequency. Therefore, the case of zero reduced frequency must be handled separately.

The zero reduced frequency transformation is derived from the frequency domain kinematic equations in Eq. (A-22a) at zero reduced frequency, shown in Eq. (A-24).

{ } [ ] [ ] [ ] ∗ − 1 ˆ ˆ ˆ { 0 } = T ∗ ˆ η + T ∗ T { η } (A-24) he b h hb hb The solution to Eq. (A-24) is not unique, but a nontrivial is defined by the nullspace, kernel ( K ), of the [ ] ˆ singular matrix, T ∗ , in Eqs. (A-25a) and (A-25b).

e h { } ∗ [ ] ˜ { η ( k ) } = T ˆ η ( k ) (A-25a) b [ ( ) ] [ ] [ ] ˜ ˆ ˆ T = 0 K T ∗ T ∗ (A-25b) he hb The kernel is defined to ensure consistency, resulting in the steady-state transformation of Eq. (A-26).

  1 0 0 0 0 0 0   0 1 0 0 0 0 0       0 0 1 0 0 0 0 [ ]   ˜   T = (A-26) 0 0 0 sec α 0 0 0     0 0 0 0 0 0 0       0 0 0 tan α 0 0 0 r 0 0 0 0 0 0 I r The fifth and sixth columns of the transformation are zero because these modal displacements are directly correlated to the angle of attack and angle of sideslip. The remaining modal displacements are independent of the body axis velocity states, ˙ η .

b References [1] 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.

doi: 10.2514/6.2005-1947 [2] 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.

doi: 10.2514/6.2010-9350 [3] 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,” 2004.

[4] Stevens, B. L., and Lewis, F. L., Aircraft Control and Simulation , John Wiley and Sons, Inc., NJ, 2003.

[5] Milne, R. D., “Dynamics of the Deformable Aeroplane,” R.&M. No. 3345, Sept. 1962.

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

doi: 10.2514/6.1974-402 [7] Roger, K. L., “Airplane Math Modeling Methods for Active Control Design,” AGARD-CP-228, April 1977, pp. 4.1–4.11.

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

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

doi: 10.2514/3.45623 [10] MSC Nastran 2013 Dynamic Analysis User’s Guide , MSC.Software Corporation, Santa Ana, CA, 2013.

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

[12] Albano, E., and Rodden, W. P., “A Doublet-Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows,” AIAA Journal , Vol. 7, No. 2, 1969, pp. 279–285.

doi: 10.2514/3.5086 [13] Tiffany, S. H., and Karpel, M., “Aeroservoelastic Modeling and Applications of Using Minumum-State Approximations of the Unsteady Aerodynamics,” NASA TM-101574, 1989.

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

doi: 10.2514/6.2006-2025 [15] Dorf, R. C., and Bishop, R. H., Modern Control Systems , Pearson, Upper Saddle River, NJ, 2008.

[16] Lin, W.-F., and Clarke, M. D., “Factors Influencing the Accuracy of Aerodynamic Hinge-Moment Prediction,” AFFDL-TR-78-72, Aug. 1978.

[17] Giesing, J. P., Kalman, T. P., and Rodden, W. P., “Correction Factor Techniques for Improving Aerodynamic Prediction Methods,” NASA-CR-144967, 1976.

[18] Britt, R. T., Volk, J. A., Dreim, D. R., and Applewhite, K. A., “Aeroservoelastic Characteristics of the B-2 Bomber and Implications for Future Large Aircraft,” RTO MP-36, 1999.

[19] Ashley, H., “Some Considerations Relative to the Prediction of Unsteady Air Loads on Lifting Configurations,” Journal of Aircraft , Vol. 8, No. 10, 1971, pp. 747–756.

doi: 10.2514/3.44302 [20] Katz, J., and Plotkin, A., Low-Speed Aerodynamics: From Wing Theory to Panel Methods , Cambridge University Press, NY, 1991.

[21] Moreno, R., Narisetti, R., von Knoblauch, F., and Taylor, P. F., “A Modification to the Enhanced Correction Factor Technique to Correlate With Experimental Data,” AIAA 2015-1421, 2015.

doi: 10.2514/6.2015-1421 [22] Hodrick, R. J., and Prescott, E. C., “Postwar U.S. Business Cycles: An Empirical Investigation,” Journal of Money, Credit, and Banking , Vol. 29, No. 1, 1997, pp. 1–16.

[23] Tibshirani, R. J., and Taylor, J., “The Solution Path of the Generalized Lasso,” The Annals of Statistics , Vol. 39, No. 3, 2011, pp. 1335–1371.

doi: 10.1214/11-AOS878 [24] Ouellette, J., “Aeroservoelastic Modeling of Body Freedom Flutter for Control System Design,” AIAA 2017-0019, 2017.

doi: 10.2514/6.2017-0019 [25] Morelli, E. A., and Klein, V., Aircraft System Identification: Theory and Practice , Sunflyte Enterprises, Williamsburg, VA, 2016.