Aeroelastic 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