AEROM: NASA’s Unsteady Aerodynamic and Aeroelastic Reduced-Order Modeling Software

The origins, development, implementation, and application of AEROM, NASA’s patented reduced-order modeling (ROM) software, are presented. Full computational fluid dynamic (CFD) aeroelastic solutions and ROM aeroelastic solutions, computed at several Mach numbers using the NASA FUN3D CFD code, are presented in the form of root locus plots in order to better reveal the aeroelastic root migrations with increasing dynamic pressure. The method and software have been applied successfully to several configurations including the Lockheed-Martin N+2 supersonic configuration and the Royal Institute of Technology (KTH, Sweden) generic wind-tunnel model, among others. The software has been released to various organizations with applications that include CFD-based aeroelastic analyses and the rapid modeling of high-fidelity dynamic stability derivatives. Recent results obtained from the application of the method to the AGARD 445.6 wing will be presented that reveal several interesting insights.


I. Introduction
In the early days, aeroelasticians typically used linear methods to compute unsteady aerodynamic responses and subsequent aeroelastic analyses. 1 These aeroelastic analyses were usually presented in the form of velocity-damping-frequency (V-g-f) plots or aeroelastic root locus plots as a function of either dynamic pressure or velocity. These plots were generated rapidly and provided significant amount of insight regarding the aeroelastic mechanisms involved.
With the subsequent development of Computational Fluid Dynamics (CFD) methods, the ability to analyze complex, nonlinear flows, and their effect on the aeroelastic response, became a reality. While CFD tools are quite powerful and provide significant insight regarding flow physics, the significant increase in computational cost (time and CPU dollars) has had an effect on how aeroelastic analyses are performed. One side-effect of the increase in computational cost is the desire to keep the number of time steps computed, and the total number of solutions generated, at a minimum. Results are, therefore, computed for a small number of dynamic pressures (per Mach number) with only a few cycles computed per dynamic pressure. A second side-effect is that the resultant time histories of each generalized coordinate cannot be directly used to identify the governing aeroelastic mechanisms at work, as was the case for the classical linear methods (V-g-f and root locus plots, for example). However, the recent development of reduced-order modeling (ROM) methods, 2-4 provides a tool for the rapid generation of traditional aeroelastic tools such as the root-locus plots.
The origin of this method started with the author's PhD dissertation 5 and related publications. 6,7 An important conceptual development first presented in these references consists of the realization that unsteady aerodynamic impulse responses do, in fact, exist and can be computed using CFD methods. This concept is an important point that is claimed to be not realizable in some of the classic aeroelastic references. The reason for this discrepancy is actually quite simple as it relates to the difference between the impulse function for a continuous-time system versus that for a discrete-time system.
For a continuous-time system, it is well known that the impulse input function is the Dirac delta function. This function serves the continuous domain well, in particular in the solution of ordinary and partial differential equations. However, its application to a discrete-time system such as a CFD-based solution, is not clear, thus the belief that an impulse input cannot be applied to a CFD code. Therefore, if an impulse input cannot be applied to a CFD code, then an unsteady aerodynamic response cannot be identified or realized.
An important contribution by the author 5 is the realization that in order to properly identify the unsteady aerodynamic impulse response using a CFD code, a discrete-time impulse input, also known as the unit sample input in discrete-time theory, is the proper function to use and not the Dirac delta function. The theory of Digital Signal Processing (DSP) demonstrates that a unit sample input is much simpler to apply and less complex to interpret than the Dirac delta function. These results proved the existence and realizability of a unit unsteady aerodynamic impulse (sample) response via a CFD code.
In the world of structural dynamics and modal identification, the concept of a structural dynamic impulse response is clear and well understood. As a result, various modal identification techniques consist of the identification of these responses and a subsequent realization of a system that captures the structural dynamic system of interest. Having familiarity with one of these methods by the name of Eigensystem Realization Algorithm (ERA) 8 /System Observer Controller Identification Toolbox (SOCIT), 9 the author applied the modal identification technique, previously limited to structural dynamic systems, to that of identifying an unsteady aerodynamic system via the identification of the unsteady aerodynamic impulse responses. Once the concept of a discrete-time unsteady aerodynamic impulse response was mathematically validated, the application of ERA/SOCIT became quite logical. 10 These results 10 represent the first time that the ERA/SOCIT algorithms were used for the identification of unsteady aerodynamic systems. It is valuable to point out that this method is now being applied at several organizations around the world. [11][12][13][14][15][16] In the area of fluid modal decompositions using, primarily, the Proper Orthogonal Decomposition (POD), the application of the ERA algorithm has become standard, with an initial appearence in the literature by Ma, Ahuja, and Rowley. 17 Following these fundamental advances, the development of linearized, unsteady aerodynamic state-space models for prediction of flutter and aeroelastic response using the CFL3Dv6 18 code was introduced. 19 Unsteady aerodynamic state-space models were generated using the ERA and coupled with a structural model within a MATLAB/SIMULINK TM environment for rapid calculation of aeroelastic responses, including the prediction of flutter. Aeroelastic responses computed directly using the aeroelastic simulation ROM showed excellent comparison with the aeroelastic responses computed using the CFL3Dv6 code.
The aerodynamic impulse responses (unit pulses) that were used to generate the unsteady aerodynamic state-space model 19 were computed by the excitation of one mode at a time. However, for more realistic cases where the number of modes can be an order of magnitude or more larger, the one-mode-at-a-time method becomes prohibitively expensive. Kim et al 20 proposed methods based on the simultaneous application of structural modes as CFD input, greatly reducing the cost of identifying the aerodynamic impulse responses from the CFD code. Silva 2 developed a method that enables the simultaneous excitation of the structural modes using orthogonal functions. Both of these methods require only a single CFD solution and the methods are independent of the number of structural modes.
A method for generating static and matched-point aeroelastic solutions, using a ROM, have also been developed. 21 These methods 2, 21 have been implemented in the FUN3D 22-25 CFD code. Methods for generating root locus plots of the aeroelastic system (combined structural and unsteady aerodynamic state-space models) have also been developed. 3 Applications of these methods include fixed-wing and launch vehicle configurations. 4 This paper will discuss the application of these methods to three configurations: a low-boom configuration, a full-span wind-tunnel model, and the AGARD 445.6 wing.
The AEROM software was granted a patent (November, 2011), Patent No. 8,060,350. The software has been distributed to the Air Force Research Laboratory, the Boeing Corporation, and the CFD Research Corporation.

II.A. FUN3D Code
The unstructured mesh solver used for this study is FUN3D. Within the code, the unsteady Navier-Stokes equations are discretized over the median dual volume surrounding each mesh point, balancing the time rate of change of the averaged conserved variables in each dual volume with the flux of mass, momentum and energy through the instantaneous surface of the control volume.
Because the CFD and computational structural mechanics (CSM) meshes usually do not match at the boundary interface where the grids are defined, CFD/CSM coupling requires a surface spline interpolation between the two domains. The interpolation of CSM mode shapes to CFD surface grid points is done as a preprocessing step. 26 Modal deflections at the CFD surface grids are first generated. Mode shape displacements located at CFD surface grid points are used in the integration of the generalized modal forces and in the computation of the deflection of the deformed surface. The final surface deformation at each time step is a linear superposition of all the modal deflections.

II.B. System Identification Method
In structural dynamics, the realization of discrete-time state-space models that describe the modal dynamics of a structure has been enabled by the development of algorithms such as the ERA 8 and the Observer Kalman Identification (OKID) 27 Algorithm. These algorithms perform state-space realizations by using the Markov parameters (discrete-time impulse responses) of the systems of interest. These algorithms have been combined into one package known as SOCIT developed at NASA Langley Research Center.
There are several algorithms within the SOCIT that are used for the development of unsteady aerodynamic discrete-time state-space models. The PULSE algorithm is used to extract individual input/output impulse responses from simultaneous input/output responses. For a four-input/four-output system, simultaneous excitation of all four inputs yields four output responses. The PULSE algorithm is used to extract the individual sixteen (all combinations of four inputs and four outputs) impulse responses that associate the response in each of the outputs due to each of the inputs. Once the individual sixteen impulse responses are available, they are then processed via the ERA in order to transform the sixteen individual impulse responses into a four-input/four-output, discrete-time, state-space model. A brief summary of the basis of this algorithm follows.
A finite dimensional, discrete-time, linear, time-invariant dynamical system has the state-variable equations where x is an n-dimensional state vector, u an m-dimensional control input, and y a p-dimensional output or measurement vector with k being the discrete time index. The transition matrix, A, characterizes the dynamics of the system. The goal of system realization is to generate constant matrices (A, B, C, D) such that the output responses of a given system due to a particular set of inputs is reproduced by the discrete-time state-space system described above.
For the system of Eqs. (1) and (2), the time-domain values of the discrete-time impulse responses of the system are also known as the Markov parameters and are defined as with A an (n x n) matrix, B an (n x m) matrix, C a (p x n) matrix, and D an (p x m) matrix. The ERA algorithm begins by defining the generalized Hankel matrix consisting of the discrete-time impulse responses for all input/output combinations. The algorithm then uses the singular value decomposition to compute the (A, B, C, D) matrices. In this fashion, the ERA is applied to unsteady aerodynamic impulse responses to construct unsteady aerodynamic state-space models.

II.C. Simultaneous Excitation Input Functions
Clearly, the nonlinear unsteady aerodynamic responses of a flexible vehicle comprise a multi-input/multioutput (MIMO) system with respect to the modal inputs and generalized aerodynamic outputs. In the situation where the goal is the simultaneous excitation of such a MIMO system, system identification techniques [28][29][30] dictate that the nature of the input functions used to excite the system must be properly defined if accurate input/output models of the system are to be generated. The most important point to keep in mind when defining these input functions is that these functions need to be different, mathematically, from each other. If the excitation inputs are identical, for example, and are applied simultaneously, it is quite difficult to separate the effects of one input from the others. This, in turn, makes it practically impossible for a system identification algorithm to extract the individual impulse responses for each input/output pair. As has already been well established, the individual impulse responses for each input/output pair are necessary ingredients towards the development of state-space models.
With respect to unsteady aerodynamic MIMO systems, these individual impulse responses correspond to time-domain generalized aerodynamic forces (GAFs), critical to understanding unsteady aerodynamic behavior. The Fourier-transformed version of these GAFs are the frequency-domain GAFs, that provide an important link to more traditional frequency-domain-based unsteady aerodynamic analyses.
Referring back to the input functions used to excite the MIMO system, the question is how different should these input functions be from each other and how can we quantify a level of difference between them? Since orthogonality (linear independence) is the most precise mathematical method for guaranteeing the difference between signals, recent developments focused on the application of families of orthogonal functions as candidate input functions. Using orthogonal functions directly provides a mathematical guarantee that the input functions are as different from each other as mathematically possible. These orthogonal input functions can be considered optimal input functions for the identification of a MIMO system.
In a previous paper, 2 four families of functions were investigated to efficiently identify a CFD-based unsteady aerodynamic state-space model. For the present paper, the Walsh family of orthogonal functions 31 are used, shown in Figure 1 for four modes. These functions are orthogonal and therefore provide a benefit in the system identification process as discussed above. Also, this family of functions consists of a combination of step functions, which have been shown to be well-suited for the identification of CFD-based unsteady aerodynamic ROMs.

III. ROM Development Processes
The ROM development process consists of two parts: the creation of the unsteady aerodynamic ROM and the creation of the structural dynamic ROM. The combination of the unsteady aerodynamic ROM with the structural dynamic ROM yields what is referred to as the aeroelastic simulation ROM.
The original unsteady aerodynamic ROM development process consisted of the excitation of one structural mode at a time per CFD solution. That approach is not practical for realistic configurations with a large number of modes. As mentioned above, an improved method has been developed and is described below.

III.A. Improved ROM Development Process
An outline of the improved simultaneous modal excitation ROM development process is as follows: 1. Generate the number of functions (from a selected family of orthogonal functions) that corresponds to the number of structural modes; 2. Apply the generated input functions simultaneously via one CFD execution resulting in GAF responses due to these inputs; these responses are computed directly from the restart of a steady rigid CFD solution (not about a particular dynamic pressure); 3. Using the simultaneous input/output responses, identify the individual impulse responses using the PULSE algorithm (within SOCIT);

Transform the individual impulse responses generated in
Step 3 into an unsteady aerodynamic statespace system using the ERA (within SOCIT); 5. Evaluate/validate the state-space models generated in Step 4 via comparison with CFD results (i.e., ROM results vs. full CFD solution results); A schematic of steps 1-4 of the improved process outlined above is presented as Figure 2.
Using modal information (generalized masses, modal frequencies, and modal dampings), a state-space model of the structure is generated. This state-space model of the structure is referred to as the structural dynamic ROM (Figure 3). Once an unsteady aerodynamic ROM and a structural dynamic ROM have been generated, they are combined to form an aeroelastic simulation ROM (see Figure 4). Then root locus plots are extracted from the aeroelastic simulation ROM.  An important difference between the original ROM process and the improved ROM process is stated in step (2) of the outline above. For the original ROM process, if a static aeroelastic condition existed, then a ROM was generated about that selected static aeroelastic condition. So a static aeroelastic condition of interest would be defined (typically based on a dynamic pressure) and that static aeroelastic condition was generated using the CFD code as a restart from a converged steady, rigid solution. Once a converged static aeroelastic solution was obtained, the development of the unsteady aerodynamic ROM process was applied about that static aeroelastic condition. This approach implies that the resultant unsteady aerodynamic ROM is limited to the neighborhood of that static aeroelastic condition. Moving too far away from that condition could result in loss of accuracy.
The reason ROMs were generated in this fashion was because no method had been defined to enable the computation of a static aeroelastic solution using a ROM. Any ROMs generated in this fashion were, therefore, limited to the prediction of dynamic responses about a static aeroelastic solution including the methods by Raveh 32 and by Kim et al. 20 The improved ROM method, however, includes a method for generating a ROM directly from a steady, rigid solution. As a result, these improved ROMs can then be used to predict both static aeroelastic and dynamic solutions for any dynamic pressure. In order to capture a specific range of aeroelastic effects (previously obtained by selecting a particular dynamic pressure), the improved ROM method relies on the excitation amplitude of the orthogonal functions to excite aeroelastic effects of interest. The details of the method for using a ROM for computing both static aeroelastic and dynamic solutions is presented in another reference by the author. 21 For the present results, all responses were computed from the restart of a steady, rigid FUN3D solution, bypassing the need (and the additional computational expense) to execute a static aeroelastic solution using FUN3D.

III.B. Error Minimization
Error minimization consists of error quantification and error reduction. Error quantification is defined as the difference (error) between the full FUN3D solution due to the orthogonal input functions used (Walsh) and the unsteady aerodynamic ROM solution due to the same orthogonal input functions. This was identified in Step 5 in the previous subsection and is shown schematically in Figure 5. The outputs shown are GAF responses per mode. Within the system identification algorithms, there are parameters that can then be used to reduce the error (error reduction). These parameters include number of states and the record length of the identified pulse responses, for example. The maximum error is the largest error encountered per mode. Using the maximum error as the figure of merit, the parameters are varied until an acceptable ROM has been obtained.

IV. Sample Results
A brief summary of results for three configurations is presented in this section. These configurations are the Lockheed-Martin N+2 low-boom supersonic configuration, the Royal Institute of Technology (KTH) generic fighter wind-tunnel model, and the AGARD 445.6 wing.

IV.A. Low-Boom N+2 Configuration
An artist's rendering of the Lockheed-Martin N+2 low-boom supersonic configuration is presented in Figure 6. This configuration has been used extensively as part of a NASA research effort to address the technologies required for a low-boom aircraft, including aeroelastic effects. Presented in Figure 7 is a comparison of the dynamic aeroelastic responses of the time histories of the first mode generalized displacements from a full FUN3D aeroelastic solution and the ROM aeroelastic solution at a Mach number of 1.7 and a dynamic pressure of 2.149 psi. Presented in Figure 8 is a comparison of the dynamic aeroelastic responses of the time histories of the second mode generalized displacements from a full FUN3D aeroelastic solution and the ROM aeroelastic solution at the same condition. As can be seen, the results indicate an excellent level of correlation between the full FUN3D solutions and the ROM solutions. Similar results are obtained for all the other modes, indicating good confidence in the ROM.   A major benefit of this ROM technology is the ability to rapidly generate an aeroelastic root locus plot that reveals the aeroelastic mechanisms occurring at that flight condition. Figure 9 presents the aeroelastic root locus plot for the low-boom N+2 configuration at M=1.70. This root locus plot clearly indicates the aeroelastic mechanisms that affect this configuration. In the root locus plot, each symbol represents the aeroelastic roots at a specific dynamic pressure. In this case, each increment in dynamic pressure corresponds to 2 psi. It is important to mention that this root locus plot is generated in seconds while multiple full FUN3D solutions would be required for each dynamic pressure of interest, with each solution requiring about two days.
The computational cost of generating these ROM solutions consists of one full FUN3D solution that is used to generate the ROM at that Mach number. This full FUN3D solution ran for three hours and consisted of 2400 time steps. Once this solution is available, a ROM can be generated and then used to generate all the aeroelastic responses at all dynamic pressures. In comparison, a full FUN3D analysis at each dynamic pressure requires two full FUN3D solutions: a static aeroelastic (∼10 hours) and a dynamic aeroelastic (∼18 hours). Therefore, full FUN3D solutions for 20 dynamic pressures would require ∼560 hours of compute time. Figure 9. Aeroelastic root locus plot for the low-boom N+2 configuration at M=1.7 with each colored marker indicating an increment of 2 psi in dynamic pressure for a given mode.

IV.B. KTH Generic Fighter
In 1985 and 1986, two wind-tunnel models of the Saab JAS 39 Gripen were designed, built, and tested in the NASA Transonic Dynamics Tunnel (TDT) for flutter clearance. One model, referred to as the stability model, was designed to be stiff, but incorporated proper scaling of both the mass and geometry. The other model, referred to as the flutter model, was also designed for proper scaling of structural dynamics, and was used for flutter testing with various external stores attached.
For the current collaboration, a generic fighter flutter-model version of these earlier models was selected. The new model, shown in Figure 10, has a similar outer mold line (OML) to the Gripen, but it has been modified into a more generic fighter configuration. Specifically, the air intakes were removed from the fuselage and the wing received an aspect ratio increase and a leading-edge sweep reduction. Details regarding the design, fabrication, and instrumentation of the wind-tunnel model can be found in the reference paper. 33 Figure 11 shows the wind-tunnel model installed in the TDT.
Using the AEROM software, aeroelastic root locus plots were generated for the KTH wind-tunnel model in air test medium for a free-air case and a solution accounting for the effects of the TDT test section via CFD modeling, 34, 35 as can be seen in Figure 12. There were three configurations tested: wing with tip stores (configuration 1), wing with tip and under-wing stores (configuration 2), and wing with tip and under-wing stores with added masses at tip stores (configuration 3). The third configuration exhibited flutter while configurations 1 and 2 did not. Presented in Figure 13 is the aeroelastic ROM root locus plot for the free-air configuration at M=0.90. For this case, the roots clearly indicate a flutter mechanism at about 8100 N/m 2 (or 169 psf) via a coalescence of modes 5 and 6. Using the ROM, any dynamic pressure can be quickly evaluated to determine the aeroelastic response, consistent with the root locus plots. At this dynamic pressure, the ROM-based flutter prediction is above the experimental flutter dynamic pressure at M=0.9 (not conservative). All results presented are for zero structural damping. Using the ROM, the effect of structural damping can be quickly evaluated as well but is not pursued in the present discussion.    Presented in Figures 14 and 15 are comparisons of the aeroelastic responses for modes 3, 4, 5, and 6 at M=0.9 and Q=7344 Pa for the FUN3D solution that includes the effect of the TDT and the ROM solution for the same configuration. As can be seen, the comparison is quite good with some variation in mode 6. Additional studies are currently underway to minimize these variations in order to reduce the error associated with the ROM.
For the CFD model that included the TDT, the ROM solution required two days whereas the full solution (for only one dynamic pressure) required five days. The ROM solution could, of course, then be used to rapidly compute the aeroelastic response due to any dynamic pressure.

IV.C. AGARD 445.6 Wing
Aeroelastic transients for the AGARD 445.6 aeroelastic wing 36 from the FUN3D full solution and from the aeroelastic ROM, for inviscid and viscous solutions, are presented in this section. The FUN3D full solution aeroelastic transients are presented at various dynamic pressures for two Mach numbers: M=0.96 and M=1.141. The aeroelastic transients generated using the FUN3D full solutions are used to validate the ROM results at specific dynamic pressures. Aeroelastic transients and root locus plots, generated using the aeroelastic ROM, are presented for the same Mach numbers.

IV.C.1. Inviscid Results
Inviscid FUN3D results are presented for both full FUN3D and ROM solutions. Figure 16 is the aeroelastic root locus plot for M=0.96 generated using the ROM method. The root locus plots consist of values at twenty dynamic pressures from zero to 114 psf. The ROM procedure generates a combined aeroelastic statespace model that consists of a state-space model of the structural dynamics and a state-space model of the unsteady aerodynamics (from FUN3D). Therefore, the root locus plots can be generated for any number, and any increment, of dynamic pressures rapidly. The aeroelastic root locus plot at this Mach number indicates a flutter mechanism dominated by the first mode with some coupling with the second mode while the third and fourth modes are stable.   A zoomed-in version of the root locus plot is presented as Figure 17. The increment in dynamic pressure for this root locus plot is 6 psf starting with 0 psf, resulting in a flutter dynamic pressure of approximately 30 psf. This flutter dynamic pressure is consistent with the FUN3D full solution flutter dynamic pressure. 36 However, at this Mach number, the inviscid result does not compare well with the experimental flutter dynamic pressure. This discrepancy is not surprising as inviscid solutions tend to have stronger shocks that are farther aft and, therefore, induce a stronger and earlier onset of flutter. With the inclusion of viscosity, the shock strength is reduced and the shock moves forward, yielding a higher flutter dynamic pressure.  Figure 18 presents the ROM aeroelastic root locus plot for M=1.141. There are two flutter mechanisms: the first flutter mechanism is an instability involving the first mode; the second flutter mechanism involves an instability of the third mode. The first mode instability occurs at a dynamic pressure of about 300 psf while the third mode is always unstable. The third mode instability was not expected as it is not mentioned in any of the references. In order to validate the accuracy of this aeroelastic root locus plot, the generalized coordinates from a FUN3D full solution are analyzed.
The aeroelastic transients for the four modes at M=1.141 and a dynamic pressure of 30 psf are presented in Figure 19. The first mode, with the largest amplitude, is clearly stable. However, the stability of the other three modes is not as obvious due to the smaller and similar amplitudes of these three modes. Figure 20 presents the generalized coordinate response of only the third mode where the unstable nature of this mode becomes obvious.   This third mode instability is not mentioned in any other publications on the flutter boundary of the AGARD 445.6 wing. Is this third mode instability present in all inviscid (Euler) solutions of the AGARD 445.6 wing presented in the literature? It appears that the first mode instability was the focus of all previous inviscid analyses at supersonic conditions. In that case, if evaluation of stability was based on a visual examination of the generalized coordinates, it is understandable how the third mode instability might have been missed. In addition, for analyses performed in the early days of computational aeroelasticity, Figure 19 would have consisted of fewer time steps (due to computational cost at the time), thereby making it difficult to visually notice the third mode instability. The authors have confirmed the existence of this third mode instability in previous solutions obtained using the NASA Langley CFL3D code as well.

IV.C.2. Viscous Results
Viscous FUN3D full and ROM solutions are presented at M=1.141 in this section. The results for FUN3D full and ROM viscous solutions at subsonic Mach numbers agree well with each other and with experiment and are not presented here. Figure 21 is the root locus plot generated using the FUN3D ROM viscous solution at M=1.141, in dynamic pressure increments of 6 psf to 114 psf. The third mode instability, noticed in the inviscid solution, has been stabilized by the inclusion of viscous effects. It is important to state a fundamental and important difference between a root locus plot and the visual, or otherwise post-processed analysis, of generalized coordinates over a short period of time. A root locus plot, by definition, exhibits the roots of a system as time approaches infinity or as the system reaches steady state. In contrast, the analysis of the initial transient response of a generalized coordinate over a short period of time can be deceiving as the response can change if the response was viewed (or analyzed) over a longer period of time. This property of root locus plots is critical for the accurate evaluation of aeroelastic stability.

V. Conclusions
The origin, implementation, and applications of AEROM, the patented NASA reduced-order modeling software, have been presented. Recent applications of the software to analyze complex configurations, including computation of the aeroelastic responses of the Lockheed-Martin low-boom N+2 configuration, the KTH (Royal Institute of Technology, Sweden) generic fighter wind-tunnel model, and the AGARD 445.6 wing were presented. Results presented demonstrate the computational efficiency and analytical capability of the AEROM software.