NASA Technical Reports Server (NTRS) 20150000841: Rapid State Space Modeling Tool for Rectangular Wing Aeroservoelastic Studies

Rapid State Space Modeling Tool
for Rectangular Wing Aeroservoelastic Studies

Peter M. Suh 1

NASA Armstrong Flight Research Center, Edwards Air Force Base, California 93523

Howard J. Conyers 2

NASA Stennis Space Center, Stennis Space Center, Mississippi 39529
Dimitri N. Mavris 3

Georgia Institute of Technology, Atlanta, Georgia 30332

This paper 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 paper is on tool presentation, verification, and validation. These
processes are carried out in stages throughout the paper. 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 damping-versus-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.

A

AIC(lk )

A 0 ...A 6
B
C
C

Ci

c

Nomenclature

state space matrix

general aerodynamic influence coefficient matrix defined for reduced frequency, k

diagonal matrix of the area of all aerodynamic panels

rational function approximation matrices

control input matrix

global damping matrix

generalized structural and aerodynamic damping

[0 O 0 0 0 0]

mean aerodynamic chord

'Aerospace Engineer, PhD, Controls Branch/4840D, AIAA Member, peter.m.suh@ nasa.gov.

2 Aerospace Technologist, PhD, Engineering and Test Directorate, AIAA Member, howard.i.conyers@nasa.gov.

3 Boeing Professor of Advanced Aerospace Systems Analysis and Director, Aerospace Systems Design Laboratory
(ASDL), School of Aerospace Engineering, 270 Ferst Drive, Fellow AIAA, dimitri.mavris@ae.gatech.edu.

1

D(ik )

exp

F z

G

9

9 wash

H

i

K

K

k

L

l

M

M

My

m

n

n a

P(t)

Q(i/c)

Q

q _

9

q _s

s

s

T

t

u c

77

“-com

Ko

^C.p.

w c . p .(ik )

w g, o
Wg.max

Wg{t )

X

X I
*2

*3 -*6
Xc.p.

Xgust

X

y

Y

z f

z

z c.p.

Pi

aerodynamic influence matrix defined at control points of aerodynamic panels for k

exponential function

normal force applied to the beam

gust design constant

damping

gust wash

gust gradient

imaginary number

global stiffness matrix

generalized structural and aerodynamic stiffness
reduced frequency

total number of selected lag coefficients
lag constant iteration number
global mass matrix

generalized structural and aerodynamic mass

torque applied to beam about the length of the beam

number of modes

number of degrees of freedom

number of aerodynamic panels

external forcing function

generalized aerodynamic force matrix defined for a reduced frequency, k

rational function approximation of generalized aerodynamic forces

vector of modal coordinates

dynamic pressure

sensor measurement

defined as ik

Laplace variable

lag constant

time

control surface rotation
control command input
freestream velocity

matrix of w c p (tfc) collected for all m mode shapes

normalwash vector defined at the control points of the aerodynamic panels for reduced

frequency, k

gust velocity

design gust velocity

maximum gust velocity

gust velocity as a function of time

state vector

state vector of generalized modal displacements
state vector of generalized modal velocities
aerodynamic lag state vectors

vector of free-stream coordinates at the control points of each aerodynamic panel defined from
the leading edge of the wing

length from the leading edge of the wing to the gust origin
chordwise axis
measurement vector
spanwise axis

matrix of z c p collected for all m mode shapes
normal axis

out-of-plane z deformation of normal mode shape defined at control points of
aerodynamic panels
I th lag constant

2

f

0 c .p.

o

(j)

()

C)

(O r

C)

Qc

Oc#

O,

C0 a #

damping ratio

slope of the normal mode shape defined at the control points of the aerodynamic panels
modal matrix
circular frequency

l 54 derivative in time of the argument

2 nd derivative in time of the argument

transpose of the argument

generalized form of the argument

control partitioning of argument

control partitioning of argument for specific #

gust partitioning of argument

gust partitioning of argument for specific #

DLM

DOF

FEM

GAF

GVT

ODE

RFA

V-f

V-g

VLM

List of Acronyms

doublet-lattice method
degrees of freedom
finite element model
generalized aerodynamic force
ground vibration test

please define; used in discussion beginning at equation 10

rational function approximation

frequency versus velocity

damping versus velocity

vortex-lattice method

I. Introduction

A EROSERVOELASTIC modeling is often performed using medium-fidelity tools such as MSC Nastran 1 (MSC
Software Corporation, Santa Ana, California) and ZAERO 2 (ZONA Technology Inc., Scottsdale, Arizona). 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 3 (NASA Langley
Research Center, Flampton, Virginia) and CFD-FASTRAN 4 (ESI Group, Paris, France). 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. 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. 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. 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) 9 ’ 10 taking
advantage of recent research on improvements to this method. The steady forces were modeled with the
vortex-lattice method (VLM). 11 To transport the force coefficient information into a state space environment for
control design, a rational function approximation (RFA) method 12 was employed.

This paper 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. 13
The process of validation adds further confirmation to Conyers et al.’s original analyses and supports their work.

3

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.

II. 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 Fig. 1.

Geometry Parameters

Aspect ratio
Half span

Spar/rib/skin dimensions
Control surface
Number/size/layout
All materials

Discretization of model (x,y)

Aerodynamic Parameters

Reduced frequency range
Lag constants
Number of aerodynamic
lag states

Actuator transfer function

Geometry

module

Geometry

specifications

Wing skin, spars,
rib finite element
layout

FEM buildup

Control surface
finite element
layout /hinge
attachments

Control-fixed
global mass,
stiffness,
damping matrices^

Mass normalized
structural modes

Control-free
global mass,
stiffness.
damping matrices >

> Control modes

Aerodynamics
module

Normal wash/
gust modes

Automatic

aerodynamic

paneling

1 GAF matrix
per reduced
frequency

Rational
function
^ approximation ^

Analysis

module

V-g/V-f analysis

State space
modeling

Control design

V

J

Figure 1. Wing model design and aeroservoelastic analysis tool.

The three modules in Fig. 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) 14 15 and the corresponding RFAs. A direct output of the aerodynamics module
is the aerodynamic steady and unsteady matrices.

4

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.

A. 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 Ref. 16. A Lagrange representation 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 Eq. (1):

Mx + Cx + Kx + qAIC(lk)x — P{t) (1)

where M £ R nxn is the global mass matrix, C £ R nxn is the global damping matrix, K £ R nxn is the global
stiffness matrix, q is the dynamic pressure, AIC (ik) is the aerodynamic influence coefficient matrix computed at the
flight condition and for reduced frequency value k — c is the mean aerodynamic chord length, <>) is the circular

2Vqo

frequency, V m is the freestream velocity, x £ M nxl is 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 Eq. (2). 18

Mx + Kx = 0 (2)

The eigenvalue solution of the system produces the natural frequencies and eigenvectors (dry or in vacuo mode
shapes, <f>) of the system. The following linear transformation of state is formed as shown in Eq. (3):

x — <t>q (3)

where <f> £ M nxm is a matrix of eigenvectors corresponding to selected elastic and control modes, q. Pre-multiplying
Eq. (1) by O t and substituting Eq. (3), the matrix equation of motion is transformed as shown in Eq. (4).

0 T M0q + 0 T C0q + <t> T K®q + q® T AIC(lk)<t>q = 0 T P(t) (4)

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

Mq + Cq + Kq + qQ{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.

1. 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 19 is utilized in this tool. In addition, the
number of aerodynamic panels in the chordwise direction is constrained to be greater than ^^/ U q for generalized

5

force matrix computations. 2 " This is a suggested aerodynamic panel criterion from Rodden et al."° 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" 1 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 Eq. (6):

W c . p .(lfc) ®c.p. ~b t - Zc.p.

where 0 c p is the interpolated slope in the chord-wise direction of the aerodynamic panel defined at the control
points and z cp 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 Eq. (7) (see Ref. 15):

Q{lk)=Z f T Dm~ 1 A p W c , p , (7)

where Zf £ M naXm is a matrix of all normal deformations centered at the force point (quarter chord) on the
aerodynamic panels, W cp £ M n “ X7n is the matrix of normalwash vectors, for all aerodynamic mode shapes,
w c p (ik), D(ik ) £ M ,laXn “ is the aerodynamic influence matrix defined at control points of all panels for reduced
frequency k, and A p £ M n “ xna 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 Fig. 1.

With tabulated Q(tk ) computed, there are several ways of computing an RFA. Elere are two: Roger’s RFA 12 and
Karpel’s minimum-state method. 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 /?;, I — 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
Eq. (8):

L

Q(s) — A 0 + sAi + s 2 A 2 + ^ A 2+ i (8)

k s+ a

where s — s = 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. 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.

2. 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. 1 " 24 Therefore by inserting Eq. (8) into Eq. (5), the equations of motion are rewritten in Laplace form as
shown in Eq. (9):

s 2 Mq + sCq + Kq + q

SC

( sc \ 2

Ao + 2^ Al +

^2Y~) J ^ 2 ^3 X 3 + A 4 x 4 + •••

q — 0

(9)

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

6

sq

; 1 = 1 ...4

( 10 )

Equation (10) may be rewritten and inverse Laplace transformed resulting in matrix ODEs of the form given in
Eq. (11).

x 2+l + —^ PiX 2+ i = q\ 1 = 1... 4 (11)

c

Like terms are grouped in Eq. (9) and the equation may be rewritten after an inverse Laplace transform as in
Eq. (12).

(if + qA 0 )q + + q ^y~ A ij Q + yM + q A 2 j q + qA 3 x 3 + qA 4 x 4 + ••• - 0 (12)

Condensing variables in Eq. (12) results in the simplified matrix ODE given in Eq. (13).

Kq + Cq + Mq + qA 3 x 3 + qA 4 x 4 + ■•• = 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 Eq. (14) (see Ref. 24).

0

I

0

0

r xp

c

ql

ql)

'Xp

*2

*3

0

^3

2V m

-ft(— ) o

0

-

*2

x 3

0

0

Vx 6 )

0

^6

0 0

2Eoo

-*<T

0

0

0 -

0

0 -

-M-\K C

C c

M c )

( Uc ) 1

-M~\K g

Cg)

0

A c i

0

i£r*

0

A ai

0

^c4

0 .

0

A g*-

(14)

The vector states x x and x 2 replace the generalized modal state displacement and velocity vectors respectively. The
remaining states x 3 , ...x 6 are four orders of aerodynamic lag states. The control surface states u c ,u c ,u 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 w g and w g ,
respectively. The following section overviews how the sensors are modeled within the tool environment.

3. Sensor Outputs

Acceleration measurements at any place on the wing model are modeled assuming four aerodynamic lag states
as shown in Eq. (15):

7

<7s = t t >( 7 = [0 O 0 0 0 0]i = Cji

(15)

Pre-multiplying Eq. (14) by C 1 results in Eq. (16).

y = q s = C i* = C 1 Ax + C x Bu com (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 Eq. (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. 7 The next
section overviews the actuator modeling capabilities of the tool.

B. Actuators

For each control surface, a 2 nd order actuator is utilized. To model command and flight computer delay a 1 st
order lag filter is multiplied with the 2 nd 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 ld
order actuator model from command to surface movement is shown in Eq. (18):

u c 1 / a> 2

u com Ts + 1 \s 2 + 2 Olcfs + 0 ) 2

where T is a scalar time constant, o> is the circular frequency, <f is the damping ratio.

The actuator used in the results of this paper is given here. The time constant is set to .02 s; <u 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 Fig. 2.

Figure 2. Actuator model.

The command lag for this case removes much of the overshoot of the 2 nd order actuator alone, as shown in Fig 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.

C. 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. 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 Eq.
(19):

dwash GXp

(19)

where x cp is the vector of stream-wise coordinates at the control points of each aerodynamic panel and x gust is the
distance from the origin to the gust. If x gust 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 to be zero. A non-zero distance may result in large spirals
in the GAF matrices. 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 Eq. (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 Eq. (13). The inputs to the state space model are driven with a temporal 1 — cos gust model.
The standard temporal variation gust model" 6 used in the tool is given as in Eq. (20):

w g (t) =

( 20 )

where w g 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 Eq. (21):

w fl (t)=^p(l-cos(Gt}) (21)

where G is a gust design constant and w g max is the maximum gust velocity. A derivative in time of Eq. (21) leads to
the gust acceleration model shown in Eq. (22).

w g (t) = ^^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 Eq. (23).

£ _ 2 W g,max
w g,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
Fig. 3.

9

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

The gust model presented in Fig. 3 can be used for a realistic control design study. The input velocity and
acceleration directly influence the state derivatives in Eq. (14). The proper maximum amplitudes of the gust depend
on the theoretical operating design conditions. The particular gust model presented in Fig. 3 is used later in this
paper for wing model gust response testing. The next section introduces the verification and validation procedures
taken during tool development.

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

A. 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 {F z — — 100A1, M y = —100N — m} at beam tip. The force F z is applied in the normal direction, and M y is a
torque applied about the length of the beam at the tip. The dimensions of the beam are 1000mm x 500 mm x
6.35 mm, with a built-in (clamped) boundary condition at one end.

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 Fig. 4.

10

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 stiffness matrices and their verification results are not
presented here, but are used later on for the analysis of a wing model, where any modeling problems, if they did
exist, would become evident. The next section overviews the verification of the RFA used in the tool.

B. 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. 13 In Conyers et al.’s work the
clamped plate was analyzed by discretizing it into a 16 X 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. 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.

For the RFA analyses of the clamped plate the lag coefficients B 1 and B 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 X 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 (<f) as described in Eq. (8). The GAF coefficients corresponding to
the 1 st bending and 1 st torsion mode, and the corresponding RFA evaluations, are plotted together in Fig. 5.

11

Real part Real part

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

For both GAF coefficients presented in Fig. 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. 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.

C. Modal Analysis Validation Study

An analytical and experimental modal frequency analysis was conducted at Duke University by Conyers et al. 13
The simulation code modal frequency predictions are compared to that from both the ANSYS 27 (ANSYS, Inc.,
Canonsburg, Pennsylvania) 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.

Table 1. Modal frequency code comparisons and experimental results.

ANSYS 4
Frequencies, Hz

Tool FEM
Frequencies, Hz

Conyers et al. GVT, Hz 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

4 Republished data with permission of Dr. Floward Conyers. Experiments conducted at Duke University.

12

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.

D. Flutter Validation Study

Conyers et al. 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 3

Tool Flutter Code

Conyers et al. Wind
Tunnel Results 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 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 1 ' 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.

E. 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 Fig. 6.

13

Velocity, m/s

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

Figure 6 shows 10 modes, but only the 1 st bending and 1 st torsion mode interact strongly and lead to flutter at the
dashed vertical line. The flutter speed is computed where the 1 st 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 Fig. 7.

Velocity, m/s

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

14

Only two interacting modes are shown in Fig. 7. The first two modes at 0 airspeed begin at their dry natural
frequencies presented in Table 1. The 1 st torsion modal frequency begins to shift with airspeed, and at the flutter
speed of 19.9 m/s, the frequency of the 1 st torsion mode (the mode which crossed through 0 damping) is
10.9 rad/s. The frequency of the 1 st 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.

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

A. 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. 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
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 Fig. 8.

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

One difference between this model and the model presented previously^ 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

15

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.

B. Mode Shapes

The mode shapes of the wing model presented in Fig. 8 were computed with an eigenvalue analysis. The mode
shapes are made orthogonal to the mass matrix, which satisfies a critical mean axis constraint. 17 Ten mode shapes
were computed; the first four mode shapes are presented in Fig. 9.

a) First bending, frequency = 16 rad/s

b) First torsion, frequency = 43.4 rad/s

E

N

c) Second bending, frequency = 85.6 rad/s d) Second torsion, frequency = 184.6 rad/s

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

The mode shapes in Fig. 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 1 st bending mode is at 16 rad/s, and the 1 st torsion mode occurs at 43.4 rad/s. The 2 nd bending mode occurs
at 85.6 rad/s, and the 2 nd 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.

C. 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 Fig. 10.

16

First bending
First torsion

Stable

speeds

Unstable

speeds

— ■ Q — O

20

40

_l_

60

Velocity, m/s

80

Higher- /
order /
modes j

i

100

120

Figure 10. V-g analysis of wing model.

Figure 10 indicates the variation of 10 modes with velocity. Only the 1 st bending and 1 st torsion mode interact
strongly and lead to flutter at the dashed vertical line. As before, the flutter speed is computed where the 1 st 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 Fig. 6.
The major difference is that divergence is not reached in the given airspeed range. The V-f analysis is presented in
Fig. 11.

Figure 11. V-f analysis of wing model.

17

As before (see Fig. 7), only the two interacting modes are shown in Fig. 11. The first two modes start off at their
dry natural frequencies as shown in Figs. 9(a) and (b). The 1 st torsion modal frequency begins to shift with airspeed,
and at the flutter speed of 76.5 m/s the frequency of the 1 st 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.

D. State Space Model Analysis

The state space model was generated for the wing model as shown in Eq. (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.

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

Real part

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

2. 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 Fig. 13.

18

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 1 st bending and 2 nd 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). Figurel3(b) shows the separate lead and lag associated
with the 1 st bending and 1 st 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 Fig. 14.

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

19

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 Fig. 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 Fig 3.

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 Fig. 15.

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.

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

20

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.

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

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

'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, S. L„ Biedron, R. T, and Rumsey, C. L. “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 25, 2014).

5 Suh, P. M., and Mavris, D. N. “Modal Filtering for Control of Flexible Aircraft,” AIAA 2013-1741 , 2013.

6 Suh, P., “Robust Modal Filtering for Control of Flexible Aircraft,” Ph.D. Dissertation, School of Aerospace Engineering,
Georgia Institute of Technology, Atlanta, Georgia, 2014.

7 Suh, P. M., Chin, A. W., and Mavris, N. D., “Virtual Defonnation Control of the X-56A Model with Simulated Fiber Optic
Sensors,” AIAA 2013-4844, 2013.

8 Szilard, R., 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, M., “A Compilation of the Mathematics Leading to the Leading to the Doublet Lattice Method,” WL-TR-92-3028,
1992.

u Katz, J., and Plotkin, A., Low-Speed Aerodynamics . Second Edition, Cambridge University Press, 2001.

12 Roger, K. L., “Airplane Math Modeling Methods for Active Control Design,” AGARD-CP-228, 1977.

13 Conyers, H. J., Dowell, E. H., and Hall, K. C., “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 Burken, J. J., “Control Surface Spanwise Placement in Active Flutter Suppression Systems,”
NASA TP-2873, 1988.

15 Abel, I., “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., Brenner, M. J., and Voelker, L. S., “Development of an Integrated Aeroservoelastic Analysis Program and
Correlation with Test Data,” NASA TP-3120, 1991.

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

ls Gawronski, W. K., Advanced Structural Dynamics and Active Control of Structures, Springer- Verlag, New York, 2004.

19 Rodden, W. P., Taylor, P. F., and McIntosh Jr., S. C., “Further Refinement of the Subsonic Doublet-Lattice Method,”
Journal of Aircraft, Vol. 35, No. 5, 1998, pp. 720-727.

20 Rodden, W. P„ Taylor, P. F., and McIntosh Jr., S. C. "’’Improvements to the Doublet-Lattice Method in MSC/NASTRAN,”
MSC 1999 Aerospace Users’ Conference Proceedings, 1999.

21 Giesing, J. P., Kalman, T. P„ and Rodden, W. P., “Subsonic Unsteady Aerodynamics for General Configurations,”
AFFDL-TR-71-5, Part II, Volume II, 1972.

' 2 Karpel, M., “Design for Active Flutter Suppression and Gust Alleviation Using State-Space Aeroelastic Modeling,” Journal
of Aircraft, Vol. 19, No. 3, 1982, pp. 221-227.

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

24 Pototzky, A. S., “Enhanced Modeling of First-Order Plant Equations of Motion for Aeroelastic and Aeroservoelastic
Applications,” AIAA 2010-7801, 2010.

25 Adams Jr., W. M., Christhilf, D. M., Waszak, M. R„ Mukhopadhyay, V., and Srinathkumar, S. “Design, Test, and
Evaluation of Three Active Flutter Suppression Controllers,” NASA TM-4338, 1992.

26 Wright, J. R., and Cooper, J. E., Introduction to Aircraft Aeroelasticity and Loads. John Wiley & Sons, Ltd, West Sussex,
England, 2007.

27 ANSYS, Inc. http://www.ansys.com/ (cited 24 September 2014).

21