Aeroelastic System Development Using Proper
Orthogonal Decomposition and Volterra
Theory
David J. Lucia
∗
and Philip S. Beran
†
Air Force Research Laboratory
Walter A. Silva
‡
NASA Langley Research Center
This research combines Volterra theory and proper orthogonal decomposition (POD)
into a hybrid methodology for reduced-order modeling of aeroelastic systems. The out-
come of the method is a set of linear ordinary differential equations (ODEs) describing
the modal amplitudes associated with both the structural modes and the POD basis
functions for the fluid. For this research, the structural modes are sine waves of varying
frequency, and the Volterra-POD approach is applied to the fluid dynamics equations.
The structural modes are treated as forcing terms which are impulsed as part of the
fluid model realization. Using this approach, structural and fluid operators are coupled
into a single aeroelastic operator. This coupling converts a free boundary fluid problem
into an initial value problem, while preserving the parameter (or parameters) of interest
for sensitivity analysis. The approach is applied to an elastic panel in supersonic cross
flow. The hybrid Volterra-POD approach provides a low-order fluid model in state-space
form. The linear fluid model is tightly coupled with a nonlinear panel model using an
implicit integration scheme. The resulting aeroelastic model provides correct limit-cycle
oscillation prediction over a wide range of panel dynamic pressure values. Time inte-
gration of the reduced-order aeroelastic model is four orders of magnitude faster than
the high-order solution procedure developed for this research using traditional fluid and
structural solvers.
Introduction
Volterra methods
1
and proper orthogonal decompo-
sition
2, 3
(POD) are two of the more prevalent reduced-
order modeling (ROM) techniques well-suited to non-
linear dynamics.
4–8
The application of ROM tech-
niques to aeroelastic systems is an active area of re-
search, motivated by the desire for faster algorithms
that are well-suited to the design environment for air-
craft. For example, transonic, fluid-structure inter-
action is a particular application of interest to both
external and internal aerodynamicists because moving
shock waves in the flow necessitate high-fidelity numer-
ical flow solvers which are too cumbersome for iterative
design analysis. Regardless of the application, when
nonlinearities are present in either the flow field or the
∗
Senior Research Aerospace Engineer, Major, USAF,
AFRL/VAS, Bldg 45, 2130 Eighth Street, Suite 1, WPAFB,
OH 45433-7542, (david.lucia@wpafb.af.mil), AIAA Member
†
Principal Research Aerospace Engineer, AFRL/VASD,
Bldg 146, 2210 Eighth Street, WPAFB, OH 45433-7531,
(philip.beran@wpafb.af.mil), Associate Fellow of AIAA
‡
Senior Research Scientist, Aeroelasticity Branch, Mail Stop
340, NASA Langley Research Center, Hampton, VA 23681-0001,
Senior Member of AIAA
Copyright
c
° 2003 by the American Institute of Aeronautics and
Astronautics, Inc. No copyright is asserted in the United States
under Title 17, U.S. Code. The U.S. Government has a royalty-
free license to exercise all rights under the copyright claimed herein
for Governmental Purposes. All other rights are reserved by the
copyright owner.
structure, established order-reduction methods that
rely on linearized dynamics are of little use.
Over the past three years, applications of POD
to the Euler equations have produced reduced order
aeroelastic models that properly capture aerodynamic
nonlinearities. A low-order POD representation of
the discrete, 2-D Euler equations
9
was coupled with
the von K´arm´an equation to simulate the dynamics
of flow over a flexible panel.
10
Subsequently, a new
approach was taken, involving domain decomposition,
that allowed LCO to be accurately simulated in the
transonic regime.
11
In that study, full-order and re-
duced order models of a small flow region containing
a moving shock were decomposed from the larger flow
domain. Both approaches enabled a physically consis-
tent treatment of the aerodynamic nonlinearity. In a
more recent paper,
12
the original POD/ROM method-
ology used for flow over an elastic panel
10
was revis-
ited to improve the temporal coupling between the
aerodynamic and structural dynamic equations. Fur-
thermore, a modal representation of the structure was
employed, which permitted a more efficient formula-
tion of the reduced-order aeroelastic system.
All of the studies mentioned above relied on a ROM
technique called subspace projection for time integra-
tion of the reduced-order model. While sufficient to
demonstrate the accuracy of the POD basis functions,
1 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
subspace projection was not an efficient way to time
integrate the low-order, aeroelastic ROM. Generally,
four orders of magnitude reduction in fluid system de-
grees of freedom (DOFs) were demonstrated in the
above studies. Time integrating these POD/ROMs
with subspace projection generally produced about
one order of magnitude improvement in compute time
to accompany a much larger drop in problem order.
The applicability of POD basis functions to nonlin-
ear problems has been documented in the literature,
but a tractable nonlinear, low-order model realiza-
tion procedure is a key missing link. Two techniques,
Galerkin projection and direct projection, have been
recently reported as having potential for obtaining
nonlinear terms for POD/ROMs.
13
However, the lin-
ear portion of these realization procedures is generally
unstable, requiring dissipation techniques that affect
model performance. The Volterra-POD approach pro-
vides a stable reduced-order equation set, and is an
important advance toward achieving stable, nonlinear
reduced-order models.
The hybrid Volterra-POD method was recently de-
veloped to replace subspace projection for time inte-
gration of POD/ROMs applied to compressible flow
fields.
14
The goal of the Volterra-POD approach was
to achieve computational savings on the order of DOF
reductions. This goal was achieved in the initial ap-
plication, where four orders of magnitude reduction
was obtained in both DOFs and compute time. To
date, the hybrid Volterra-POD method has only been
applied to subsonic flow-fields characterized by linear
behavior, with fixed boundaries. The product of the
technique was a linear, state-space system of ODEs
governing the dynamics of modal coefficients corre-
sponding to a small number of POD basis functions.
The state-space realization was obtained from a set of
impulse responses that were processed using the Eigen-
system Realization Algorithm (ERA).
15, 16
This research will extend the Volterra-POD ap-
proach to supersonic flow-fields with dynamic bound-
ary behavior. The POD-Volterra method will be ap-
plied to a two-dimensional elastic panel in inviscid,
supersonic cross-flow. The Volterra-POD approach
will be used to identify a low-dimensional, linear
POD/ROM for the fluid. The POD/ROM will be
tightly coupled to a low-dimensional, nonlinear model
of the von K´arm´an plate equation.
17
The aeroelas-
tic response will be obtained using an implicit time-
integration scheme.
The Volterra-POD technique involves procedures
that require the selection of parameters such as im-
pulse size, several data windowing lengths, and im-
pulse sampling frequency. The choice of POD basis
affects performance as well. Some considerations for
generating the POD basis include choice of base flow
(the POD/ROM determines the perturbation to this
base flow), snapshot collection method and sampling
frequency (the method of snapshots
18
will be discussed
in the next section).
The research will consider two base flow cases, and
two snapshot collection methods. Both uniform flow
at free stream conditions, and steady-state flow over
a static panel deflection will be considered as base
flow cases. An aeroelastic POD basis will be gener-
ated by sampling a small portion of the time history
for a baseline LCO case, which was the approach in
recent applications using subspace projection for this
problem.
12
In addition, we will investigate using the
impulse response of the fluid system to generate a POD
basis. The full-system impulse response is collected as
part of the Volterra-POD approach, and the impulse
responses can be sampled as snapshots in lieu of the
LCO time-history. Finally, we will apply POD to the
structural dynamics, couple the structural POD/ROM
with the fluid POD/ROM and examine performance.
We will record the various parameter settings used to
generate the aeroelastic POD/ROMs for each case.
The linearity of the supersonic flow-field will be ex-
amined as part of the ROM analysis. The principle
of superposition applies in a linear flow-field, which
enables a host of linear order-reduction techniques, in-
cluding the Volterra-POD technique detailed in this
paper. While the supersonic, aeroelastic flow-field is
well represented by a linear fluid model, we will demon-
strate that the supersonic flow-field itself is not linear
in general.
The performance of the Volterra-POD aeroelastic
ROMs will be quantified in accuracy, order reduction,
and computational savings. A high-order, full-system
representation of the problem is required for snapshot
collection. The flow field and panel response for the
full-system model will serve as the baseline for per-
formance comparison. Accuracy will be quantified by
comparing LCO panel response, flow-field pressure dis-
tribution on the elastic panel, LCO frequency, and
LCO phase for a variety of panel dynamic pressure
values. Finally, the robustness of the Volterra-POD
method for predicting LCO response across a broad
parameter space will also be addressed.
Formulation
This section describes the full-system aeroelastic
model, introduces POD, and overviews Volterra meth-
ods. In addition, we fully develop the Volterra-POD
approach and the synthesis of aeroelastic ROMs.
Structural Dynamics Equations
Two-dimensional flow over a semi-infinite, pinned
panel of length L is considered. Panel dynamics are
computed with von K´arm´an’s large-deflection plate
equation, which is placed in nondimensional form us-
ing aerodynamic scales L and u
∞
19
(0 < x < 1):
µ
λ
∂
4
w
∂x
4
− N
x
∂
2
w
∂x
2
+
∂
2
w
∂t
2
= µ
µ
1
γM
2
∞
− p
¶
, (1)
2 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
N
x
≡
6µ
λ
µ
h
L
¶
−2
¡
1 − ν
2
¢
Z
1
0
µ
∂w
∂x
¶
2
dx. (2)
The nonlinear, in-plane load in Eqn. (2), serves to limit
panel deflections w(x, t) induced by fluid-structure in-
teraction. Here, the load is assumed to be distributed
uniformally over the panel.
17
Equation (1) is compa-
rable to similar formulations found in the literature,
17
although w
d
and t
d
are scaled by h and
¡
ρ
s
hL
4
¢
1/2
,
respectively. Two pinned boundary conditions are en-
forced at the panel’s endpoints: w = 0 and
∂
2
w
∂x
2
= 0.
A modal solution for the deflection w(x, t) is as-
sumed:
w(x, t) =
m
s
X
i=1
a
i
(t) sin(iπx), (3)
where m
s
is the number of structural modes retained,
and the modal amplitudes a
i
vary in time and are col-
located in the array a. The Galerkin method is used
to obtain a low-order set of ordinary-differential equa-
tions describing the behavior of a
i
.
17
First, Eqn. (3)
is substituted into Eqn. (1). The resulting expres-
sion is then integrated, following pre-multiplication by
sin(iπx), to yield (i = 1, ..., m
s
)
1
2
¨a
i
+
µ(iπ)
4
2λ
a
i
+
6µ
λ
µ
h
L
¶
−2
¡
1 − ν
2
¢
α
(iπ)
2
2
a
i
= µP
i
,
(4)
where α ≡
P
r
a
2
r
(rπ)
2
2
and
P
i
≡
Z
1
0
µ
1
γM
2
∞
− p
¶
sin(iπx)dx. (5)
The projected pressure components, P
i
, are integrated
from the aerodynamic solution with the midpoint rule,
using flowfield pressures obtained at grid points on the
panel surface.
12
The aerodynamic equations, their dis-
cretization, and their solution are discussed in later
sections. While equivalent to other formulations in
the literature, Eq. (4) has two notable distinctions.
First, the different form of scaling described above al-
ters equation coefficients, and, second, an expression
relating p to the state of the panel is not assumed.
12
The structural dynamics equation Eqn. (4) is placed
in first-order form by introducing a mode speed array,
b, such that ˙a
i
= b
i
,
˙
b
i
= −
"
µ(iπ)
4
λ
+
6µ
λ
µ
Liπ
h
¶
2
¡
1 − ν
2
¢
α
#
a
i
+ 2µP
i
.
(6)
The mode speeds and amplitudes are collocated into
a structural solution array, Y
s
, leading to a general
form of the structural equation:
Y
s
= [b, a]
T
(7a)
˙
Y
s
= R
s
(Y
s
, P ; µ, λ, h/L). (7b)
Fluid Dynamics
The dynamics of inviscid fluid flows are governed by
the Euler equations. The two-dimensional Euler equa-
tions are given below in strong conservation form:
20
∂U
∂t
+
∂E(U )
∂x
+
∂F (U )
∂y
= 0 , (8a)
U(x, t) =
ρ
m
x
m
y
E
T
, (8b)
where ρ, m
x
, m
y
and E
T
are functions of space and
time. Since we assume an ideal gas for our applica-
tions, this equation set can be closed using the ideal
gas law.
The solution of the Euler equations can be approx-
imated using either finite-difference, finite-volume, or
finite-element techniques. To do this, the spatial do-
main is discretized, and the flow variables in U (x, t) at
each discrete location are collocated into a column vec-
tor U(t). Time integration across the computational
mesh is used to obtain flow solutions.
Since the Euler equations are linear in the time
derivative, and quasi-linear in the spatial deriva-
tive,
20, 21
the spatial derivatives and the time deriva-
tives in Eqn. (8a) can be separated to form an evolu-
tionary system. To accomplish this, the spatial deriva-
tives of the flux terms
∂E
∂x
and
∂F
∂y
are grouped to form
a nonlinear operator R acting on the set of fluid vari-
ables. The fluid dynamics from Eqn. (8a) can then be
expressed as
dU(x, t)
dt
= R(U(x, t)) . (9)
When discretized this expression takes the form
dU (t)
dt
= R(U (t)) . (10)
Equation (10) is referred to as the full-system dynam-
ics.
A finite-volume scheme was the basis for the full-
order solver used in this research, which approximated
the integral form of the Euler equations:
d
dt
Z
V
UdV +
Z
∂V
(Ebı + F b) · d
S = 0 . (11)
The grid points in the computational mesh described
earlier were used to form corners for cells. For each
cell, the integral form of the Euler equations reduced
to the following, assuming no grid deformation:
d
dt
U
i,j
+
X
sides
(E
i,j
bı + F
i,j
b) ·
d
S
i,j
dA
i,j
= 0 . (12)
The flux terms E
i,j
and F
i,j
were computed using
second-order Roe averaging,
20
and the flow variables
3 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
U
i,j
were evaluated as cell averages. Time integration
across the computational mesh was used to obtain flow
solutions. This was accomplished with a first-order-
accurate, forward Euler approximation.
External boundaries were handled with ghost cells.
The fluid values for the ghost cells at the far
field boundaries were determined using characteristic
boundary conditions.
20
The bump-surface was mod-
eled using a transpiration approximation.
22
The finite-
volume fluid solver and the transpiration boundary
condition were validated using a combination of the-
ory and experimental data. Subsonic performance was
validated using wind-tunnel data.
23
Supersonic per-
formance was validated using oblique shock theory.
24
Time-accurate performance was validated by correctly
predicting the vortex shedding frequency from a cylin-
der in cross flow.
Time Integration of the Coupled Full-Order
Equations
The systems of discretized fluid dynamic equations,
U (t), and modal structural equations, Y
s
, are com-
bined into a single time-dependent system represen-
tative of the complete interaction between structure
and inviscid flow. Time integration proceeds in two
steps, assuming an O(∆t) lag in the synchronization
of fluid and structure. First, the structural variables
are updated from time level n to n + 1 using a Crank-
Nicolson procedure to be described below (but limited
here to only structural variables). During this step, the
pressures known at grid points on the panel surface are
considered frozen. In the second step, the aerodynamic
variables are explicitly updated using only structural
variables defined at time level n.
Grid Generation and Time Step
The flow is simulated over a physical domain of
length D
L
, centered about x = 0, and height D
H
. The
domain is discretized using I nodes in the streamwise
direction and J nodes normal to the panel. Indices i
(1 ≤ i ≤ I) and j (1 ≤ j ≤ J) are used to denote grid
points comprising corners of cells for the finite-volume
scheme. Grid points are clustered in the direction nor-
mal to the panel at the panel surface, with minimum
spacing denoted by ∆
wall
. The spacing of grid points is
specified to grow geometrically with j from the panel
boundary. In the streamwise direction, the node spac-
ing is chosen to be uniform over the deforming panel
segment (coincident with the structural grid), while
growing geometrically upstream of the leading edge
(positioned at i = I
LE
) and downstream of the trail-
ing edge (positioned at i = I
T E
). Calculations are
carried out with a baseline grid given by the following:
I = 141, J = 116, D
L
= 50, D
H
= 25, I
LE
= 45,
I
T E
= 97, and ∆
wall
= 0.0125.
Proper Orthogonal Decomposition
POD is a technique to identify a small number of
basis functions that adequately describe the behav-
ior of the full-system dynamics (Eqn. (10)) across
some parameter space of interest. A summary of POD
as it applies to a spatially-discretized flow field fol-
lows. A detailed description of POD is available in the
literature.
4, 8
For simplicity, consider only one fluid
variable, w(x, t), which when spatially descritized us-
ing N nodes is denoted w(t). For this fluid variable,
the full-system dynamics in Eqn. (10) is expressed as
dw
dt
= R
w
(w) . (13)
Spectral methods approximate the solution w(X, t) as
w(x, t) ≈
M
X
k=1
a
k
(t)φ
k
(x) . (14)
When the domain is spatially discretized, φ
k
(x) be-
comes a vector φ
k
, and the following relation applies:
w(t) ≈
M
X
k=1
a
k
(t)φ
k
. (15)
The set of vectors {φ
k
} are discrete basis functions
corresponding to the computational mesh defined for
the numerical solver. The set {a
k
} are the modal coef-
ficients, and Eqn. (15) can be represented using matrix
algebra. The fluid modes comprise columns of a modal
matrix Φ, and the coefficients are collocated into a
column vector
b
w(t). POD produces a linear transfor-
mation Φ between the full-order solution, w, and the
reduced-order solution,
b
w:
w(t) = W
0
+ Φ
b
w(t) . (16)
The reduced-order variable
b
w(t) represents deviations
of w(t) from a base solution W
0
. The subtraction of
W
0
will result in zero-valued boundaries for the POD
modes wherever constant boundary conditions occur
on the domain.
Φ is constructed by collecting observations of the
solution w(t)−W
0
at different time intervals through-
out the time integration of the full-system dynamics.
These observations are called snapshots
18
and are gen-
erally collected to provide a good variety of flow field
dynamics while minimizing linear dependence. The
snapshot generation procedure is sometimes referred
to as POD training.
8
A total of Q snapshots are collected from the full-
system dynamics. These are vectors of length N. The
set of snapshots describe a linear space that is used
to approximate both the domain and the range of the
nonlinear operator R
w
. The linear space is defined
by the span of the snapshots.
25
POD identifies a new
basis for this linear space that is optimally convergent
4
4 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
in the sense that no other set of basis functions will
capture as much energy in as few dimensions as the
POD basis functions. To identify the POD basis, the
snapshots are compiled into an N ×Q matrix S, known
as the snapshot matrix. The mapping function Φ is
then developed using
S
T
SV = V Λ , (17a)
Φ = SV . (17b)
Here V is the matrix of eigenvectors of S
T
S, and Λ is
the corresponding diagonal matrix of eigenvalues. To
eliminate redundancy in the snapshots, the columns
of V corresponding to very small eigenvalues in Λ are
truncated. The matrix of eigenvalues Λ is also resized
to eliminate the rows and columns corresponding to
the removed eigenvalues. If Q − M columns of V are
truncated, the resulting reduced order mapping Φ will
be an N × M matrix. The reduced-order mapping
represents an isomorphism between N and M dimen-
sional space. It determines the coordinates of w(t) in
terms of the M remaining basis functions, φ
k
.
The reduced-order mappings for each fluid variable
are developed separately, and individual S and V ar-
rays are collocated as blocks into a larger set of arrays,
also denoted S and V , to form
U (t) ≈ U
0
+ Φ
b
U (t) , (18a)
Φ = SV . (18b)
These versions of Eqn.’s (16) and (17b), respectively,
apply to the entire set of fluid variables.
POD of the discrete, panel position vector w(x, t) →
w(t) and panel velocity vector s(t) =
˙
w(t) is accom-
plished in a similar manner as described above for the
fluid system. Unlike the fluid POD basis functions,
there is no base term subtracted from the snapshots
when generating a structural POD basis.
Volterra Methods
Consider time-invariant, nonlinear, continuous-
time, systems. Of interest is the response of the system
about an initial state X(0) = X
0
due to an arbi-
trary input u(t) (we take u as a real, scalar input)
for t ≥ 0. As applied to these systems, Volterra the-
ory
1, 26–28
yields the response
X(t) = h
0
+
Z
t
0
h
1
(t − τ )u(τ)dτ (19)
+
Z
t
0
Z
t
0
h
2
(t − τ
1
, t − τ
2
)u(τ
1
)u(τ
2
)dτ
1
dτ
2
+
N
X
n=3
Z
t
0
..
Z
t
0
h
n
(t − τ
1
, .., t − τ
n
)
u(τ
1
).. u(τ
n
)dτ
1
.. dτ
n
.
The Volterra series can be accurately truncated be-
yond the second-order term when a weakly nonlinear
formulation is considered:
X(t) = h
0
+
Z
t
0
h
1
(t − τ )u(τ)dτ
+
Z
t
0
Z
t
0
h
2
(t − τ
1
, t − τ
2
)u(τ
1
)u(τ
2
)dτ
1
dτ
2
. (20)
The assumption of a weakly nonlinear system is con-
sistent with the emergence of limit-cycle oscillation of
a 2-D aeroelastic system in transonic flow through a
supercritical Hopf bifurcation.
29
For linear systems,
only the first-order kernel is non-trivial, and there are
no limitations on input amplitude.
The first- and second-order kernels are presented be-
low in final form:
5
h
1
(τ
1
) = 2X
0
(τ
1
) −
1
2
X
2
(τ
1
) , (21)
h
2
(τ
1
, τ
2
) =
1
2
(X
1
(τ
1
, τ
2
) − X
0
(τ
1
) − X
0
(τ
2
)) .(22)
In (21), X
0
(τ
1
) is the time response of the system to
a unit impulse applied at time 0 and X
2
(τ
1
) is the
time response of the system to an impulse of twice
unit magnitude at time 0. These response functions
represent the memory of the system. If the system is
linear, then X
2
= 2X
0
and h
1
= X
0
, which is why
the first-order kernel is referred to as the linear unit
impulse response. The identification of the second-
order kernel is more demanding, since it is dependent
on two parameters. Assuming τ
2
> τ
1
in (22), X
0
(τ
2
)
is the response of the system to an impulse at time τ
2
.
Time is discretized with a set of time steps of equiv-
alent size. Time levels are indexed from 0 (time 0) to
n (time t), and the evaluation of X at time level n is
denoted by X[n]. The convolution in discrete time is
X[n] = h
0
+
N
X
k=0
h
1
[n − k]u[k] (23)
+
N
X
k
1
=0
N
X
k
2
=0
h
2
[n − k
1
, n − k
2
]u[k
1
]u[k
2
] .(24)
The linearized and nonlinear Volterra kernels can
be transformed into linearized and nonlinear (bilinear)
state-space systems that can be easily implemented
into other disciplines such as controls and optimiza-
tion.
5, 27
For linear dynamics, state-space realization
using the Eigensystem Realization Algorithm (ERA)
has been used to generate linear, aeroelastic systems.
16
Nonlinear system realization is an active area of re-
search.
System Realization
The ERA method
15
identifies a discrete, linear,
time-invariant state-space realization of the form,
X[n + 1] = AX[n] + Bu[n]
Y [n] = CX[n] , (25)
5 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
