NASA/TM—2015–218875

Rapid State Space Modeling Tool for 

Rectangular Wing Aeroservoelastic Studies 

Peter M. Suh

Armstrong Flight Research Center, Edwards, California

Howard J. Conyers

Stennis Space Center, Mississippi 

Dimitri N. Mavris

Georgia Institute of Technology, Atlanta, Georgia 

 

Click here: Press F1 key (Windows) or Help key (Mac) for help 

September 2015
NASA STI Program ... in Profile

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space 

science. The NASA scientific and technical 

information (STI) program plays a key part in 

helping NASA maintain this important role.

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NTRS Registered

and its public interface, the NASA Technical 

Reports Server, thus providing one of the largest 

collections of aeronautical and space science STI 

in the world. Results are published in both nonNASA channels and by NASA in the NASA STI 

Report Series, which includes the following report 

types:

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA Programs and include extensive data 

or theoretical analysis. Includes compilations of significant scientific and technical 

data and information deemed to be of 

continuing reference value. NASA counterpart of peer-reviewed formal professional 

papers but has less stringent limitations on 

manuscript length and extent of graphic 

presentations.

 TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 

preliminary or of specialized interest, 

e.g., quick release reports, working 

papers, and bibliographies that contain 

minimal annotation. Does not contain 

extensive analysis.

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

 CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or 

co-sponsored by NASA.

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest.

 TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and 

feeds, providing information desk and personal 

search support, and enabling data exchange 

services.

For more information about the NASA STI 

program, see the following:

 Access the NASA STI program home page 

at http://www.sti.nasa.gov

 E-mail your question to help@sti.nasa.gov

 Phone the NASA STI Information Desk at 

757-864-9658

 Write to:

NASA STI Information Desk

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

This page is required and contains approved text that cannot be changed.
NASA/TM—2015–218875

Rapid State Space Modeling Tool for 

Rectangular Wing Aeroservoelastic Studies 

Peter M. Suh

Armstrong Flight Research Center, Edwards, California

Howard J. Conyers

Stennis Space Center, Mississippi 

Dimitri N. Mavris

Georgia Institute of Technology, Atlanta, Georgia 

 Insert conference information, if applicable; otherwise delete 

Click here: Press F1 key (Windows) or Help key (Mac) for help 

National Aeronautics and

Space Administration

Armstrong Flight Research Center

Edwards, CA 93523-0273

September 2015
 Enter acknowledgments here, if applicable. 

Available from:

NASA STI Program National Technical Information

Service Mail Stop 148 5285 Port Royal Road

NASA Langley Research Center Springfield, VA 22161

Hampton, VA 23681-2199 703-605-6000

Click here: Press F1 key (Windows) or Help key (Mac) for help 

This report is available in electronic form at http://www.sti.nasa.gov/ or http://ntrs.nasa.gov
1

Abstract

This report introduces a modeling and simulation tool for aeroservoelastic analysis of 

rectangular wings with trailing-edge control surfaces. The inputs to the code are planform design

parameters such as wing span, aspect ratio, and number of control surfaces. Using this 

information, the generalized forces are computed using the doublet-lattice method. Using Roger’s 

approximation, a rational function approximation is computed. The output, computed in a few 

seconds, is a state space aeroservoelastic model which can be used for analysis and control 

design. The tool is fully parameterized with default information so there is little required interaction 

with the model developer. All parameters can be easily modified if desired. 

The focus of this report is on tool presentation, verification, and validation. These processes

are carried out in stages throughout the report. The rational function approximation is verified 

against computed generalized forces for a plate model. A model composed of finite element plates 

is compared to a modal analysis from commercial software and an independently conducted 

experimental ground vibration test analysis. Aeroservoelastic analysis is the ultimate goal of this 

tool, therefore, the flutter speed and frequency for a clamped plate are computed using dampingversus-velocity and frequency-versus-velocity analysis. The computational results are compared

to a previously published computational analysis and wind-tunnel results for the same structure.

A case study of a generic wing model with a single control surface is presented. Verification of 

the state space model is presented in comparison to damping-versus-velocity and 

frequency-versus-velocity analysis, including the analysis of the model in response to a 1-cos 

gust.

Nomenclature

A state space matrix

AIC(ik̂ ) general aerodynamic influence coefficient matrix defined for reduced frequency, k

Ap diagonal matrix of the area of all aerodynamic panels

A0 ... A6 rational function approximation matrices

B control input matrix

c̅ mean aerodynamic chord

C global damping matrix

Ĉ generalized structural and aerodynamic damping 

C1

[0 Φ 0 0 0 0]

D(îk) aerodynamic influence matrix defined at control points of aerodynamic panels for k

exp exponential function

Fz normal force applied to the beam 

g damping

gwash gust wash

G gust design constant

H gust gradient

î imaginary number

k reduced frequency

K global stiffness matrix

K̂ generalized structural and aerodynamic stiffness

l lag constant iteration number

L total number of selected lag coefficients

m number of modes
2

M global mass matrix

M̂ generalized structural and aerodynamic mass

My torque applied to beam about the length of the beam

n number of degrees of freedom

na number of aerodynamic panels

P(t) external forcing function

q vector of modal coordinates

q̅ dynamic pressure

qs sensor measurement

Q(ik̂ ) generalized aerodynamic force matrix defined for a reduced frequency, k

Q̂ rational function approximation of generalized aerodynamic forces

s̅ defined as îk

s Laplace variable

t time

T lag constant

uc control surface rotation

ucom control command input 

V∞ freestream velocity

wc.p.

(ik̂ ) normalwash vector defined at the control points of the aerodynamic panels for reduced 

frequency, k

wg gust velocity

wg,0 design gust velocity

wg,max maximum gust velocity

wg(t) gust velocity as a function of time

Wc.p. matrix of wc.p.

(ik̂ ) collected for all m mode shapes

x state vector

x1 state vector of generalized modal displacements 

x2 state vector of generalized modal velocities

x3 ... x6 aerodynamic lag state vectors 

xc.p. vector of free-stream coordinates at the control points of each aerodynamic panel 

defined from the leading edge of the wing

xgust ̅ length from the leading edge of the wing to the gust origin

X chordwise axis

y measurement vector

Y spanwise axis

zc.p. out-of-plane z deformation of normal mode shape defined at control points of 

aerodynamic panels

Zf matrix of zc.p. collected for all m mode shapes

Z normal axis

βll

th

lag constant

ξ damping ratio

θc.p. slope of the normal mode shape defined at the control points of the aerodynamic panels

Φ modal matrix

ω circular frequency

(∙)

̇ 1

st derivative in time of the argument

(∙)

̈ 2

nd derivative in time of the argument

(∙)

T

transpose of the argument

(̃∙) generalized form of the argument
3

(∙)c control partitioning of argument

(∙)c# control partitioning of argument for specific #

(∙)g gust partitioning of argument

(∙)g# gust partitioning of argument for specific #

List of Acronyms

DLM doublet-lattice method

DOF degrees of freedom

FEM finite element model

GAF generalized aerodynamic force

GVT ground vibration test

ODE ordinary differential equation

RFA rational function approximation

V-f frequency versus velocity

V-g damping versus velocity

VLM vortex-lattice method

Introduction

Aeroservoelastic modeling is often performed using medium-fidelity tools such as MSC 

Nastran (MSC Software Corporation, Santa Ana, California) (ref. 1) and ZAERO (ZONA 

Technology Inc., Scottsdale, Arizona) (ref. 2). In transonic regimes, the properties of the flow are 

often modeled by high-fidelity computational fluid dynamics tools such as the National 

Aeronautics and Space Administration (NASA) -developed CFL3D (NASA Langley Research 

Center, Hampton, Virginia) (ref. 3) and CFD-FASTRAN (ESI Group, Paris, France) (ref. 4). The 

aforementioned tools are well known solution finders for academic and industry aeroelastic

problems. The aeroservoelastic tool development presented herein approaches the level of 

fidelity of the MSC Nastran aeroservoelastic code. The code is currently restricted to the analysis 

of rectangular wing models; its original purpose was to model fiber optic sensor feedback in a 

flutter-unstable aeroservoelastic environment (ref. 5). Fiber optic sensors can be used to feed 

back 1,000s of strain or shape measurements on a wing to a flutter suppression or shape 

controller. Investigating the nature of this sensor-control paradigm was a major thrust of Suh’s 

PhD Thesis (ref. 6). This research was also conducted in support of concepts to be tested on the 

Lockheed Martin (Bethesda, Maryland) X-56A Multi-utility Aeroelastic Demonstrator unmanned 

air vehicle (ref. 7).

For iterative analytical investigations, a tool was developed which could rapidly convert 

information such as the number of control surfaces, the wing span, and the chord length to a 

medium-fidelity finite-element-based aeroservoelastic model. The model is able to capture

low-speed aerodynamic structural interactions with good detail in unsteady aerodynamics and 

finite-element model stiffness, mass, and damping interactions. The tool models the aerodynamic 

effect of trailing-edge control surfaces. 

This medium-fidelity tool supports 12 degrees of freedom (DOF) plate- (ref. 8) and 6-DOF 

beam-based isotropic finite elements and was coded into the MATLAB® (MathWorks, Inc., Natick, 

Massachusetts) environment. The unsteady aerodynamic forces were derived from the

doublet-lattice method (DLM) (refs. 9 and 10) taking advantage of recent research on 

improvements to this method. The steady forces were modeled with the vortex-lattice method 
4

(VLM) (ref. 11). To transport the force coefficient information into a state space environment for 

control design, a rational function approximation (RFA) method (ref. 12) was employed.

This report includes an overview of the tool functionality and a discussion of the components 

in the tool as well as the verification of some of the tool components. The beam model verification 

is given first, followed by verification of the RFA and a comparison of the model outputs to 

experimental data published by Conyers et al (ref. 13). The process of validation adds further 

confirmation to Conyers et al.’s original analyses and supports their work. After validation a wing 

model is presented and qualitatively compared to the flutter analysis methods used in the tool. 

The presentation of this wing model also demonstrates the breadth and depth of the tool. It is 

hoped that the tool will become available to any student of aeroservoelasticity. The following 

section provides an overview of the functionality of the tool.

Tool Functionality

The tool has three purposes. The first purpose is to make it easier for someone new to the 

subject to study aeroservoelasticity with a higher level of detail than that offered by an airfoil 

model. The second purpose of the tool is to have a compilation of the required tools in a single 

coding environment. The third purpose of the tool is to support rapid model deployment by

speeding up the model development process. 

Whereas many structural modeling tools require extensive model design and setup time, this 

tool is fully parameterized. By inputting only two parameters, aspect ratio and half span, a full 

aeroservoelastic state space model can be generated within seconds. Anticipating that design 

changes, such as the number and placement of control surfaces, the range of reduced 

frequencies, or the actuator parameters, may be desired, these detailed parameters, while set to 

default parameters, are also modifiable. The model can then be rerun with small changes to study 

the effects of the changes in the final state space design. 

The tool development is deployed in three key modules: geometry, aerodynamics, and 

analysis. The inputs and full functional capability of the aeroservoelastic analysis tool for 

rectangular wings are presented in figure 1.
5

Figure 1. Wing model design and aeroservoelastic analysis tool.

The three modules in figure 1 showcase the aeroservoelastic design tool. Arrows point to the 

flow of computation in the tool. The geometry module is highly parameterized. Two spars of 

adjustable material type are assumed to lie along the one-fourth-chord and at control surface-wing 

connection points. The cross-sectional dimensions can be tapered down toward the wing tips if it 

is so desired. The ribs are spaced at every panel width; cross-sectional dimensions can also be 

tapered. Control surface finite beam elements can also be adjusted in the same manner. The 

material properties of the wing are adjustable.

The aerodynamics module is parameterized by a reduced frequency range and number of lag 

states, utilizing the geometry of the wing specified in the finite element model (FEM) module for 

aerodynamic panel design. The lag constants are nominally selected from the low-frequency 

range. The aerodynamics module generates aerodynamic modes from the normal modes and a 

gridded aerodynamic panel layout. The aerodynamics module computes the generalized 

aerodynamic force (GAF) (refs. 14 and 15) and the corresponding RFAs. A direct output of the 

aerodynamics module is the aerodynamic steady and unsteady matrices.

The analysis module is useful for several purposes. Using the RFA, the V-g and V-f analysis 

are computed for the resulting aeroservoelastic model. These plots are useful for predicting the 

flutter speed and frequency. These charts are useful for the predicting the properties of the state 

space matrices. The state space matrices, which are parameterized by the actuator dynamics, 

can be used for most forms of control design and simulation with accelerometer feedback.

The aforementioned modules are initialized with default values, making this a rapid design 

process and removing much of the modeling effort toward the case of one being purely interested 

in understanding aeroservoelastic phenomena. For modeling a specific wing, the default 

parameters can be adjusted accordingly. The following section details the aeroservoelastic state 

space model development which the simulation tool was designed to support. The sections below

also discuss specific modeling choices in the tool, where appropriate. 
6

Aeroservoelastic State Space Model

The aeroservoelastic model equations of motions used in the tool are built up using notation

and method similar to that found in reference 16. A Lagrange representation (ref. 17) was chosen 

to model the wing perturbation dynamics. The aeroelastic matrix equation of motion used to model 

the structure is presented as shown in equation (1):

Mẍ+ Cẋ + Kx + q̅AIC(ik̂ )x = P(t) (1)

where M ∈ Rn×nis the global mass matrix, C ∈ Rn×nis the global damping matrix, K ∈ Rn×nis 

the global stiffness matrix, q̅ is the dynamic pressure, AIC(îk) is the aerodynamic influence

coefficient matrix computed at the flight condition and for reduced frequency value k =

ωc̅

2V∞

, c̅is 

the mean aerodynamic chord length, ω is the circular frequency, V∞ is the freestream velocity, 

x ∈ Rn×1is the physical displacement vector of rotational and translational DOF, and P(t) is the 

external forcing function. 

Equation (1) must be transformed to a modal representation in order to reduce the number of 

DOF of the system representation for large finite-element models. This is accomplished from an 

eigenvector-based approach as follows. The free vibration of the unforced undamped system in 

equilibrium is given as shown in equation (2) (ref. 18).

Mẍ+ Kx = 0 (2)

The eigenvalue solution of the system produces the natural frequencies and eigenvectors (dry or 

in vacuo mode shapes, Φ) of the system. The following linear transformation of state is formed 

as shown in equation (3):

x = Φq (3)

where Φ ∈ Rn×m is a matrix of eigenvectors corresponding to selected elastic and control modes, 

q. Pre-multiplying equation (1) by ΦT and substituting equation (3), the matrix equation of motion 

is transformed as shown in equation (4).

ΦTMΦq̈+ ΦTCΦq̇ + ΦTKΦq + q̅ΦTAIC(ik̂ )Φq = ΦTP(t) (4)

Equation (4) is rewritten in the generalized form shown in equation (5).

M̃q̈+ C̃q̇ + K̃q + q̅Q(ik̂ )q = P̃(t) (5)

Equation (5) is in an inappropriate form for state space modeling. Aerodynamic forces are 

represented at specific frequencies and external forces are represented as varying with time. The 

aerodynamic forces are transformed through an RFA as described in the next section.

Generalized Aerodynamic Forces and Rational Function Approximation

To model the generalized aerodynamic forces, an RFA or conversion to s plane of the 

tabulated GAF matrices, Q(ik̂ ), must be completed. The GAF matrices are calculated at specific 

reduced frequencies k through a doublet-lattice procedure with aerodynamic mode shapes and 

modal shape derivatives. To improve the accuracy of the RFA for higher Mach numbers, 

Rodden’s quartic kernel approximation (ref. 19) is utilized in this tool. In addition, the number of 
7

aerodynamic panels in the chordwise direction is constrained to be greater than 50k

πc̅

⁄ for 

generalized force matrix computations (ref. 20). This is a suggested aerodynamic panel criterion 

from Rodden et al. (ref. 20) for the case of lifting surface theory.

Before computing the RFA directly, the normalwash at the control points of the aerodynamic

panels must be computed. The normalwash (ref. 21) is required for doublet lattice and GAF 

computations. The normalwash is computed for each mode from deflections and rotations at the 

control points (centered at 3/4 of the length of the aerodynamic panel in the chordwise direction 

from the leading edge). The vector normalwash is defined at the control points of all aerodynamic

panels as shown in equation (6):

wc.p.

(ik̂ ) = θc.p. + î

2k

c̅

zc.p.

(6)

where θc.p.is the interpolated slope in the chord-wise direction of the aerodynamic panel defined 

at the control points and zc.p.is the interpolated deformation defined of the aerodynamic panel 

defined at the control points.

To compute the normalwash, interpolation of the mode shapes at control points of aerodynamic

panels is required since deformations and rotations are defined only at the structural nodes. This 

interpolation is completed here with a distance-weighted average approach using the deformation 

points at the structural nodes. The GAF matrices are next computed at each reduced frequency 

using the relationship given in equation (7) (see ref. 15):

Q(îk) = Zf

TD(îk)

−1ApWc.p.

(7)

where Zf ∈ Rna×m is a matrix of all normal deformations centered at the force point (quarter chord) 

on the aerodynamic panels, Wc.p. ∈ Rna×m is the matrix of normalwash vectors, for all 

aerodynamic mode shapes, wc.p.

(ik̂ ), D(îk) ∈ Rna×na is the aerodynamic influence matrix defined

at control points of all panels for reduced frequency k, and Ap ∈ Rna×na is a diagonal matrix of

the area of all aerodynamic panels. The matrix D(ik̂ ) has both a steady component which is here 

computed by VLM and an unsteady component computed with DLM as shown in figure 1.

With tabulated Q(îk) computed, there are several ways of computing an RFA. Here are two: 

Roger’s RFA (ref. 12) and Karpel’s minimum-state method (ref. 22). For the current development

and its simplicity, Roger’s RFA is chosen for rendering the RFA from the GAF matrices. To 

implement the RFA, the lag constants βl,l = 1 ... L are automatically selected from the lower 

range of the predefined set of reduced frequencies. For fitting of the generalized force matrices, 

the tool allows for one to four lag states.

The RFA of the generalized unsteady aerodynamic forces in the Laplace domain may be 

written as shown in equation (8):

Q̂(s̅) = A0 + s̅A1 + s̅

2A2 +∑

s̅

s̅+ βl

A2+l

L

l=1

(8)
8

where s̅= s

c̅

2V∞

= î

ωc̅

2V∞

= ik̂ , c̅is the mean aerodynamic chord, and s is the Laplace variable. The 

matrix coefficients may be found through a least-squares approximation of a set of GAF matrices. 

Oftentimes, constraints must also be applied to get good matching at lower reduced frequencies

(ref. 23). Constraints are not currently used in this tool, restricting it to high-frequency analysis.

The next section gives an overview of how the RFA is incorporated into the state space model.

Aggregate State Space Model

In order to get to a time-based representation of the state space model, the RFA must be 

replaced in the equations of motion (refs. 16 and 24). Therefore by inserting equation (8) into

equation (5), the equations of motion are rewritten in Laplace form as shown in equation (9):

s

2M̃q̈+ sC̃q + K̃q + q̅ [A0 +

sc̅

2V∞

A1 + (

sc̅

2V∞

)

2

A2 + A3x3 + A4x4 + ⋯ ] q = 0 (9)

where the lag states are defined as shown in equation (10).

x2+l ≜

sq

s +

2V∞

c̅

βl

; l = 1 ... 4

(10)

Equation (10) may be rewritten and inverse Laplace transformed resulting in matrix ODEs of the 

form given in equation (11).

ẋ2+l +

2V∞

c̅

βlx2+l = q̇; l = 1 ... 4 (11)

Like terms are grouped in equation (9) and the equation may be rewritten after an inverse 

Laplace transform as in equation (12).

(K̃ + q̅A0)q + (C̃ + q̅

c̅

2V∞

A1) q̇ + (M̃ + q̅ (

c̅

2V∞

)

2

A2) q̈+ q̅A3x3 + q̅A4x4 + ⋯ = 0 (12)

Condensing variables in equation (12) results in the simplified matrix ordinary differential equation 

(ODE) given in equation (13).

K̂q + Ĉq̇ + M̂q̈+ q̅A3x3 + q̅A4x4 + ⋯ = 0 (13)

Equations (11) and (13) represent a set of matrix ODEs and may be converted to state space 

form. Before doing so, it is important to understand that the modal coordinates may include control 

modes, rigid body modes, elastic modes, and gust modes. The mass, stiffness, aerodynamic lag 

and damping matrices are partitioned accordingly. Assuming only the presence of elastic modes, 

control modes, and gust modes in the analysis, the exact matrix state space equation used in the 

tool assuming four lag states are modeled is given in equation (14) (see ref. 24).
9

ẋ =

{

 

 

 

 

ẋ1

ẋ2

ẋ3

⋮

ẋ6}

 

 

 

 

=

[

 

 

 

 

 

 

0 I 0 ... 0

−M̂−1

(K̂ Ĉ q̅I ... q̅I)

0 A3 −β1(

2V∞

c̅

) 0 0

⋮ ⋮ 0 ⋱ 0

0 A6 0 0 −β4(

2V∞

c̅

)

]

 

 

 

 

 

 

{

 

 

 

 

x1

x2

x3

⋮

x6}

 

 

 

 

+

[

 

 

 

 

0 0 0

−M̂−1

(K̂c Ĉ

c M̂c)

0 Ac1 0

⋮ ⋮ ⋮

0 Ac4 0 ]

 

 

 

 

{

uc

u̇c

üc

} +

1

V∞

[

 

 

 

 

0 0

−M̂−1

(K̂g Ĉ

g)

0 Ag1

⋮ ⋮

0 Ag4]

 

 

 

 

{

wg

ẇg

}

(14)

The vector states x1 and x2 replace the generalized modal state displacement and velocity vectors 

respectively. The remaining states x3, ... x6 are four orders of aerodynamic lag states. The control 

surface states uc, u̇c, üc may be modeled with actuator transfer function dynamics. This also is 

required to convert these input states into a control command form required for control design. 

The gust velocity and acceleration are represented by wg and ẇg, respectively. The following 

section overviews how the sensors are modeled within the tool environment. 

Sensor Outputs

Acceleration measurements at any place on the wing model are modeled assuming four

aerodynamic lag states as shown in equation (15):

q̈s = Φq = [0 Φ 0 0 0 0]ẋ = C1ẋ (15)

Pre-multiplying equation (14) by C1 results in equation (16).

y = q̈s = C1ẋ = C1Ax + C1Bucom (16)

There also exists the capability to model modal coordinate outputs directly, in the case that fiber 

optic sensors and a modal filter (see ref. 5) are utilized as in equation (17).

y = [I 0 0 0 0 0]x (17)

The significance of this type of sensor output is that it supports the concept of virtual deformation 

control (ref. 7). The next section overviews the actuator modeling capabilities of the tool. 

Actuators

For each control surface, a 2nd order actuator is utilized. To model command and flight 

computer delay a 1st order lag filter is multiplied with the 2nd order actuator transfer function. This 

method has an additional effect of removing direct feedthrough from the sensor output matrix, 

when accelerometers are used for feedback. The 3

rd order actuator model from command to 

surface movement is shown in equation (18):
10

uc

ucom

=

1

Ts + 1

(

ω

2

s

2 + 2ωξs + ω2

) (18)

where T is a scalar time constant, ω is the circular frequency, and ξ is the damping ratio. 

The actuator used in the results of this report is given here. The time constant is set to .02 s; 

ω was set to 74 rad/s; and ξ is set to 0.58. A step of the actuator with and without command lag 

is given in figure 2.

Figure 2. Actuator model.

The command lag for this case removes much of the overshoot of the 2nd order actuator alone,

as shown in figure 2, allowing the tool to model some conservatism in control. Additionally, the

high-frequency roll-off allows for more control authority at higher bandwidth, which is useful for 

structural control. The following section describes the gust modeling capabilities of the tool.

Gust Model

Gusts are important to test the wing model for performance in disturbance rejection. The gust 

wash is found by modeling the phase lag between individual panels and the gust origin (ref. 25). 

For a coordinate system with origin at the leading edge of the model where x increases in the 

trailing-edge direction, the gust wash is given as shown in equation (19):

gwash = −exp (î

2k

c̅

(xc.p. − xgust ̅ )) (19)

where xc.p.is the vector of stream-wise coordinates at the control points of each aerodynamic 

panel and xgust ̅ is the distance from the origin to the gust. If xgust ̅ is specified to be zero, the gust 

begins to build intensity at the leading edge of the wing. For gust simulations here, it is specified 
11

to be zero. A non-zero distance may result in large spirals in the GAF matrices (ref. 2). This spiral 

nature is difficult to capture with the RFA and is especially problematic for models attempting to 

capture large reduced frequencies.

Once the gust GAFs are computed with the gust wash in place of the normalwash as in 

equation (7), the gust RFA is computed as before (see eq. (8)) in a separate calculation. In this 

way, the gust model is added to the state space equations as shown in equation (13). The inputs 

to the state space model are driven with a temporal 1 − cos gust model. The standard temporal 

variation gust model (ref. 26) used in the tool is given as in equation (20):

wg(t) =

wg,0

2

(1 − cos (

πV∞

H

t)) (20)

where wg,0 is the design gust velocity, and H is the gust gradient, or half the distance of the total 

gust span. A maximum gust velocity and acceleration can be specified at the design operating 

condition. The gust velocity is specified as shown in equation (21):

wg(t) =

wg,max

2

(1 − cos(Gt)) (21)

where G is a gust design constant and wg,max is the maximum gust velocity. A derivative in time 

of equation (21) leads to the gust acceleration model shown in equation (22).

ẇg(t) =

wg,max

2

sin(Gt)G (22)

By inspection, it is clear that the maximum value occurs when sin(Gt) = 1, and thus the constant 

G is chosen as shown in equation (23).

G = 2

ẇg,max

wg,max

(23)

The ratio of the maximum acceleration and velocity determines the frequency of the gust. The

ratios are defined in the above way for physical intuition during gust testing. Using this model of 

a gust, a maximum velocity of 1 m/s and a maximum acceleration of 1 m/s

2 can be specified. 

The gust model resulting from these settings is given in figure 3.
12

Figure 3. 1 − cos(x) gust model for control design and simulation.

The gust model presented in figure 3 can be used for a realistic control design study. The input 

velocity and acceleration directly influence the state derivatives in equation (14). The proper 

maximum amplitudes of the gust depend on the theoretical operating design conditions. The 

particular gust model presented in figure 3 is used later in this report for wing model gust response 

testing. The next section introduces the verification and validation procedures taken during tool 

development.

Simulation Verification and Validation

Some aspects of the aeroservoelastic modeling tool are tested against published material or 

theoretical relationships. First, it is shown that the deflected FEM of the beam matches cantilever 

beam deflection theory results. Second, the RFA is verified using GAF information for the plate 

model. Lastly, the modal frequencies and flutter speed prediction modules are compared with that 

of an experimental model. First, the beam verification results are presented.

Cantilever Beam Verification

To build a wing model with control surfaces, a wing box or spars and ribs must be incorporated. 

It was decided that this information could be captured with 6-DOF isotropic beams. Here it is 

demonstrated that the FEM stiffness relations of the modeled beams match cantilever beam 

theory results. 

For the test, two beams with identical physical properties are simulated under the same static 

loading conditions of {Fz = −100N, My = −100N − m} at beam tip. The force Fzis applied in the 

normal direction, and My is a torque applied about the length of the beam at the tip. The 

dimensions of the beam are 1000mm x 500mm x 6.35 mm, with a built-in (clamped) boundary 

condition at one end. 
13

The first beam is analyzed as a simple continuous isotropic cantilever beam. The second beam 

is analyzed as a finite element beam, using the tool’s FEM. For this analysis the beam was 

discretized lengthwise into 30 finite isotropic beam elements. The results of the test are shown in 

figure 4.

Figure 4. Beam model verification, {force, moment} = {-100N, -100N-m} at beam end.

Figure 4 confirms that the FEM beam deflections and rotations overlaid the theoretical

cantilever beam deflections and rotations, giving confidence that the beam stiffness properties of 

the structure are modeled correctly. Correct modeling is very important in the generation of 

accurate mode shapes for the modal filtering process. The FEM beam mass matrices are 

modeled in a similar way. The next section overviews the verification of the RFA used in the tool.

Rational Function Approximation Verification

Verification of the RFA used here is completed for a clamped 0.3048m x 0.1524m x 0.001588m

polybicarbonate plate with published findings on its aeroelastic characteristics (ref. 13). In

Conyers et al.’s work the clamped plate was analyzed by discretizing it into a 16 × 16 layout of 

shell elements. It is similarly discretized with the tool presented here; the aerodynamic panel 

spans coincide with the FEM plate element spans. 

What is different from Conyers et al.’s analysis is that the number of aerodynamic panels in 

the chord-wise direction is allowed to vary with reduced frequency, as discussed previously

(ref. 20). This method represents a significant difference in the computational analysis published 

previously on the test article. Another difference is that the aerodynamic-structural matrix 

computations required for converting normal modes to normalwash are computed with 

distance-weighted averaging rather than surface interpolation. 
14

For the RFA analyses of the clamped plate the lag coefficients Β1 and Β2 were chosen by the 

tool to be 0.1889 and 0.3778, corresponding to two aerodynamic lag states. The generalized 

forces were computed at a set of pre-selected reduced frequencies given by {0, 0.2, 0.4, 0.6, 0.8, 

1.0}. For these reduced frequencies, a set of corresponding m × m GAF matrices were generated, 

where m is the number of modes in the model. The GAF matrices Q(ik̂ ) (see eq. (7)) were fit to 

a RFA, Q̂(s̅) as described in equation (8). The GAF coefficients corresponding to the 1st bending 

and 1st torsion mode, and the corresponding RFA evaluations, are plotted together in figure 5.

Figure 5. Rational function approximation: a) bending GAF coefficient and RFA; and b) torsion 

GAF coefficient and RFA.

For both GAF coefficients presented in figure 5, the RFA overlay is very good. The error has 

the characteristics of a least-squares error, which is typical since it is solved as a least-squares 

problem (ref. 12). The good fit does not verify that the GAFs themselves are correct, however;

this topic is addressed in the next section. But the overlay does indicate that enough aerodynamic

lag states have been modeled to capture the GAF variation with reduced frequency and that the 

RFA has been computed reasonably well. This is a very important step before state space 

modeling.

Modal Analysis Validation Study

An analytical and experimental modal frequency analysis was conducted at Duke University 

by Conyers et al. (ref. 13) The simulation code modal frequency predictions are compared to that 

from both the ANSYS (ANSYS, Inc., Canonsburg, Pennsylvania) (ref. 27) code (which they set 

up and computed) and experimental ground vibration test (GVT) measurements. For the current 

tool, the modal frequencies were computed using a standard eigenvalue analysis for the first five 

modes of the experimental plate model. A comparison of the results is presented in table 1.
15

Table 1. Modal frequency code comparisons and experimental results.

ANSYS (ref. 4)1

frequencies, Hz

Tool FEM

frequencies, Hz

Conyers et al. 

GVT, Hz (ref. 3)

Mode # 1 3.99 3.99 4.13

Mode # 2 16.96 16.97 17.24

Mode # 3 24.86 24.89 24.38

Mode # 4 55.33 55.40 54.25

Mode # 5 69.84 69.92 69.00

 1 Republished data with permission of Dr. Howard Conyers. Experiments conducted at Duke University.

The frequencies of the modes computed with each method were nearly identical in table 1. 

The frequencies computed by the tool nearly matched those computed by the ANSYS code, within 

0.2 percent error. The GVT results from Conyers et al. show frequencies that were slightly higher 

for the first two modes, but still within 4 percent error. The GVT results were slightly lower than 

the tool frequencies for the last three modes, within 3 percent error. No significant differences 

were noted overall. This experiment was considered critical since the skin of the wing and control 

surface was modeled this way. The FEM build-up technique is also verified in this way. The next 

section overviews the most critical validation step for the tool: the flutter evaluation.

Flutter Validation Study

Conyers et al. (ref. 13) also conducted an analytical and experimental investigation of the flutter 

properties of the plate. To perform this test, they developed their own flutter code, and to validate 

their code, they conducted an experimental flutter test in a wind tunnel. The flutter speed and 

frequency comparisons are presented in table 2.

Table 2. Flutter code comparisons and experimental results.

Conyers et al. 

flutter code (ref. 3)

Tool 

flutter code

Conyers et al. 

wind tunnel results (ref. 3)

Flutter speed, m/s 20.8 19.9 20.05

 Flutter frequency, Hz 10.3 10.9 11.50

Table 2 shows that the theoretical flutter frequency and speed calculated from this tool is very 

close to both Conyers et al.’s flutter code and that from their wind-tunnel experiment. The tool 

flutter speed was computed within 1 percent error of the experimental results. Another interesting 

aspect of the tool is that its computed flutter speed and frequency were slightly closer to the 

experimental results than those of the in-house flutter code that Conyers et al. used. It is 

postulated that one reason could be due to the increased aerodynamic paneling at higher reduced 

frequencies used in the tool. Another reason could be differences in computing the normalwash 

and structural-to-aero transference. This tool used a distance-weighted averaging scheme,

whereas Conyers et al. used surface interpolation.

The closeness of the flutter prediction values indicates that doublet-lattice and vortex-lattice 

computations in the model are good. The resulting GAF computations must also be reasonably 
16

good. Although only plate elements were used for this analysis, the difference of adding structure 

with beams will not change the aerodynamic model. Only the normalwash will change, since the 

deflections of the mode shapes will change. This validates the tool not only for plate analysis, but

also for wing model analysis, where beam structure is incorporated.

A discussion of how the numbers in table 2 were computed is needed. Conyers et al. computed 

their flutter speed and frequency with the V-g and V-f analysis (refs. 13 and 26) respectively, but 

did not present the plots. This information is presented in the next section to demonstrate the 

tool’s analysis methods while assisting with verification of their analytical work. 

Computation of the Flutter Speed and Frequency

Within Conyers et al.’s work, 10 modes were modeled in the aeroelastic analysis to find the 

theoretical flutter speed and frequencies of the rectangular plate without a hole, but a V-g and V-f 

plot was not given there. The V-g plot for the plate model computed using the RFA is given in 

figure 6.

Figure 6 shows 10 modes, but only the 1st bending and 1st torsion mode interact strongly and 

lead to flutter at the dashed vertical line. The flutter speed is computed where the 1st torsion mode 

crosses the 0 damping line. The flutter speed is computed in this way to be 19.9 m/s. At 26 m/s, 

a divergent mode appears, and this is verified with the V-f analysis presented in figure 7. 

Figure 6. V-g plot of theoretical plate model modes.
17

Figure 7. V-f plot of theoretical plate model modes.

Only two interacting modes are shown in figure 7. The first two modes at 0 airspeed begin at 

their dry natural frequencies presented in table 1. The 1st torsion modal frequency begins to shift 

with airspeed, and at the flutter speed of 19.9 m/s, the frequency of the 1st torsion mode (the 

mode which crossed through 0 damping) is 10.9 rad/s. The frequency of the 1st bending mode 

actually goes to 0 at 26 m/s, indicating that this is a divergence speed (see ref. 26). The model 

flutters before this speed, however, so this occurrence is inconsequential. 

The results indicated in this section support Conyers et al.’s experimental work. The next major 

section overviews a case study, wherein the tool is qualitatively verified for control design in 

aeroservoelasticity. 

Computational Wing Case Study

Much of the material presented so far has been for simple models with either theoretical or 

experimental results available. The following work does not have the luxury of available 

experimental or theoretical data validation case studies, so it is presented in such a way that each 

result confirms the previous result. There is some independence between each step, therefore 

this process is referred to as a qualitative verification procedure. First, the geometrical layout of 

the wing model is presented.

Geometrical Layout

The wing model is taken from a previous research paper, where its detailed spar, rib, and wing 

skin structural specifications can be found (ref. 5). The computational model is a 3.354m by 

0.838m half-span wing made completely from 6061-T6 aluminum. It is discretized with 16 finite 

element plates in the Y direction and 10 plates in the X direction. The single control surface is 

made up of three plates in the Y direction and three plates in the X direction. To reinforce the 

structure skin, two spars are added: one at the quarter-chord and one at the connection point of 

the control surface. Sixteen ribs runs through the entire wing chord-wise. The control surface is 
18

also supported with beam finite elements and is attached with two springs to the structure. The 

wing is clamped at all structural nodes at the wing root and free at all over points. The wing model 

is presented in figure 8.

Figure 8. Half-span wing model with one control surface and one accelerometer.

One difference between this model and the model presented previously (ref. 5) is that this 

model only uses one trailing-edge control surface. It also uses only a leading-edge wing tip

acceleration measurement. This technique renders this model a single-input single-output (SISO) 

system. In fact, many control surfaces and locations can be chosen. Sensors can be placed at 

any structural node point and multiple-input multiple-output (MIMO) systems can be formed. The 

mode shapes for the present structure are presented in the next section.

Mode Shapes

The mode shapes of the wing model presented in figure 8 were computed with an eigenvalue 

analysis. The mode shapes are made orthogonal to the mass matrix, which satisfies a critical 

mean axis constraint (ref. 17). Ten mode shapes were computed; the first four mode shapes are 

presented in figure 9.
19

Figure 9. Modal representation of wing model: a) 1st bending; b) 1st torsion; c) 2nd bending; and 

d) 2nd torsion.

The mode shapes in figure 9 take the usual form of what one would expect from the first four 

modes of a rectangular clamped structure. The frequency of each mode progressively increases 

with mode order as expected. The 1st bending mode is at 16 rad/s, and the 1st torsion mode 

occurs at 43.4 rad/s. The 2nd bending mode occurs at 85.6 rad/s, and the 2nd torsion mode occurs

at 184.6 rad/s. 

The normal modes are computed in a control-fixed position. That is, the torsional stiffness of 

the springs connecting the respective control surface to the wing is set to a very large value so 

that the control surface does not rotate. This is similar to clamping the control surface to the wing 

during a GVT.

A control-free analysis is completed to compute the control mode shapes. Control-free refers 

to reducing the torsional stiffness of the springs connecting the wing to the control surface before 

modal analysis. This method reduces the stiffness coupling due to control surface rotation which 

could lead to unwanted twisting motion of the wing during control surface rotation. The control 

mode is defined by a 1-deg boundary condition which is imposed on the control surface spring

attachment to the wing (see ref. 24).

With both control modes and normal modes defined in this way, the aerodynamic coupling of 

the modes can be analyzed. The next section presents the aerodynamic coupling of modes within

the structure. 

Wing Model Flutter Analysis

The modes of the wing model structure are coupled through aerodynamic forces, which is

apparent in both the V-g and V-f analyses. First the V-g analysis of the wing model is presented 

in figure 10.
20

Figure 10. V-g analysis of wing model.

Figure 10 indicates the variation of 10 modes with velocity. Only the 1st bending and 1st torsion 

mode interact strongly and lead to flutter at the dashed vertical line. As before, the flutter speed 

is computed where the 1st torsion mode crosses the 0 damping line. The flutter speed is found in 

this way to be 76.5 m/s. The characteristics of the aerodynamic damping changes are reminiscent 

of the V-g plot of the composite plate model shown in figure 6. The major difference is that 

divergence is not reached in the given airspeed range. The V-f analysis is presented in figure 11.

Figure 11. V-f analysis of wing model.
21

As before (see fig. 7), only the two interacting modes are shown in figure 11. The first two 

modes start off at their dry natural frequencies as shown in figures 9(a) and (b). The 1st torsion 

modal frequency begins to shift with airspeed, and at the flutter speed of 76.5 m/s the frequency 

of the 1st torsion mode is 28 rad/s. The V-g and V-f analysis results must now be confirmed with 

the corresponding state space model analysis presented in the next section.

State Space Model Analysis

The state space model was generated for the wing model as shown in equation (14). The 

following sections overview first the stability analysis of the model. The control characteristics, 

which are important for control design, are then overviewed. Finally the wing model is subjected 

to the gust model developed above. 

Stability Analysis

The stability analysis is conducted by plotting the plant poles versus airspeed. This process is 

critical because the poles of the plant affect stability. Figure 12 shows the wing model plant pole 

migration. 

Figure 12. Wing model plant pole migration with airspeed from 20 m/s to 78 m/s.

The pole migration with airspeed as shown in figure 12 indicates that the plant becomes 

neutrally stable (at the flutter boundary) at 76.5 m/s as was predicted in the V-g analysis shown 

in figure 10. Also as expected, the torsion mode is the one which crosses the stability boundary 

first. The bending mode interacts and its margin of stability is increased, typical of a flutter 

mechanism. The chart (see fig. 12) confirms that the process of implementing the RFA (see eq.

(8)), mass and stiffness matrices of the plant into a state space model was successful. The control 

characteristics of the wing model are presented in the next section.
22

Control Characteristics

Control design often takes place by shaping the loop gain and phase of a system. It is common 

to plot the magnitude and phase of the plant transfer function against frequency. The gain and 

phase plots of a system elicit properties of a system which are useful for control design. 

For the wing model, the gain and phase of two systems at different speeds are overlaid. The 

first system is modeled at an airspeed which is 0.5 m/s past the flutter speed. The second system 

is modeled at a much lower airspeed of 40 m/s where aerodynamic interaction is minimal. Both 

gain and phase of the plant model at the two different speeds are plotted versus frequency as

shown in figure 13.

Figure 13. Wing model plant mode for pre- and post-flutter: a) magnitude; and b) phase.

Figure 13(a) shows two plots. The plant at 40 m/s shows expected amplification at specific 

modal frequencies. The 1st bending and 2nd bending amplifications occur at their respective dry 

structural frequencies. The frequency at which amplification of the 1st torsion mode occurs is 

slightly shifted downward in frequency from 43 rad/s to 40 rad/s due to slight aerodynamic

coupling (see fig. 11). Figure 13(b) shows the separate lead and lag associated with the 

1

st bending and 1st torsion modes at lower airspeed. At the flutter speed, this lead and lag couple 

to produce an unfavorable coupling condition, or flutter. 

This frequency-based analysis agrees with the intuition of what might be expected from a 

model subject to flutter. A time history of the model from a control impulse gives important 

information about the system as well and is presented in figure 14.
23

Figure 14. Wing model leading edge wing tip accelerometer response to control impulse.

Figure 14 shows that the two system responses to the impulse were quite different. The model 

at lower airspeed experiences an initial perturbation and within 5 s it is suppressed. The model 

slightly above the flutter speed exhibits divergent oscillatory characteristics that one might expect 

from flutter. Its initial response was also much stronger, which is to be expected, as the damping 

is near 0 at the flutter airspeed. The time history indicates that after 1 s the flutter mode dominates 

the response and higher modes die out quickly. After this period of time, oscillations occur near 

the expected flutter frequency of 28 rad/s, as predicted from figure 11. The divergent 

characteristics of the time history make it apparent that control is required to suppress the 

instability. The next section shows the reaction of the model to the gust model which was given 

in figure 3.

Gust Response

For control design, a gust response can be useful to ascertain disturbance rejection in the 

system. The reaction of the state space model to the 1 − cos(x) gust model is shown in figure 15. 
24

Figure 15. Wing model leading edge wing tip accelerometer response to gust model.

The acceleration measurements at the leading-edge wing tip accelerometer move as expected

from a gust input. The average acceleration moves upward with the shape of the gust velocity 

input (see fig. 3) over the period of 3 s. The model at the flutter speed also exhibited less damping,

as expected. An interesting occurrence, however, was the sharp perturbation of the acceleration 

of the leading edge after the gust expired. This action is due to the wing moving back into position 

from its displaced position. The gust model is not validated with experimental data yet, but its 

performance appears adequate for disturbance rejection tests.

Conclusions and Future Work

An aeroservoelastic state space modeling tool for rectangular wing models was presented. 

The inputs to the tool are basic. The tool is highly parameterized, requiring very little design 

experience by the user; however, the outputs of the tool are of medium fidelity, giving the user a 

feel for realistic aeroservoelastic phenomena for complex structures. This tool allows for rapid 

investigation of aeroservoelastic models, improved understanding of the effects of changes, and 

most importantly, improved learning. It is hoped that this tool will become available to students 

and researchers who are interested in the field of aeroservoelasticity.

Verification and validation studies were presented for this tool. Toward this effort, unsteady 

aerodynamic code was validated against experimental data for a clamped plate. A generic wing 

model was also presented, its characteristics analyzed first with traditional flutter analyses. The 

traditional analyses were then compared to corresponding state space models and found to be in 

agreement. 

Future simulation work is planned to include incorporation of rigid body degrees of freedom. 

Future validation work is planned to include wing procurement based on the tool finite element 

model and further validation with wind-tunnel and ground vibration test measurements.
25

Acknowledgments

This research was funded by the Aeronautics Research Mission Directorate Subsonic Fixed 

Wing Project and was approved by Steve Jacobson, NASA Armstrong Controls Branch Chief.

Marty Brenner, NASA Armstrong researcher, has been a helpful advisor in the development of 

these results.

References

1. MSC Software Corporation. www.mscsoftware.com (cited 24 September 2014).

2. ZONA Technology Inc., ZAERO Theoretical Manual: Engineer’s Toolkit for Aeroelastic 

Solutions (ZAERO. Version 8.5), ZONA Technology, Inc., Scottsdale, Arizona, 2011. Available 

from publisher.

3. Krist, Sherrie L., Robert T. Biedron, and Christopher L. Rumsey, “CFL3D User's Manual 

(Version 5.0),” NASA/TM-1998-208444, 1998.

4. ESI Group. https://www.esi-group.com/software-services/virtual-environment/cfd-multiphysics/

ace-suite/cfd-fastran (accessed August 24, 2015).

5. Suh, Peter M., and Dimitri N. Mavris, “Modal Filtering for Control of Flexible Aircraft,” AIAA 

2013-1741, 2013.

6. Suh, Peter, “Robust Modal Filtering for Control of Flexible Aircraft,” Ph.D. Dissertation, School 

of Aerospace Engineering, Georgia Institute of Technology, Atlanta, Georgia, 2014.

7. Suh, Peter M., Alexander W. Chin, and Dimitri N. Mavris, “Virtual Deformation Control of the 

X-56A Model with Simulated Fiber Optic Sensors,” AIAA 2013-4844, 2013.

8. Szilard, Rudolph, Theories and Applications of Plate Analysis, John Wiley & Sons, Inc., New 

Jersey, 2004.

9. Albano, E., and Rodden, W. P., “A doublet lattice method for calculating lift distributions on 

oscillating surfaces in subsonic flows," AIAA Paper 68-73, 1968.

10. Blair, Max, “A Compilation of the Mathematics Leading to the Leading to the Doublet Lattice 

Method,” WL-TR-92-3028, 1992.

11. Katz, Joseph, and Allen Plotkin, Low-Speed Aerodynamics, Second Edition, Cambridge 

University Press, 2001.

12. Roger, Kenneth L., “Airplane Math Modeling Methods for Active Control Design,” AGARDCP-228, 1977.

13. Conyers, Howard J., Earl H. Dowell, and Kenneth C. Hall, “Aeroelastic Studies of a 

Rectangular Wing with a Hole: Correlation of Theory and Experiment,” 2010 Aerospace 

Systems Conference, National Society of Black Engineers, Los Angeles, California, 2010.

14. Nissim, E., and John J. Burken, “Control Surface Spanwise Placement in Active Flutter 

Suppression Systems,” NASA TP-2873, 1988.
26

15. Abel, Irving, “An Analytical Technique for Predicting the Characteristics of a Flexible Wing 

Equipped with an Active Flutter-Suppression System and Comparison with Wind-Tunnel 

Data," NASA TP-1367, 1979.

16. Gupta, K. K., M. J. Brenner, and L. S. Voelker, “Development of an Integrated 

Aeroservoelastic Analysis Program and Correlation with Test Data,” NASA TP-3120, 1991.

17. Schmidt, David K., Modern Flight Dynamics, McGraw-Hill, New York, 2012, pp. 129-154.

18. Gawronski, Wodek K., Advanced Structural Dynamics and Active Control of Structures, 

Springer-Verlag, New York, 2004.

19. Rodden, William P., Paul F. Taylor, and Samuel C. McIntosh Jr., “Further Refinement of the 

Subsonic Doublet-Lattice Method,” Journal of Aircraft, Vol. 35, No. 5, 1998, pp. 720-727.

20. Rodden, William P., Paul F. Taylor, and Samuel C. McIntosh Jr., “Improvements to the 

Doublet-Lattice Method in MSC/NASTRAN,” MSC 1999 Aerospace Users’ Conference 

Proceedings, 1999.

21. Giesing, J. P., T. P. Kalman, and W. P. Rodden, “Subsonic Unsteady Aerodynamics for

General Configurations,” AFFDL-TR-71-5, Part II, Volume II, 1972.

22. Karpel, Mordechay, “Design for Active Flutter Suppression and Gust Alleviation Using StateSpace Aeroelastic Modeling,” Journal of Aircraft, Vol. 19, No. 3, 1982, pp. 221-227.

23. Tiffany, Sherwood H., and Mordechay Karpel, “Aeroservoelastic Modeling and Applications 

Using Minimum-State Approximations of the Unsteady Aerodynamics,” NASA TM-101574, 

1989.

24. Pototzky, Anthony S., “Enhanced Modeling of First-Order Plant Equations of Motion for 

Aeroelastic and Aeroservoelastic Applications,” AIAA 2010-7801, 2010.

25. Adams Jr., William M., David M. Christhilf, Martin R. Waszak, Vivek Mukhopadhyay, and 

S. Srinathkumar, “Design, Test, and Evaluation of Three Active Flutter Suppression 

Controllers,” NASA TM-4338, 1992.

26. Wright, Jan R., and Jonathan E. Cooper, Introduction to Aircraft Aeroelasticity and Loads, 

John Wiley & Sons, Ltd, West Sussex, England, 2007.

27. ANSYS, Inc. http://www.ansys.com/ (accessed August 24, 2015).