- Please check and comment entries here.
Aeroelasticity Methods in Turbomachinery
Aeroelastic phenomena in turbomachinery are one of the most challenging problems to model using computational fluid dynamics (CFD) due to their inherent nonlinear nature, the difficulties in simulating fluid–structure interactions and the considerable computational requirements. Nonetheless, accurate modelling of self-sustained flow-induced vibrations, known as flutter, has proved to be crucial in assessing stability boundaries and extending the operative life of turbomachinery. Flutter avoidance and control is becoming more relevant in compressors and fans due to a well-established trend towards lightweight and thinner designs that enhance aerodynamic efficiency.
2. Aeroelasticity Methods in Turbomachinery
The aim of aeroelastic methods for flutter prediction in turbomachinery is to solve the generic motion Equation:
where x is the vector containing physical coordinates of the structural model degrees of freedom (d.o.f.); M, D and K are the mass, damping and stiffness matrices; and Fa(t) is the time-dependent vector of aerodynamic forces acting on degrees of freedom at the domain interfaces. When different effects are included, such as nonlinear forces due to friction dampers, the related terms Fnl(t) are usually added to aerodynamic forces. Computational approaches to model flutter in turbomachinery could be divided into two broad categories: classical and integrated methods, the former ignoring the interaction between the fluid and the structure, and the latter attempting to model it.
2.1. Classical Methods
Classical aeroelasticity methods are those where the fluid and structural domains are uncoupled, so that the fluid flow does not affect the structural response which is pre-determined. Such methods thus split an inherently coupled nonlinear phenomenon into two separate uncoupled analyses.
One of the key concepts in turbomachinery aeroelasticity, introduced by Lane , is the interblade phase angle (IBPA): the individual blades in a cascade are assumed to vibrate with the same amplitude and mode shape, but the maximum is reached with a constant phase lag, known as IBPA. With non-zero IBPAs, the cascade is assumed to vibrate in the so-called travelling wave mode, whereas a zero IBPA implies a constant phase between all blade displacements and a standing wave pattern. Figure 2 provides a visual representation of different nodal diameters (ND) in a vibrating bladed disk, where colours indicate blade displacements. IBPA is related to the number of NDs (diameters of the disk where blade displacements are zero) according to the Equation
Figure 2. Vibration patterns with different NDs in a compressor bladed disk.
In classical methods, the structure and the fluid are decoupled so that a free vibration problem can be solved first. Normal mode-shapes and frequencies can either be calculated in a free vibration problem (in vacuum) or taking into account gyroscopic effects, centrifugal stiffening and steady airloads of a given operative condition. Then, the predicted mode-shapes are used with an arbitrary amplitude to calculate the unsteady aerodynamic forces, from which a measure of stability is inferred. Typical outputs of classical imposed-oscillation analysis include real and imaginary components of pressure distribution on blades’ surfaces, as shown in Figure 3. Different formulations of classical methods have attained high popularity thanks to three main factors:
Figure 3. Real and imaginary part of unsteady pressure on a 2D cascade using different CFD techniques, from Ekici .
- Simplicity of models and relatively low computational costs, due to lack of interexchanging information between the two domains at each substep.
- The assumption of uncoupling between the fluid and structural domain is usually realistic for many situations, especially at low speeds, so that the motion of the structure is well defined by a certain mode shape and frequency, which are substantially unaffected by the unsteady flow.
- The information regarding asymptotic stability of given modes and IBPAs, though conservative, is sufficient in many cases, especially for preliminary design purposes.
The main drawback of this method is the inherent simplification of flow nature given by the linearisation of unsteady flows, especially true when nonlinear phenomena such as shock waves are present and the mass ratio is low. A common way to partially circumvent this aspect is achieved by varying the prescribed (small) oscillation amplitude and assessing the linearity of flow response on the blades. Another direct consequence of uncoupling is the inability to predict limit cycle behaviour in blade oscillation patterns, which could prevent a complete understanding of the blade behaviour under actual loading conditions.
2.1.1. Actuator Disk Theory
Actuator disk representations are simplified forms of linearized cascade theories, based on the fundamental assumptions of low reduced frequency and small interblade phase angle (Whitehead ). Tanida and Okazaki  developed a variant where the hypothesis of small frequencies was relaxed, and the method was improved to be used with 2D compressible flows by Adamczyc . The actuator disk model calculated the unsteady aerodynamic forces as a function of the steady flow field and the dynamic response of the cascade. Linearized potential equations were generally adopted, and viscous and stalling effects were included via experimentally derived quasi-steady loss coefficients. In principle, different viscous loss coefficients should be derived for each new geometry; however, it was not uncommon to use the original coefficients for a wide range of blades. The main limitation is that the method requires ‘tuning’ to each particular type of problem, and its performance would be poorer where viscous effects are known to be important, such as in transonic stall flutter. The early version of actuator disk theory, as developed by Whitehead, could not properly simulate non-zero IBPA and used the assumption of infinitesimal chord length and spacing, so that flutter boundaries in strictly limited conditions could be computed. Later methods such as the one by Adamczyc adopted compressible transonic flows and quasi-steady total pressure loss coefficients to predict flutter boundaries of four high-speed rotors, with slightly conservative results. The steep variation in viscous forces associated with a rapid rise in the tip blade-element was identified as the cause of flutter onset.
2.1.2. Indicial Method
The indicial method was proposed by Ballhaus and Goorjian  and further explored by Ueda and Dowell  and Stark . It used transonic small disturbances (TSD) or linearized potential equations to compute the unsteady forces acting on an airfoil in transonic flow. The flow response due to a step-change in the airfoil position was computed first, and the response to any other prescribed motion was obtained via the convolution integral from this initial result. The cascade motion was assumed to be simple harmonic, so that the aeroelastic equations of motion could be solved by standard matrix techniques once the aerodynamic forcing terms were known. The indicial method by Ueda and Dowell proved suitable in predicting flutter boundaries of two-dimensional pitching cascades in a wide range of reduced frequencies with a very small amount of computation time compared to time-marching solutions. The method was capable of predicting nonlinear limit cycle flutter confirmed by time-marched simulations, but generally, the accuracy of the describing function method decreased as the amplitude of the motion increased. Moreover, the neglection of the component in the upwash distribution due to the angular velocity of the airfoil motion caused a fictitious flutter instability at high frequencies. However, the use of convolution integrals and excitation signals presents the advantage of predicting the flutter boundary over a range of conditions from only one set of time-marching fluid calculations. It is also the basis of another method named aerodynamic influence coefficients, frequently adopted as a ROM and described later on in the paper.
2.1.3. Aeroelastic Eigensolution Method
This method is based on obtaining the linearized harmonic unsteady aerodynamic coefficients for the motion of a freely vibrating structure. The unsteady aerodynamics can be provided by anything from empirically determined airfoil lift and moment coefficients, through linearized potential methods, to nonlinear RANS codes. The structural model is usually represented as a two-degree-of-freedom section, but the formulation is general, so that 3D blade descriptions can also be accommodated. Once determined, the aeroelastic forces are expressed in the frequency domain, either directly if analytical theories are used or by Fourier analysis if the forces were first calculated in the time domain. The resulting aeroelastic equations of motion are obtained by adding the aerodynamic contributions to the mass and/or stiffness matrices of the structural equation. However, these new system matrices may well become a function of the frequency, in which case, the eigenproblem is no longer linear. In such situations, an approximate solution can be found using iterative techniques. The stability of the system is then assessed by considering the amount of damping in each aeroelastic mode. Many applications of the methods can be found in the open literature, among those [25,26,27]. Early models employing the aeroelastic eigensolution method used potential linearized CFD codes, computing the airloads as a linear function of the blade displacement and assuming the superposition of effect due to different modes. The flutter boundaries obtained are essentially qualitative; nevertheless, the dependency of vibration eigenmodes on airloads and a resonance at nominal speed for an experimental fan was captured by Klose . More advanced techniques employing multiple-passage domains often rely on the assumption of airloads which are linearly dependent to blade displacement and velocity. The model of Tateishi  correctly predicted flutter boundaries in a range of conditions except low speed lines, where the effect of grid size and the accurate simulation of shock impingement appear to be critical factors.
2.1.4. Energy Method
The energy method, first developed by Carta , is based on calculating the sum of the work done by unsteady aerodynamic forces on blades, which oscillate with a prescribed motion of pitching or plunging (2D sections), or determined from a free-vibration analysis (when 3D models are used). It is inherently assumed that the flutter behaviour is linear: positive aerodynamic work, i.e., energy transfer to the structure, indicates instability, while a stable situation is characterized by a negative net aerodynamic work. The aerodynamic damping δ
is defined as the work of aerodynamic forces on blades over a cycle W, usually normalized with respect to the blade vibration mechanical energy Ke, as indicated in Equation (3).
The stability of the bladed disk is assessed by calculating the aerodynamic damping for each IBPA for every single mode-shape oscillating at its natural frequency, as shown in Figure 4, from Mahler . The method was historically implemented calculating the work per unit of a span on a number of 2D sections, using linearised aerodynamic theories. The process of assessing the overall stability by performing this worksum procedure was found to be ill-conditioned, the result being indicated by a small number, obtained as the algebraic sum of large numbers; the accuracy of the calculated damping was thus largely affected by the number of sections chosen. Today, with more computational power available, this approach and its many variants are implemented with 3D simulations and represent the state of the art in turbomachinery flutter analysis in the aerospace industry because of their simplicity. However, the fluid and structure domains are still treated as two separate media; thus, considerations and limitations of classical methods hold true. Examples of energy method applications are found in many works, including simulations of stall flutter and flows with large, separated regions, e.g., in . The method is also used to compare the accuracy of coupled methods to well-known uncoupled formulations, as in [26,30,31,32,33]. Figure 5, taken from Mahler , plots aerodynamic damping versus IBPA for two different coupling algorithms.
Figure 4. Aerodynamic damping over IBPAs for different modes.
Figure 5. Aerodynamic damping over IBPAs for a fixed mode using different aerodynamic codes.
2.2. Integrated Aeroelasticity Methods
Integrated aeroelasticity methods are those which do not uncouple the flow motion from that of the structure but attempt to treat the problem in one continuous medium . The need for such an approach arises from the inherent nonlinear behaviour of unsteady high-speed viscous fluid flows, especially when shock waves and separation phenomena are considered, e.g., in transonic flutter or in stalled conditions. The existence of nonlinear structural effects such as friction damping in the blade roots and shroud interfaces further justifies the implementation of such methods. The mathematical formulation of integrated methods is devised in such a way that the flow behaviour is allowed to modify the structural motion, and vice versa, as in the true physical phenomenon. The most striking difference between the classical methods and integrated ones is that the former can only predict the onset of flutter as a sudden change from a stable to an unstable region, while the latter are capable of predicting limit-cycle behaviour, a constant amplitude oscillation due to nonlinear aerodynamic phenomena, as shown in Figure 6. Experimental evidence suggests that flutter can occur at different limit cycle amplitudes influenced by both aerodynamics and structural damping, and the prediction of that amplitude is crucial in assessing the life of turbomachinery components. Integrated algorithms can also predict nonlinear instabilities, i.e., situations where the equilibrium point is locally stable, but for sufficiently strong perturbation, the frictional dissipation cannot bound the flutter vibrations, as proved by Berthold . As with classical flutter analysis methods, different formulations are possible for the integrated methods: one key factor in categorizing different algorithms is the numerical treatment of the physical domains (i.e., whether fluid and structure are discretized into a single domain or not). Another parameter is the accuracy in assessing fluid–structure equilibrium at each physical timestep in time-marching simulations.
Figure 6. Total energy of a cascade and deflection angle at three different conditions: stable mode, limit cycle behaviour and unstable mode, from Sadeghi .
2.2.1. Partially Integrated Method
In the partially integrated methods, the solutions of the fluid and structural equations are calculated separately, but the information is exchanged at each time step, so that the solution from one domain is used as a boundary condition for the other domain. In other words, a new blade position is calculated at each time step using the fluid forces of the previous time step, and this new position is used as the new fluid–structure boundary for the next time step when the aerodynamic forces are computed first. Many examples of partially integrated methods can be found in the open literature. In the works by Sayma , Sadeghi , Debrabandere , Zheng  and Li , the eigenmodes were computed using constant or Rayleigh structural damping model, and the domain was reduced to one or two passages to contain the computational effort, with a limitation on the explorable IBPA range as a payoff. Coupled simulations showed that the normal ’in-vacuum’ eigenmodes can be significantly different from the aeroelastic mode, especially at low mass ratios, leading to inaccurate results when assumed modes and the energy method are employed, as evidenced in Figure 7. Furthermore, a nonlinear fluid–structure interaction can lead to the oscillation of blades around a displaced position, with two different harmonics and a pattern of alternated chocked and unchocked flows in neighbouring passages, as evidenced by Sadeghi . Fluid–structure coupling also permits the simulation of a forced response ; the influence of the passing rotor on stator blades is captured, and the Fourier transformation of induced vibrations showed that they are more related to the natural stator frequencies than to the rotor passing frequency. Various works adopted multiple mode superposition in coupled simulations, but some evidence suggested that the influence of modes next to the first is often negligible. In the papers by Liang  and Zhang , the structural reduction to a few modes (including centrifugal stiffening and the Coriolis effect) was coupled to the three-dimensional full annulus fluid model of transonic compressors, relaxing the limitations on permitted IBPA and allowing the simulation of complex phenomena such as rotating stall and nonsynchronized vibrations. The computational model of Liang  showed that rotating stall cells have about 30% the rotor speed (see Figure 8) and an energy mainly distributed in the first two orders of excitation. Some resonances were found between blade modes and the stalled flow cell velocity; a novel algorithm was built to predict those resonances. In the paper by Zhang , the rotating stall was simulated showing a limit-cycle oscillation of all blades as the stall line was approached, followed by flutter as the system became dynamically unstable, with all the blades vibrating at the same mode and frequency but with different amplitudes. A sinusoidal travelling wave pulsation was observed during flutter, but with a non-constant IBPA. Im  simulated flutter for the NASA Rotor 67 at different operative conditions using a fully coupled partially integrated model with both a four-passage domain and the full annulus, showing comparable results for the two cases in unstalled flutter conditions.
Figure 7. Reduced frequency at neutral stability over mass ratio. From Sadeghi .
Figure 8. Axial velocity contours of a bladed disk during rotating stall. From Liang .
2.2.2. Periodic Mode Updating Method
This method is a hybrid between classical and partially integrated aeroelasticity methods. Similar to the aeroelastic eigensolution method, a free vibration problem is solved in the frequency domain, and the mode shape of interest is used to describe the blade motion for one period using a 3D time-marching nonlinear CFD code. After the first period is simulated, the mode shape is updated using aerodynamic coefficients which were calculated in the previous period, via a frequency-domain structural calculation based on modal projection. The cycle of time-domain aerodynamics, followed by frequency-domain structural dynamics, is repeated until the mode calculated in the frequency domain converges. An example of such a procedure can be found in the work of Gerolymos . The model adopted a rather coarse mesh and was used for illustrative purposes. It was found that the coupling slightly affected the mean accumulated power on blades, but this difference was sufficient to lower the aerodynamic damping by 30%.
2.2.3. Fluid–Structure Coupling Algorithms
Since each of the components of the coupled aeroelastic problem has different mathematical and numerical properties, well-established but distinct numerical solvers and readily available commercial software, the simultaneous solution of the equations by a monolithic scheme is, in general, computationally challenging, mathematically and economically suboptimal and unmanageable software-wise. On the other hand, the solution of the aeromechanical system through different integration schemes and the exchange of information between the two domains presents various issues and requires some precautions.
The simplest form of the integration algorithm is depicted in Figure 9. The flow solver is advanced for one physical time step, then the fluid load is sent to the structural solver that computes the corresponding deformation and a new iteration starts. As the equilibrium between the structure and the fluid is not ensured at the end of each iteration, this is called a weak-coupling method. Examples of weakly coupled algorithms can be found in [34,35]. More accurate representations of fluid–structure interactions ensure the equilibrium of quantities such as force, displacement and energy at each physical substep by subiterating the governing equations until the prescribed tolerance on the chosen quantities is reached: those methods are named fully coupled and were adopted in [10,31,37,38,39].
Figure 9. Conventional serial staggered algorithm. The circled numbers indicate steps order.
Since the typical time integration steps needed for accuracy and stability reasons may differ by orders of magnitude between the structural and fluid problem, in almost all aeroelastic problems, it is the integration of the fluid state vector that requires a much finer temporal resolution than the structural vibration. For this reason, some algorithms were developed that subiterate the fluid flow with a factor equal to the structural/fluid timestep ratio, reducing the needed advancements of the structural domain and I/O transfers, as shown in Figure 10. However, in certain cases, the cost of subiterations can offset the benefits of the larger timestep they can allow.
Figure 10. Conventional serial staggered algorithm with fluid subiterations.
More advanced methods try to achieve an accuracy comparable to that of a monolithic or fully coupled scheme but with a reduced computational cost: among these are the displacement and velocity conserving staggered algorithm by Farhat and Lesoinne  and employed by Carstens . As schematized in Figure 11, the computations of structural displacements and fluid state vectors are shifted by t/2, while the CFD-mesh deflections at the blades’ surfaces are calculated from the structural displacements. Furthermore, the mesh deformation is computed by Newmark’s algorithm estimating the constant acceleration in the time interval [t,t+Δt] from the unsteady aerodynamic forces at time t+t/2.
Figure 11. Improved serial staggered procedure.
2.3. Fully Integrated Methods
Fully integrated methods differ from their partially integrated counterparts in the discretization and subsequent numerical treatment of the governing equations. They are based on formulating the structural and fluid dynamics together so that they can be solved at each time step using the same integrator. The two domains are discretized into one arbitrary Lagrangian–Eulerian (ALE) space; as a consequence, the motion of the grid becomes an integral part of the equations of motion and does not have to be handled separately. The direct method uses an explicit temporal discretization which is integrated using a Runge–Kutta scheme, and upwind differencing is used for the spatial discretization of the arbitrary Lagrangian–Eulerian formulation of the aeroelasticity equations. The structural equations are formulated on a local nodal level, which enables them to be discretized using the same integrator as the aerodynamics. Thus, both sets of equations are handled simultaneously at the same time step without distinguishing between the fluid and the structure. This model is claimed to calculate the energy transfer between the structure and the fluid more accurately than similar partially integrated schemes since there is no possibility of a time lag between the structure and the fluid. The method was developed by Bendiksen  and applied to both wing and cascade flutter using quasi three-dimensional models, predicting divergent behaviours and limiting cycle flutter in a range of conditions.
The entry is from 10.3390/aerospace8090242
- Collar, A.R. The expanding domain of aeroelasticity. J. R. Aeronaut. Soc. 1947, 51, 1–34.
- Marshall, J.G.; Imregun, M. A review of aeroelasticity methods with emphasis on turbomachinery. J. Fluids Struct. 1996, 10, 237–267.
- Lane, F. System mode shapes in the flutter of compressor blade rows. J. Aeronaut. Sci. 1956, 23, 54–66.
- Ekici, K.; Hall, K.C. Time-linearized Navier—Stokes analysis of flutter in multistage turbomachines. In Proceedings of the 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 10–13 January 2005; AIAA: Reston, VA, USA; Volume 836, pp. 1–19.
- Whitehead, D.S. The vibration of cascade blades treated by actuator disk methods. Proc. Inst. Mech. Eng. 1959, 173, 555–557.
- Tanida, Y.; Okazaki, T. Translatory vibration of cascade blades as treated by semi-actuator disk methods—Parts 1 and 2. Trans. JSME 1963, 6, 744–758.
- Adamczyk, J.J. Analysis of Supersonic Stall Bending Flutter in Axial-Flow Compressor by Actuator Disk Theory; NASA Technical Paper No. NASA-TP-1345; NASA: Washington, DC, USA, 1978.
- Ballhaus, W.F.; Goorjian, P.M. Computation of unsteady transonic flows by the indicial method. AIAA J. 1978, 16, 117–124.
- Ueda, T.; Dowell, E.H. Flutter analysis using nonlinear aerodynamic forces. AIAA J. Aircr. 1984, 21, 101–109.
- Stark, V.J.E. General equations of motion for an elastic wing and methods of solution. AIAA J. 1984, 22, 1146–1153.
- Klose, A.H. Advanced Ducted Engines: Impact of Unsteady Aerodynamics on Fan Vibration Properties; ASME: New York, NY, USA, 1992.
- Mahler, A.; Placzec, A. Efficient coupling strategies for the numerical prediction of the aeroelastic stability of bladed disks. In Proceedings of the International Symposium on Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines, Tokyo, Japan, 8–11 September 2012.
- Tateishi, A.; Watanabe, T.; Himeno, T.; Aotsuka, M.; Murooka, T. Verification and Application of Fluid-Structure Interaction and a Modal Identification Technique to Cascade Flutter Simulations. Int. J. Gas Turbine Propuls. Power Syst. 2016, 8, 20–28.
- Carta, F.O. Coupled blade-disc-shroud flutter instabilities in turbojet engine rotors. ASME J. Eng. Power 1967, 89, 419–426.
- Clark, W.S.; Hall, K.C. A time-linearized Navier-Stokes analysis of flutter. In Proceedings of the ASME 1999 International Gas Turbine and Aeroengine Congress and Exhibition, Indianapolis, IN, USA, 7–10 June 1999.
- Nowinski, M.; Panowsky, J. Flutter Mechanisms in Low Pressure Turbine Blades. Trans. ASME J. Eng. Gas Turbines Power 2000, 122, 82–88.
- Sadeghi, M.; Liu, F. Computation of cascade flutter by uncoupled and coupled methods. Int. J. Comput. Fluid Dyn. 2005, 19, 559–569.
- Su, D.; Zhang, W.; Ye, Z. A reduced order model for uncoupled and coupled cascade flutter analysis. J. Fluids Struct. 2016, 61, 410–430.
- Lian, B.; Hu, P.; Chen, Y.; Zhu, X.; Du, Z. Insight on aerodynamic damping of the civil transonic fan blade. In Proceedings of the ASME Turbo Expo 2021 Turbomachinery Technical Conference and Exposition GT2021, Virtual, 7–11 June 2021.
- Berthold, C.; Gross, J.; Frey, C.; Krack, M. Development of a fully-coupled harmonic balance method and a refined energy method for the computation of flutter-induced Limit Cycle Oscillations of bladed disks with nonlinear friction contacts. J. Fluids Struct. 2021, 102, 103233.
- Sayma, A.I.; Vahdati, M.; Imregun, M. An integrated nonlinear approach for turbomachinery forced response prediction. Part I: Formulation. J. Fluids Struct. 2000, 14, 87–101.
- Debrabandere, F.; Tartinville, B.; Hirsch, C.; Coussement, G. Fluid—Structure interaction using a modal approach. J. Turbomachinery-Trans. ASME 2012, 134, 051043.
- Zheng, Y.; Yang, H. Coupled Fluid-structure Flutter Analysis of a Transonic Fan. Chin. J. Aeronaut. 2011, 24, 258–264.
- Li, J.; Yang, X.; Hou, A.; Chen, Y.; Li, M. Aerodynamic Damping Prediction for Turbomachinery Based on Fluid-Structure Interaction with Modal Excitation. Appl. Sci. 2019, 9, 4411.
- Liang, F.; Xie, Z.; Xia, A.; Zhou, M. Aeroelastic simulation of the first 1.5-stage aeroengine fan at rotating stall. Chin. J. Aeronaut. 2020, 33, 529–549.
- Zhang, M.; Hou, A.; Zhou, S.; Yang, X. Analysis on Flutter Characteristics of Transonic Compressor Blade Row by a Fluid-Structure Coupled Method. In Proceedings of the ASME Turbo Expo 2012 GT2012, Copenhagen, Denmark, 11–15 June 2012.
- Im, H.; Chen, X.; Zha, G. Detached Eddy Simulation of Transonic Rotor Stall Flutter Using a Fully Coupled Fluid-Structure Interaction. In Proceedings of the ASME Turbo Expo 2011 GT2011, Vancouver, BC, Canada, 6–10 June 2011.
- Gerolymos, G.A. Coupled 3-D Aeroelastic Stability Analysis of Bladed Disks; ASME paper 92-GT-171; ASME: New York, NY, USA, 1992.
- Farhat, C.; Lesoinne, M. Higher-Order Staggered and Subiteration Free Algorithms for Coupled Dynamic Aeroelasticity Problems. In Proceedings of the 36th Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 12–15 January 1998.
- Carstens, V.; Kemme, R.; Schmitt, S. Coupled simulation of flow-structure interaction in turbomachinery. Aerosp. Sci. Technol. 2003, 7, 298–306.
- Bendiksen, O.O. A New approach to computational aeroelasticity. In Proceedings of the 31st Structural Dynamics and Materials Conference, Long Beach, CA, USA, 2–4 April 1991. AIAA Paper 91-0939.