Table Of ContentAeroelastic System Development Using Proper
Orthogonal Decomposition and Volterra
Theory
David J. Lucia⁄ and Philip S. Berany
Air Force Research Laboratory
Walter A. Silvaz
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 difierential equations (ODEs) describing
the modal amplitudes associated with both the structural modes and the POD basis
functions for the (cid:176)uid. For this research, the structural modes are sine waves of varying
frequency, and the Volterra-POD approach is applied to the (cid:176)uid dynamics equations.
The structural modes are treated as forcing terms which are impulsed as part of the
(cid:176)uid model realization. Using this approach, structural and (cid:176)uid operators are coupled
into a single aeroelastic operator. This coupling converts a free boundary (cid:176)uid 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
(cid:176)ow. ThehybridVolterra-PODapproach providesa low-order (cid:176)uidmodel in state-space
form. The linear (cid:176)uid 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 (cid:176)uid and
structural solvers.
Introduction structure, established order-reduction methods that
Volterramethods1 andproperorthogonaldecompo- rely on linearized dynamics are of little use.
sition2,3(POD)aretwoofthemoreprevalentreduced- Over the past three years, applications of POD
order modeling (ROM) techniques well-suited to non- to the Euler equations have produced reduced order
linear dynamics.4{8 The application of ROM tech- aeroelastic models that properly capture aerodynamic
niques to aeroelastic systems is an active area of re- nonlinearities. A low-order POD representation of
search, motivated by the desire for faster algorithms the discrete, 2-D Euler equations9 was coupled with
that are well-suited to the design environment for air- the von K¶arm¶an equation to simulate the dynamics
craft. For example, transonic, (cid:176)uid-structure inter- of (cid:176)ow over a (cid:176)exible panel.10 Subsequently, a new
action is a particular application of interest to both approach was taken, involving domain decomposition,
externalandinternalaerodynamicistsbecausemoving that allowed LCO to be accurately simulated in the
shockwavesinthe(cid:176)ownecessitatehigh-fldelitynumer- transonic regime.11 In that study, full-order and re-
ical(cid:176)owsolverswhicharetoocumbersomeforiterative duced order models of a small (cid:176)ow region containing
design analysis. Regardless of the application, when a moving shock were decomposed from the larger (cid:176)ow
nonlinearitiesarepresentineitherthe(cid:176)owfleldorthe domain. Both approaches enabled a physically consis-
tent treatment of the aerodynamic nonlinearity. In a
⁄Senior Research Aerospace Engineer, Major, USAF, morerecentpaper,12 theoriginalPOD/ROMmethod-
AFRL/VAS, Bldg 45, 2130 Eighth Street, Suite 1, WPAFB, ology used for (cid:176)ow over an elastic panel10 was revis-
OH45433-7542,([email protected]),AIAAMember
yPrincipal Research Aerospace Engineer, AFRL/VASD, ited to improve the temporal coupling between the
Bldg 146, 2210 Eighth Street, WPAFB, OH 45433-7531, aerodynamic and structural dynamic equations. Fur-
([email protected]),AssociateFellowofAIAA thermore, a modal representation of the structure was
zSeniorResearchScientist,AeroelasticityBranch,MailStop
employed, which permitted a more e–cient formula-
340,NASALangleyResearchCenter,Hampton,VA23681-0001,
tion of the reduced-order aeroelastic system.
SeniorMemberofAIAA
Copyright(cid:176)c 2003bytheAmericanInstituteofAeronauticsand AllofthestudiesmentionedabovereliedonaROM
Astronautics, Inc. No copyright is asserted in the United States
under Title 17, U.S. Code. The U.S. Government has a royalty- technique called subspace projection for time integra-
freelicensetoexerciseallrightsunderthecopyrightclaimedherein tion of the reduced-order model. While su–cient to
for Governmental Purposes. All other rights are reserved by the
copyrightowner. demonstrate theaccuracyofthePODbasisfunctions,
1 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
subspace projection was not an e–cient way to time frequency(themethodofsnapshots18willbediscussed
integrate the low-order, aeroelastic ROM. Generally, in the next section).
four orders of magnitude reduction in (cid:176)uid system de- The research will consider two base (cid:176)ow cases, and
grees of freedom (DOFs) were demonstrated in the two snapshot collection methods. Both uniform (cid:176)ow
above studies. Time integrating these POD/ROMs at free stream conditions, and steady-state (cid:176)ow over
with subspace projection generally produced about a static panel de(cid:176)ection will be considered as base
one order of magnitude improvement in compute time (cid:176)ow cases. An aeroelastic POD basis will be gener-
to accompany a much larger drop in problem order. ated by sampling a small portion of the time history
The applicability of POD basis functions to nonlin- for a baseline LCO case, which was the approach in
ear problems has been documented in the literature, recent applications using subspace projection for this
but a tractable nonlinear, low-order model realiza- problem.12 In addition, we will investigate using the
tion procedure is a key missing link. Two techniques, impulseresponseofthe(cid:176)uidsystemtogenerateaPOD
Galerkin projection and direct projection, have been basis. The full-system impulse responseis collected as
recently reported as having potential for obtaining part of the Volterra-POD approach, and the impulse
nonlinear terms for POD/ROMs.13 However, the lin- responses can be sampled as snapshots in lieu of the
earportionoftheserealizationproceduresisgenerally LCO time-history. Finally, we will apply POD to the
unstable, requiring dissipation techniques that afiect structuraldynamics,couplethestructuralPOD/ROM
modelperformance. TheVolterra-PODapproachpro- with the (cid:176)uid POD/ROM and examine performance.
vides a stable reduced-order equation set, and is an We will record the various parameter settings used to
important advance toward achieving stable, nonlinear generate the aeroelastic POD/ROMs for each case.
reduced-order models. The linearity of the supersonic (cid:176)ow-fleld will be ex-
The hybrid Volterra-POD method was recently de- amined as part of the ROM analysis. The principle
veloped to replace subspace projection for time inte- of superposition applies in a linear (cid:176)ow-fleld, which
gration of POD/ROMs applied to compressible (cid:176)ow enablesahostoflinearorder-reductiontechniques,in-
flelds.14 The goal of the Volterra-POD approach was cluding the Volterra-POD technique detailed in this
toachievecomputationalsavingsontheorderofDOF paper. While the supersonic, aeroelastic (cid:176)ow-fleld is
reductions. This goal was achieved in the initial ap- wellrepresentedbyalinear(cid:176)uidmodel,wewilldemon-
plication, where four orders of magnitude reduction strate that the supersonic (cid:176)ow-fleld itself is not linear
was obtained in both DOFs and compute time. To in general.
date, the hybrid Volterra-POD method has only been The performance of the Volterra-POD aeroelastic
applied to subsonic (cid:176)ow-flelds characterized by linear ROMs will be quantifled in accuracy, order reduction,
behavior, with flxed boundaries. The product of the and computational savings. A high-order, full-system
technique was a linear, state-space system of ODEs representation of the problem is required for snapshot
governing the dynamics of modal coe–cients corre- collection. The (cid:176)ow fleld and panel response for the
sponding to a small number of POD basis functions. full-system model will serve as the baseline for per-
The state-space realization was obtained from a set of formance comparison. Accuracy will be quantifled by
impulseresponsesthatwereprocessedusingtheEigen- comparingLCOpanelresponse,(cid:176)ow-fleldpressuredis-
system Realization Algorithm (ERA).15,16 tribution on the elastic panel, LCO frequency, and
LCO phase for a variety of panel dynamic pressure
This research will extend the Volterra-POD ap-
values. Finally, the robustness of the Volterra-POD
proach to supersonic (cid:176)ow-flelds with dynamic bound-
method for predicting LCO response across a broad
ary behavior. The POD-Volterra method will be ap-
parameter space will also be addressed.
plied to a two-dimensional elastic panel in inviscid,
supersonic cross-(cid:176)ow. The Volterra-POD approach
Formulation
will be used to identify a low-dimensional, linear
POD/ROM for the (cid:176)uid. The POD/ROM will be This section describes the full-system aeroelastic
tightly coupled to a low-dimensional, nonlinear model model,introducesPOD,andoverviewsVolterrameth-
of the von K¶arm¶an plate equation.17 The aeroelas- ods. In addition, we fully develop the Volterra-POD
tic response will be obtained using an implicit time- approach and the synthesis of aeroelastic ROMs.
integration scheme.
Structural Dynamics Equations
The Volterra-POD technique involves procedures
Two-dimensional (cid:176)ow over a semi-inflnite, pinned
that require the selection of parameters such as im-
panel of length L is considered. Panel dynamics are
pulse size, several data windowing lengths, and im-
computed with von K¶arm¶an’s large-de(cid:176)ection plate
pulse sampling frequency. The choice of POD basis
equation, which is placed in nondimensional form us-
afiects performance as well. Some considerations for ing aerodynamic scales L and u 19 (0<x<1):
1
generating the POD basis include choice of base (cid:176)ow
(the POD/ROM determines the perturbation to this „@4w @2w @2w 1
¡N + =„ ¡p ; (1)
base (cid:176)ow), snapshot collection method and sampling ‚ @x4 x@x2 @t2 (cid:176)M2
(cid:181) 1 ¶
2 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
Fluid Dynamics
6„ h ¡2 1 @w 2 Thedynamicsofinviscid(cid:176)uid(cid:176)owsaregovernedby
N · 1¡”2 dx: (2) theEulerequations. Thetwo-dimensionalEulerequa-
x
‚ L @x
(cid:181) ¶ Z0 (cid:181) ¶ tions are given below in strong conservation form:20
¡ ¢
Thenonlinear,in-planeloadinEqn.(2),servestolimit
@U @E(U) @F(U)
panel de(cid:176)ections w(x;t) induced by (cid:176)uid-structure in- + + =0 ; (8a)
@t @x @y
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 U(x;t) = 2 mx 3; (8b)
although w and t are scaled by h and ‰ hL4 1=2, my
d d s 6 E 7
respectively. Two pinned boundary conditions are en- 6 T 7
forced at the panel’s endpoints: w =0 and¡ @2w =¢0. where ‰; m ; m and4E ar5e functions of space and
@x2 x y T
A modal solution for the de(cid:176)ection w(x;t) is as- time. Since we assume an ideal gas for our applica-
sumed: tions, this equation set can be closed using the ideal
ms
gas law.
w(x;t)= a (t)sin(i…x); (3)
i
The solution of the Euler equations can be approx-
i=1
X imated using either flnite-difierence, flnite-volume, or
where m is the number of structural modes retained,
s flnite-element techniques. To do this, the spatial do-
andthemodalamplitudesa varyintimeandarecol-
i mainisdiscretized,andthe(cid:176)owvariablesinU(x;t)at
located in the array a. The Galerkin method is used
eachdiscretelocationarecollocatedintoacolumnvec-
to obtain a low-order set of ordinary-difierential equa-
tor U(t). Time integration across the computational
tions describing the behavior of a .17 First, Eqn. (3)
i mesh is used to obtain (cid:176)ow solutions.
is substituted into Eqn. (1). The resulting expres-
Since the Euler equations are linear in the time
sionisthenintegrated,followingpre-multiplicationby
derivative, and quasi-linear in the spatial deriva-
sin(i…x), to yield (i=1;:::;m )
s tive,20,21 the spatial derivatives and the time deriva-
1 „(i…)4 6„ h ¡2 (i…)2 tives in Eqn. (8a) can be separated to form an evolu-
a˜i+ ai+ 1¡”2 fi ai =„Pi; tionarysystem. Toaccomplishthis,thespatialderiva-
2 2‚ ‚ L 2
(cid:181) ¶ tivesofthe(cid:176)uxterms @E and @F aregroupedtoform
¡ ¢ (4) @x @y
where fi· a2(r…)2 and a nonlinear operator R acting on the set of (cid:176)uid vari-
r r 2 ables. The(cid:176)uiddynamicsfromEqn. (8a)canthenbe
P 1 1 expressed as
P · ¡p sin(i…x)dx: (5)
i (cid:176)M2 dU(x;t)
Z0 (cid:181) 1 ¶ =R(U(x;t)) : (9)
dt
Theprojectedpressurecomponents,P ,areintegrated
i
fromtheaerodynamicsolutionwiththemidpointrule, When discretized this expression takes the form
using(cid:176)owfleldpressuresobtainedatgridpointsonthe
panelsurface.12 Theaerodynamicequations,theirdis- dU(t) =R(U(t)) : (10)
cretization, and their solution are discussed in later dt
sections. While equivalent to other formulations in Equation (10) is referred to as the full-system dynam-
the literature, Eq. (4) has two notable distinctions. ics.
First, the difierent form of scaling described above al- A flnite-volume scheme was the basis for the full-
ters equation coe–cients, and, second, an expression ordersolverusedinthisresearch,whichapproximated
relating p to the state of the panel is not assumed.12 the integral form of the Euler equations:
ThestructuraldynamicsequationEqn.(4)isplaced
in flrst-order form by introducing a mode speed array, d
UdV + (E{+F|)¢dS =0 : (11)
b, such that a_i =bi, dtZV Z@V
„(i…)4 6„ Li… 2 The grid points in the compbutatiobnal mesh described
b_i =¡ + 1¡”2 fi ai+2„Pi: earlier were used to form corners for cells. For each
‚ ‚ h
" (cid:181) ¶ # cell, the integral form of the Euler equations reduced
¡ ¢ (6)
to the following, assuming no grid deformation:
The mode speeds and amplitudes are collocated into
a structural solution array, Ys, leading to a general d dSi;j
U + (E {+F |)¢ =0 : (12)
form of the structural equation: i;j i;j i;j
dt dA
i;j
sides
X
Ys = [b;a]T (7a) The (cid:176)ux terms E abnd F bwere computed using
i;j i;j
Y_ s = Rs(Ys;P;„;‚;h=L): (7b) second-order Roe averaging,20 and the (cid:176)ow variables
3 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
Ui;j were evaluated as cell averages. Time integration Proper Orthogonal Decomposition
acrossthecomputationalmeshwasusedtoobtain(cid:176)ow POD is a technique to identify a small number of
solutions. This was accomplished with a flrst-order- basis functions that adequately describe the behav-
accurate, forward Euler approximation. ior of the full-system dynamics (Eqn. (10)) across
someparameterspaceofinterest. AsummaryofPOD
External boundaries were handled with ghost cells.
as it applies to a spatially-discretized (cid:176)ow fleld fol-
The (cid:176)uid values for the ghost cells at the far
lows. AdetaileddescriptionofPODisavailableinthe
fleld boundaries were determined using characteristic
literature.4,8 For simplicity, consider only one (cid:176)uid
boundary conditions.20 The bump-surface was mod-
variable, w(x;t), which when spatially descritized us-
eledusingatranspirationapproximation.22 Theflnite-
ing N nodes is denoted w(t). For this (cid:176)uid variable,
volume (cid:176)uid solver and the transpiration boundary
the full-system dynamics in Eqn. (10) is expressed as
condition were validated using a combination of the-
oryandexperimentaldata. Subsonicperformancewas dw
validated using wind-tunnel data.23 Supersonic per- =Rw(w) : (13)
dt
formance was validated using oblique shock theory.24
Time-accurateperformancewasvalidatedbycorrectly Spectralmethodsapproximatethesolutionw(X;t)as
predictingthevortexsheddingfrequencyfromacylin-
M
der in cross (cid:176)ow.
w(x;t)… a (t)` (x) : (14)
k k
k=1
X
Time Integration of the Coupled Full-Order
Equations When the domain is spatially discretized, `k(x) be-
comes a vector ` , and the following relation applies:
k
The systems of discretized (cid:176)uid dynamic equations,
U(t), and modal structural equations, Y , are com- M
s
bined into a single time-dependent system represen- w(t)… ak(t)`k : (15)
tative of the complete interaction between structure kX=1
and inviscid (cid:176)ow. Time integration proceeds in two
The set of vectors f` g are discrete basis functions
k
steps, assuming an O(¢t) lag in the synchronization
corresponding to the computational mesh deflned for
of (cid:176)uid and structure. First, the structural variables
thenumericalsolver. Thesetfa garethemodalcoef-
k
are updated from time level n to n+1 using a Crank-
flcients,andEqn. (15)canberepresentedusingmatrix
Nicolsonproceduretobedescribedbelow(butlimited
algebra. The(cid:176)uidmodescomprisecolumnsofamodal
heretoonlystructuralvariables). Duringthisstep,the
matrix ', and the coe–cients are collocated into a
pressuresknownatgridpointsonthepanelsurfaceare column vector w(t). POD produces a linear transfor-
consideredfrozen. Inthesecondstep,theaerodynamic mation ' between the full-order solution, w, and the
variables are explicitly updated using only structural reduced-order sbolution, w:
variables deflned at time level n.
w(t)=W +'w(t) : (16)
0b
Grid Generation and Time Step
The reduced-order variable w(t) represents deviations
b
The (cid:176)ow is simulated over a physical domain of of w(t) from a base solution W . The subtraction of
0
lengthDL,centeredaboutx=0,andheightDH. The W0 will result in zero-valuedbboundaries for the POD
domain is discretized using I nodes in the streamwise modes wherever constant boundary conditions occur
direction and J nodes normal to the panel. Indices i on the domain.
(1•i•I) and j (1•j •J) are used to denote grid ' is constructed by collecting observations of the
pointscomprisingcornersofcellsfortheflnite-volume solutionw(t)¡W atdifierenttimeintervalsthrough-
0
scheme. Gridpointsareclusteredinthedirectionnor- out the time integration of the full-system dynamics.
mal to the panel at the panel surface, with minimum Theseobservationsarecalledsnapshots18 andaregen-
spacingdenotedby¢ . Thespacingofgridpointsis erally collected to provide a good variety of (cid:176)ow fleld
wall
specifled to grow geometrically with j from the panel dynamics while minimizing linear dependence. The
boundary. Inthestreamwise direction, thenodespac- snapshot generation procedure is sometimes referred
ing is chosen to be uniform over the deforming panel to as POD training.8
segment (coincident with the structural grid), while A total of Q snapshots are collected from the full-
growing geometrically upstream of the leading edge system dynamics. These are vectors of length N. The
(positioned at i = I ) and downstream of the trail- set of snapshots describe a linear space that is used
LE
ing edge (positioned at i = I ). Calculations are to approximate both the domain and the range of the
TE
carriedoutwithabaselinegridgivenbythefollowing: nonlinear operator R . The linear space is deflned
w
I = 141, J = 116, D = 50, D = 25, I = 45, by the span of the snapshots.25 POD identifles a new
L H LE
I =97, and ¢ =0:0125. basisforthislinearspacethatisoptimallyconvergent4
TE wall
4 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
in the sense that no other set of basis functions will formulation is considered:
capture as much energy in as few dimensions as the
t
POD basis functions. To identify the POD basis, the X(t)=h + h (t¡¿)u(¿)d¿
0 1
snapshotsarecompiledintoanN£QmatrixS,known Z0
as the snapshot matrix. The mapping function ' is t t
then developed using + h (t¡¿ ;t¡¿ )u(¿ )u(¿ )d¿ d¿ : (20)
2 1 2 1 2 1 2
Z0 Z0
STSV = V⁄ ; (17a) The assumption of a weakly nonlinear system is con-
' = SV : (17b) sistent with the emergence of limit-cycle oscillation of
a 2-D aeroelastic system in transonic (cid:176)ow through a
Here V is the matrix of eigenvectors of STS, and ⁄ is supercritical Hopf bifurcation.29 For linear systems,
the corresponding diagonal matrix of eigenvalues. To only the flrst-order kernel is non-trivial, and there are
eliminate redundancy in the snapshots, the columns no limitations on input amplitude.
of V corresponding to very small eigenvalues in ⁄ are Theflrst-andsecond-orderkernelsarepresentedbe-
truncated. The matrix of eigenvalues ⁄ is also resized low in flnal form:5
to eliminate the rows and columns corresponding to
1
the removed eigenvalues. If Q¡M columns of V are h1(¿1) = 2X0(¿1)¡ X2(¿1) ; (21)
2
truncated, theresultingreducedordermapping'will
1
be an N £ M matrix. The reduced-order mapping h2(¿1;¿2) = (X1(¿1;¿2)¡X0(¿1)¡X0(¿2)()22:)
2
represents an isomorphism between N and M dimen-
sional space. It determines the coordinates of w(t) in In (21), X0(¿1) is the time response of the system to
terms of the M remaining basis functions, `k. a unit impulse applied at time 0 and X2(¿1) is the
The reduced-order mappings for each (cid:176)uid variable time response of the system to an impulse of twice
are developed separately, and individual S and V ar- unit magnitude at time 0. These response functions
raysarecollocatedasblocksintoalargersetofarrays, represent the memory of the system. If the system is
also denoted S and V, to form linear, then X2 = 2X0 and h1 = X0, which is why
the flrst-order kernel is referred to as the linear unit
U(t) … U0+'U(t) ; (18a) impulse response. The identiflcation of the second-
' = SV : (18b) order kernel is more demanding, since it is dependent
b on two parameters. Assuming ¿ >¿ in (22), X (¿ )
2 1 0 2
These versions of Eqn.’s (16) and (17b), respectively, istheresponseofthesystemtoanimpulseattime¿ .
2
apply to the entire set of (cid:176)uid variables. Timeisdiscretizedwithasetoftimestepsofequiv-
PODofthediscrete,panelpositionvectorw(x;t)! alent size. Time levels are indexed from 0 (time 0) to
w(t) and panel velocity vector s(t) = w_(t) is accom- n (time t), and the evaluation of X at time level n is
plished in a similar manner as described above for the denoted by X[n]. The convolution in discrete time is
(cid:176)uid system. Unlike the (cid:176)uid POD basis functions,
there is no base term subtracted from the snapshots N
X[n] = h + h [n¡k]u[k] (23)
when generating a structural POD basis. 0 1
k=0
X
Volterra Methods N N
Consider time-invariant, nonlinear, continuous- + h [n¡k ;n¡k ]u[k ]u[k ](2:4)
2 1 2 1 2
time,systems. Ofinterestistheresponseofthesystem kX1=0kX2=0
about an initial state X(0) = X due to an arbi-
0 The linearized and nonlinear Volterra kernels can
trary input u(t) (we take u as a real, scalar input)
betransformedintolinearizedandnonlinear(bilinear)
for t ‚ 0. As applied to these systems, Volterra the-
state-space systems that can be easily implemented
ory1,26{28 yields the response
into other disciplines such as controls and optimiza-
t tion.5,27 For linear dynamics, state-space realization
X(t) = h0+ h1(t¡¿)u(¿)d¿ (19) using the Eigensystem Realization Algorithm (ERA)
Z0 hasbeenusedtogeneratelinear,aeroelasticsystems.16
t t
+ h (t¡¿ ;t¡¿ )u(¿ )u(¿ )d¿ d¿ Nonlinear system realization is an active area of re-
2 1 2 1 2 1 2
Z0 Z0 search.
N t t
System Realization
+ :: h (t¡¿ ;::;t¡¿ )
n 1 n
n=3Z0 Z0 The ERA method15 identifles a discrete, linear,
X
u(¿ ):: u(¿ )d¿ :: d¿ : time-invariant state-space realization of the form,
1 n 1 n
The Volterra series can be accurately truncated be- X[n+1] = AX[n]+Bu[n]
yond the second-order term when a weakly nonlinear Y[n] = CX[n] ; (25)
5 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
using data from a complete ensemble of impulse re- where0 and0 asthenullmatricesoforderpandm
p m
sponses. Initial state responses can be used in lieu respectively, and I and I are the identity matrices
p m
of impulse responses, but we only consider impulse of order p and m.
response data in this overview for simplicity. The sys- Since the discrete time step ¢t = t ¡t is con-
k+1 k
tems realization procedure takes measurement data stant, the continuous form of the discrete state-space
Y[n] from the free response of the system and pro- realization (Eqn. (25)) is easily obtained. The contin-
duces a minimal state-space model A;B; and C such uous from, shown below, may require additional state
that functions Y are accurately reproduced. dimensionality when the discrete realization has real,
Thefreepulseresponseoflinear,time-invariant,dis- negative poles:
crete systems is given by a function known as the
X_ (t) = AX(t)+Bu(t)
Markov parameter,
Y(t) = CX(t) : (31)
Y [n]=CAn¡1B : (26)
m
Aeroelastic ROM Development
The superposition principle states that a system re-
The full-order vector of (cid:176)uid variables U(t) repre-
sponse to any arbitrary input can be obtained from
sents the spatially discretized (cid:176)ow fleld obtained from
a linear combination of impulse responses from that
the full-system (cid:176)ow solver. POD provides a trans-
system. The generalized Hankel matrix of impulse
formation ' that maps U(t) to a low order vector of
responses is related to the Markov parameter by the
modalcoe–cientsU(t)(fromEqn. 18a). Thereduced-
superposition principle. The Hankel matrix is formed
order (cid:176)uid variable U(t) will be denoted Y , which is
by windowing the impulse response data. A total f
the vector of outpubts Y(t) (Eqn. (31)) .
of K data points are provided at discrete time steps
A state-space modbel for Y can be obtained from
n = 1;:::K, and the r £ s matrix H is formed as f
rs
impulse responses using the ERA method. Impulses
follows,
for the (cid:176)uid system use the plate position and ve-
Ym[n] ::: Ym[n+s¡1] locity coe–cients (Ys) as the forcing terms. Each
Ym[1+n] ::: Ym[1+n+s¡1] structural term is impulsed, and the (cid:176)uid system re-
Hrns¡1 =2 .. .. .. 3 (27) sponse is generated using the full-order model. The
. . .
66Ym[r¡1+n] ::: Ym[r¡1+n+s¡1]77 time history of the impulse response is projected onto
4 5 eachofthePODbasisfunctionstoobtaintheimpulse
where s is the total size of the data window, and r is response of the reduced-order (cid:176)uid vector Y . POD
f
the number of time steps used to shift the data win- basisfunctionsareobtainedusingthemethodofsnap-
dow. The choice of r and s is arbitrary as long as shots as described previously. The process is repeated
r+s+n•K+2. foreachstructuralmode,andthecollectionofimpulse
The ERA method eliminates redundant data by us- responsesisusedtogeneratealinearstate-spacemodel
ing Singular Value Decomposition (SVD) on H0 , for the reduced order (cid:176)uid system,
rs
H0 =PDQT : (28) X_ = A X +B Y ; (32a)
rs f f f f s
Y = C X ; (32b)
Unwanted state dimensionality is eliminated by trun- f f f
cating the elements of P;D; and Q associated with
where X is the state vector from the ROM realiza-
very small singular values of H0 . The number of f
rs tion, and Y represents the modal coe–cients for the
s
states is reduced to a minimal number q. The number
structural deformation.
of observations p and the number of forcing terms m
The structural model from Eqn. (7b) is coupled to
are known from the problem formulation. The dimen-
thereduced-order(cid:176)uidmodel(Eqns. (32a)and(32b))
sionoftheMarkovparameterY [n]isp£m. Algebra
m through the projected pressure term P. The reduced-
isusedtorecastEqn. (26)intermsofthetimeshifted
order variables, Y , are obtained from the dynamic
Hankel matrix H1 , and the elements P;D; and Q. s
rs statesX bythelinearmappingC . Fluidpressureon
f f
The state-space realization (cid:176)ows from this manipula-
thepanelisextractedfromY usingtheportionofthe
s
tion, and is as follows:
reduced order mapping (Eqn. (18a)) that pertains to
A = D¡12PTHr1sQD¡21 ; (29a) the moving boundary. Equation (5) is used to project
the pressure values into P. The mapping of reduced
B = D12QTEm ; (29b) order (cid:176)uid states to projected pressure is denoted f ,
P
C = EpTPD21 : (29c) P =f (X ) : (33)
P f
ET and ET are deflned below:
p m Equation(33)isusedtocouplethestructureand(cid:176)uid
ET =[I ;0 ;:::;0 ] ; (30a) dynamic state variables,
p p p p
ET =[I ;0 ;:::;0 ] ; (30b) Y_ =R (Y ;f (X );„;‚;h=L) : (34)
m m m m s s s P f
6 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
Equations (32a) and (34) comprise the low order, this fashion since any forcing function can be assem-
aeroelastic ROM with linear (cid:176)uid dynamics, and non- bled from a series of impulses.
linear plate dynamics. The supersonic (cid:176)ow fleld was essentially linear10
The (cid:176)uid and structural terms are grouped into ar- once LCO was fully developed. Shock waves formed
rays Y and R as follows: at the ends of the panel, and although they varied in
strengthdynamicallywiththeLCO,theywerestation-
X
Y = f ; (35) ary,andthe(cid:176)owfleldbetweentheshocks(directlyover
Y
s the panel) was linear for this case. However, impuls-
• ‚
A X +B Y ing the uniform, supersonic (cid:176)ow fleld with structural
R = f f f s : (36)
R (Y ;f (X );„;‚;h=L) modes(andmodalvelocities)producedverynonlinear
s s P f
• ‚
transientbehavior. Figure1showsdensitycontoursof
b
The reduced-order, aeroelastic system is denoted as the (cid:176)ow 1:10802 nondimensional time units from the
simply, impulse of the fundamental structural velocity mode.
The sudden appearance of a velocity proflle on the
Y_ = R(Y) : (37)
Time Integration of the Aebroelastic System 2
Density
1.00011
The aeroelastic ROM (Eqn. (37)) is integrated in 1.5 11..0000000069
time with the two-time level, second-order accurate, 1.00004
1.00002
Y 1 0.999997
Crank-Nicolson method: 0.999975
0.999953
0.99993
0.5 0.999908
Yn+1¡Yn 1
= R^n+1+R^n ; (38)
¢t 2 0-2 -1 0 X 1 2 3
‡ ·
¢t ¢t
R·Yn+1¡ R^n+1¡Yn¡ R^n =0: (39) 2
2 2 L1e0vel1D.0e0n0s1it1y
At each time level, Y ·Yn+1 is computed from Eqn. 1.5 987 111...000000000000964
6 1.00002
(b3ia9n) usiIng¡a¢cthJo^ord teYckh+n1iq¡ueYwkith=a¡tRim(Ye-kfr)o;zen Ja(4c0o)- Y 0.501-2 54321 00000.....99999-9999919999999999975307538 67676767675065767544X64 56256471663154646544568453857526175566243626545 5 553
2
(cid:181) ¶
¡ ¢ Fig. 1 Density contours 1:10802 time units after
where k is a subiteration index and J^ is the Jacobian impulse
o
of the reduced order aeroelastic system, evaluated for
the base (cid:176)ow condition and Y =0. A suitable num- boundary(anditssuddenremovalonetimesteplater)
s
producedashockwaveofvaryingstrengthrunningthe
berofsubiterationsarecomputedateachtimestepto
obtain a good approximation to Yn+1; typically, 1-2 lengthofthepanel. Earlyduringthetransientperiod,
this shock welled upward away from the panel, and
subiterationsaregenerallysu–cienttodriveRtonear
convected downstream. The convection of the shock,
machine zero. Since peak panel de(cid:176)ection is no more
combined with the patterns of varied intensity pro-
than 2% of panel length for the cases considered, the
duce the odd (but physical) spatial oscillation above
chord method is rapidly convergent. Prior to subiter-
thepanelshowninFig. 1. Afterabout2:5 timeunits,
ation, Y is predicted from the explicit formula
this pattern had both convected well beyond the end
Y =Yn+¢tR^n: (41) ofthepanel,anddifiusedintoamorebenign(cid:176)owpat-
tern. After 25 time units uniform (cid:176)ow was restored.
Results The(cid:176)owdynamicswereessentiallylinearaftertheini-
tial 2:5 time unit transient.
The results that follow consider supersonic free
The results that follow will detail the usefulness of
stream(cid:176)owconditionsatMach1:2,withsealevelcon-
such impulse responses for generating a reduced-order
ditions. The Galerkin panel model contains 4 modes,
(cid:176)uidmodel. Thelinearportionoftheimpulseresponse
for a total of 8 DOFs.
timeintegrationscontainedsomeoftheLCO(cid:176)owfleld
Impulse Response of a Supersonic Fluid characteristics, otherwise none of the ROM options
Impulsing the forcing term (or terms) of a truly lin- would have reproduced LCO. However, the impulse
ear system produces a response that is the building responses themselves would not produce POD basis
block necessary to recreate the system output from functions capable ofcorrectlymodelling theLCO(cid:176)ow
any arbitrary forcing function. Linear superposition fleld. This suggests that the supersonic (cid:176)ow fleld was
allows the response of the system to be constructed in not truly linear.
7 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
Identiflcation of Fluid Modes data was windowed using s=192 and r =100. Every
flfth data point was used in the realization algorithm
Fluid modes were obtained using POD as outlined
providing data at the rate of dt=0:07716. These val-
previously. Aeroelastic (cid:176)uid modes were obtained
ues of dt;s and r were chosen by trial and error to
from a set of 100 snapshots. Snapshots of the full-
produce realizations whose impulse responses closely
order, aeroelasticsystemweretakenatequallyspaced
matched the data. The value of m was 8 to match
intervals, from start-up through 25 time units, using
the number of forcing terms, and p=8 was chosen to
timeintegrationofthefullsystem,with‚=25. Forall
matchthenumberofROMcoe–cients. Thecollection
time integration cases that follow, the aeroelastic sys-
ofimpulseresponsesformedan8£8Markovparameter
tem was initialized with the base (cid:176)ow condition, and
Y [n]functionfromEqn. (26). Thenumberofstates,
a small perturbation (height of 0:0001) in the funda- m
q = 8, was chosen to match the number of ROM co-
mental panel position mode (denoted a in Eqn. (7a).
1
e–cients, so SVD on the Hankel matrix formed from
Snapshots were taken of primitive (cid:176)uid variables ‰, u,
Y[n = 1] was used to truncate all but the largest 8
v, and E , not the conserved variables given in Eqn.
T
singular values, yielding the matrices P; D; and Q.
(8b). Primitive variables enable Galerkin projection
Equations (29a, 29b, 29c) were then used to generate
for the (cid:176)uid as a possible means for obtaining non-
linear terms13 in future analysis. At Mach 1.2, LCO a linear state-space model for the reduced-order (cid:176)uid
system,
required about 300 time units to become fully devel-
oped,andthesmall,25timeunittrainingwindowwas x[n+1] = Ax[n]+B fl u[n] ; (42a)
shown to be adequate in previous work.10,12
y[n] = Cx[n] ; (42b)
The results that follow refer to two cases. For the
flrstcase,thebase(cid:176)owtermU fromEqn. (18a)con- wherefl was ascaling parameter that wasusedtocal-
0
sisted of uniform, free-stream conditions everywhere ibrate the forcing amplitude.
throughout the domain (referred to as slug (cid:176)ow). The The value of the scaling parameter fl = 650 was
second case considered steady state (cid:176)ow over the ini- set by tuning the ROM results to the snapshot data.
tial panel perturbation as the base (cid:176)ow term U . For Theoretically, fl should have been the inverse of the
0
both cases, the flrst two modes for each (cid:176)uid variable impulse size (fl = 1=0:1 = 10). The need for an or-
contained over 98% of the energy content, and sys- der magnitude increase in fl reveals an ine–ciency in
tem realization was performed using a total of M =8 projecting the impulsed (cid:176)ow-fleld onto the aeroelastic
(cid:176)uid modes (2 modes per (cid:176)uid variable). We also at- modes. Evidently, the impulsed (cid:176)ow-fleld contained
temptedtousetheimpulseresponsedataassnapshots, structures not adequately represented in the aeroe-
in lieu of aeroelastic time integration. The same sys- lastic modes, and a signiflcant amount of (cid:176)ow energy
tems realization procedure was repeated to produce was not captured in the projections used to compute
a Volterra-POD ROM for this third case, but this the modal-impulse behavior. However, enough linear,
ROMdidnotcorrectlyproduceLCO.Wedocumented aeroelastic information was resident in the impulsed
ourobservations,andrecommendationsregardingthis (cid:176)ow-fleld for the Volterra-POD realization to produce
third approach in a separate section. correct results (with fl properly adjusted).
System Realization 0.0015
We considered 8 POD modes, the dimensionality of 0.001 Fullsystem
1 POD/ROM
Y , which produced eight impulse responses for each de0.0005
f o
M
forcingfunction. With8forcingtermsinYs,thetotal sity 0
number of impulse responses numbered 64. Realiza- en-0.0005
D
tion via the ERA process for each of the aeroelastic -0.001
cases is detailed below. -0.00150 5 10 15
Uniform Base Flow 0.0015
A state-space realization of the form in Eqn. (25) 20.001
was obtained using ERA for the slug-(cid:176)ow base case. ode0.0005
M
The impulse amplitude was arbitrarily chosen to be sity 0
0:1. Thefull-systemresponsetothisimpulsewassam- en-0.0005
D
pled over 30 non-dimensional time units at a rate of -0.001
dt = 0:015432 for a total of K = 1944 discrete data -0.00150 5 10 15
points. The (cid:176)uid system impulse response was gen-
Fig.2 Responseofdensitymodestovelocityterm
erated using the full-order model. The time history
impulse, uniform base (cid:176)ow case
of the full-order impulse response was projected onto
eachofthePODbasisfunctionstoobtaintheimpulse The impulse response of Eqns. (42a and 42b) was
response of the reduced-order (cid:176)uid variable Y . The obtained with fl = 1. The impulse responses of the
f
8 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
reduced-ordersystemwereingoodagreementwiththe ROM Time History
impulse responses from the full-order system used by
Both the slug-(cid:176)ow base case and the steady-state
ERA. The response of the two density modes to an
basecaseVolterra-PODROMs,describedabove, were
impulse in the flrst panel velocity term b within Y
1 s
time-integrated using the aeroelastic training condi-
(deflned in Eqn. (7a)) are shown in Fig. 2. The con-
tions of Mach 1:2 and ‚ = 25. The results are shown
tinuous form of Eqns. (42a and 42b) were obtained
in Fig. 4. Both cases correctly predicted LCO, but
via a function call in MATLAB, which provided the
matrices A , B and C for time integration of the
f f f
aeroelastic model given in Eqn. (37).
Steady-state Base Flow Chord) 00..34 F8uMlloSdyestReOmM,Slug-FlowBaseCase
4 0.2
3/
( 0.1
l(cid:176)aosTwtihcceomnsoaddmiteiesonwpredoreecseccdoruimbreepduwteaeadsrlruieesprin.egaFtaeodrs,ttehbaiudstyr-teshataelitzeaaetbriaoosnee-, elDeflection ---000...0321 TrRaOinMing
n
every20th data point wasusedinthe realization algo- Pa -0.4
0 100 T2im00e 300 400
rithm providing data at the rate of dt = 0:3086. The
asanmdethimepfuullls-esyasmtepmlitiumdpeuolsfe0:r1eswpaosnsuesewdasfosratmhpislecdasaet, Chord) 00..34 8FuMlloSdyestReOmM,Steady-StateBaseCase
the same rate. However, the impulses were added to 3/4 0.2
( 0.1
ttdhhaeetavsatweluaasedywo-fisntmadtoewwaepsdan8uestliondgme(cid:176)saet=ccthi4ot7nheafonnrdutmrhb=iser2ca0os.fef.AorgcTaiinhnge, elDeflection ---000...0321 TrRaOinMing
n
terms, and p = 8 was chosen to match the number of Pa -0.4
0 100 T2im00e 300 400
ROM coe–cients. The scaling parameter from Eqn.
(42a) was fl =800, and was determined by tuning the Fig.4 Panelde(cid:176)ection(wd=h)timehistory, ‚=25,
ROM to the snapshot data. Mach 1:2
the steady-state base case was more accurate in am-
0.004
plitude, frequency and phase than the slug-(cid:176)ow base
FullSystem
e1 0.002 POD/ROM case. The 3=4 chord panel-amplitude was muted by
d
o
M 0 9% for the steady-state base case, and magnifled by
y
sit 25% for the slug-(cid:176)ow base case. The phase and fre-
en -0.002
D quency error were negligible for the steady-state base
-0.004 case, but the larger amplitudes of the slug-(cid:176)ow base
5 10 15
caseintroducedasmallincreaseinLCOfrequency(re-
0.004 sulting in an accumulating phase error).
2 0.002 We suspect the improvement in performance asso-
e
od ciated with the steady-state base (cid:176)ow was most likely
M 0
sity due to the choice of data windowing parameters used
en -0.002 in the ERA realization. We noticed no substantive
D
-0.004 difierences in either the POD modes, or the impulse
5 Time 10 15 response of the full system (cid:176)ow-fleld between either
case. Data windowing parameters were selected to
Fig.3 Responseofdensitymodestovelocityterm
provide a realization whose impulse response closely
impulse, steady-state base (cid:176)ow case
matched the original impulse data. As noted previ-
ously,reducingthesizeofdtintroducedhigh-frequency
The larger value of dt eliminated much of the high- dynamics into the realization. The slug-(cid:176)ow base case
frequency transient, and focused ERA on the low- was formed using a very small value of dt. Conse-
frequency portion of the impulse response. This is quentlytheimpulseresponseofthemodel(seeFig. 2)
re(cid:176)ected in the impulse response accuracy shown in matched the initial transient in the data better than
Fig. 3,whichconsiderstheresponseofthetwodensity the steady-state base (cid:176)ow case (i.e. Fig. 3), which
modestoanimpulseintheflrstpanelvelocitytermb used a much larger dt. However, the high-frequency
1
within Y (deflned in Eqn. (7a)). The flrst 3 seconds data in the impulse response was not germane to the
s
of the impulse response curve was not matched very large-time behavior of the aeroelastic system, and its
closely by the ROM; however, these initial transients inclusionresultedinalessaccurateROMunderaeroe-
were not important to the LCO (cid:176)ow fleld. lastic conditions.
9 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
ROM Robustness structural model. Subspace projection relied on the
Both aeroelastic ROMs were time-integrated using Galerkin panel model for time integration. The panel
a variety of panel dynamic pressure values (‚). The position and velocity from the Galerkin panel model
intent was to explore the predictive accuracy of the were projected onto the POD basis functions at ev-
Volterra-POD ROMs across a nonlinear parameter ery step in the time integration. Subspace projection
space. Both ROMs were trained at ‚ = 25 (as de- demonstratedtheadequacyofthePODmodesatcap-
scribed earlier), and robustness was deflned as the turing the dynamic panel behavior, while maintaining
ROM’s ability to predict panel amplitude (at the 3=4 the nonlinearity of the panel dynamics.
chord position) in fully developed LCO across the pa-
FullSystem
rameter space, including non-LCO cases. hord) 0.3 888MMMooodddeeeRRROOOMMM,,,442MDDOOodFFeSS,tt8rruuDccOttuuFrreeG((aPPlOOerDDki))nStructure
C 0.2
4
(3/ 0.1
1 MMaacchh11..22((FGuolrldSnyiestreamn)dVisbal(2000)) ction 0
0.9 MMaacchh11..22((88--MMooddeeRROOMM,,SStluegadFyloSwtaBteasBeasCeasCea)se) nelDefle --00..21 TrRaOinMing
ord)0.8 Pa -0.3 100 T2im00e 300 400
h
C
3/40.7 hord)
w/h(d0.6 (3/4C 0.3
Deflection,000...345 anelDeflection 00..12
anel0.2 P 320 340 Time 360 380
P
0.1 Fig. 6 Panel response (wd=h) with aeroelastic
structural modes
0
0 10 20 30 40 50
Thereduced-orderstructuralmodelwastightlycou-
NondimensionalDynamicPressure,l
pled with the steady-state base case (cid:176)uid ROM and
Fig. 5 Panel response verse dynamic pressure time-integrated using the parameter value ‚ = 25.
Figure 6 compares the results with the full system
The ROM results are compared with full system re- response, and the POD/ROM results from Fig. 4
sults, and results from the literature30 in Fig. 5. The (for the steady-state base case). For clarity, the top,
steady-statebasecaseROMwasbettersuitedtolarge right-hand corner of the entire time response is ex-
LCOamplitudeswherethelargerpanelamplitudesex- panded in the lower portion of Fig. 6. Two POD
cited a panel nonlinearity that corrupted the results modes per structural variable (4 DOFs total) yielded
of the slug-(cid:176)ow base case. Conversely, the use of slug essentially identical results to the 4 mode (8 DOF)
(cid:176)ow as the base-(cid:176)ow term permited a more accurate Galerkin result. Further order-reduction greatly de-
prediction of LCO onset at the lower values of ‚. A creased the panel response. For this problem, the
non-LCO solution for the slug-(cid:176)ow base case occurred full-system structural model was very low order, and
when Yf = [0] and Ys = [0], but the use of steady- theadditionalorderreductionfromPODwasimmate-
state(cid:176)owoveraninitial,non-zerovalueofYsrequired rial. However,futureapplicationofthistechniquewill
a non-zero value of Yf to produce Ys = [0]. The involve very high-order, nonlinear structural models
steady-state base case had di–culty producing this requiring order-reduction along with the (cid:176)uid model.
result. In addition, the small panel amplitudes near These results demonstrate that a single training event
LCO-onset did not excite the high-frequency errors in can produce adequate POD modes for both the (cid:176)uid
theslug-(cid:176)owbasecasethatwereevidentatlargerval- and structure.
ues of ‚.
ROMs Using Impulse Response Modes
Aeroelastic Structural Modes As an excursion, the aeroelastic (cid:176)uid modes were
Aeroelastic structural modes were generated using replaced with POD modes derived from the impulse
100snapshotsofthestructuralresponseobtaineddur- responsesofthefullsystem. Fortysnapshotswerecol-
ing the training of the (cid:176)uid ROM. The structural lected from each of the eight impulse responses gener-
snapshots corresponded exactly in time with the set ated by impulsing the elements of Y . The snapshots
s
of 100 snapshots used to construct the (cid:176)uid ROMs. were taken at even intervals over a time integration
Snapshots were taken of the panel position and ve- lasting 30 time units. Since there were 8 impulse re-
locity vectors (w(t) and w_(t) respectively), and sub- sponses,atotalof320snapshotsweregenerated. Over
space projection8 was used to form a reduced-order 98% of the (cid:176)ow energy was contained in the flrst 4
10 of 12
American Institute of Aeronautics and Astronautics Paper 2003-1922
