UFDC Home  myUFDC Home  Help 



Full Text  
COMPUTATIONAL AEROELASTICITY USING A PRESSUREBASED SOLVER By RAMJI KAMAKOTI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by Ramji Kamakoti This dissertation is dedicated to my parents and sister. ACKNOWLEDGMENTS I would like to express my sincere gratitude to Professor Wei Shyy for his constant support and guidance throughout this work. Equally, I would like to thank Dr. Bhavani Sankar, Dr. Andrew Kurdila, Dr. Renwei Mei, Dr. Nagaraj Arakere, and Dr. Michael Frank for serving on my committee and providing their support in completing this work. I would like to extend my sincere gratitude to Dr. Siddarth Thakur and other members of the computational thermofluids laboratory for making the work environment very lively and enjoyable to work in. Lastly, I would like to acknowledge the support given by my family throughout my career. TABLE OF CONTENTS page A C K N O W LE D G M EN T S ............................... ......... ................................................... iv LIST OF TA BLES ........... ... ......................................... .................. viii LIST OF FIGURES ......... .................................. ............ ........... ix A B ST R A C T ...........................................................................................xiii CHAPTER 1 INTRODUCTION .................... .................... ....... ............. Aeroelasticity and the FluidStructure Interaction Problem........................................1 P problem Statem ent ................................................................... .. ................ .4 2 L ITE R A TU R E R E V IE W ...................................................................... ..................10 Aerodynamic Models ............. .................................10 Physical M odels................................................... 10 R educedO rder M odels ..................... ........................... ............... ....12 Review of Coupled Computational Aeroelasticity (CAE) Models .......................... 13 F u lly cou pled A n aly sis ............................................................. ................. .. 14 Loosely and Closely Coupled Analysis.................................... ............... 16 L oosely coupled analysis ........................................ ......... ............... 16 Closely coupled analysis ................................................... ....... ........ 17 Review of M oving Boundary M odels .................................. ............ .................. 22 Review of Geometric Conservation Law ....................................... ............... 24 R eview of Interfacing Techniques........................................................... ... .......... 26 3 GOVERNING EQUATIONS AND OVERVIEW OF ALGORITHM.........................32 G ov ern in g E qu action s ....................................................................... ....................32 F lo w M o d u le ................................................................................................. 3 2 NavierStokes equations.................................................. 32 Transformation to curvilinear coordinates ....................................... 33 Geometric Conservation Law ..............................................37 Turbulence M odeling ................................................ .............................. 38 The ketransport equations ........................................ ....... ............... 40 Filterbased turbulence model for unsteady ReynoldsAveraged Navier Stokes (RAN S) computations..... .......... ........................................42 B oundary conditions .............................................................................. 43 W all treatm ent ............................................................................. 43 Structural D ynam ics M odel ...................................................... .... ........... 44 M oving G rid M odule................................................. .............................. 48 O overview of A lgorithm .............................................................................. .............50 D iscretized Form of E equations ............... ..... ........... ......... .....................50 Evaluation of Contravariant velocities on Nonstaggered Grid ........................52 PressureBased Flow Solver (SemiImplicit Method for PressureLinked E qu ation s, SIM P L E ) ................ ...................... .................. ......................55 PressureImplicit Splitting of Operators (PISO) Algorithm for unsteady com putations ............................... ............................... ............... 58 Updating Jacobian values for moving boundary treatment..............................60 First order Im plicit Schem e:...................................................... .............61 Firstorder timeaveraged scheme: ............................. ................... 62 Second order im plicit schem e .............. ......................... ............... .... 62 Second order timeaveraged evaluation of Jacobian............... .......... 63 Newmark Integration Method for Structure Solver........................................64 4 COMPUTATIONAL PROCEDURE AND CODE VALIDATION.............................66 C om putational P procedure .................................................................. .....................66 Geometry definition and Computational Grids .................................. ...............67 G eom etry D definition ................................................................. ....................67 C om putational G rids ...................... ........ ................................ ................. 68 Computational fluid dynamic (CFD) grid............................... ...............68 Computational structural dynamic (CSD) grid ........................................69 Coupling and Interfacing Procedure.................................. .............................. ........ 70 C o d e V alid atio n ...................................... ............................. ................ 7 4 Steadystate CFD Com putations ............................................... .....................75 Unsteady Computations using PISO Algorithm ................................................77 Effect of number of stages on accuracy and stability of PISO algorithm ....80 Momentum Interpolation Techniques for Computing Contravariant Velocities 84 G eom etric C conservation L aw ......................................................... ................. .88 Twodimensional channel flow: First order backward Euler.......................89 Twodimensional channel flow: PISO algorithm.......................................94 Threedimensional elastic wing: AGARD 445.6 ........................................95 M oving Boundary M odule ................................................... ............... ... 100 Structure Solver .................. ....................................... .. ............ 102 5 RESULTS AND DISCU SSION ........................................................ ............... 105 Coupled Simulation for Incompressible Flow Conditions ......................................105 Comparison of PISO and SIMPLE Algorithms...................................................111 Coupled Simulation for Compressible Flow Conditions .............. .... ...............112 Time Scales and Choice of Time Step Size for the Coupled Problem ............13 Flutter Boundary Prediction for AGARD Wing at a Transonic Mach Numberl 16 Flutter Computations Using a FilterBased Turbulence Model (M=0.96)........124 Summary of Flutter Boundary Prediction for AGARD Wing.......................... 128 6 CONCLUSIONS AND FUTURE WORK .............. .............................................132 C o n clu sio n s................................................... .................. 13 2 Future D directions ............................................ ........ ........ .......... .. 137 LIST OF REFEREN CES ........................................................... .. ............... 138 BIOGRAPHICAL SKETCH ............................................................. ............... 144 LIST OF TABLES Table pag 21. Description and key results of a few fullycoupled analysis methods .....................15 22. Description of CAE simulations using CAPTSD, ENS3DAE and CFL3DAE ........19 23. Summary of work with a moving mesh algorithm .......................... .....................21 24. Summary of work related to ALE formulation .......................................................22 25. Comparison of moving mesh algorithms ..................................................24 26. Summary of representative interface techniques..................................................28 27. Summary of boundary element methods........ ....... ..... ............... ............... 30 41. Error Norm versus grid velocity for the four GCL schemes for 3D wing case using Backward Euler m ethod .................................. ............. ....................95 42. Error Norm versus grid velocity for the four GCL schemes for 3D wing case ........98 43. Tip deflection at two different time instants for different GCL schemes for 3D w ing case ........ ............................................................................................. 99 44. Comparison of wing mode shapes between 10 element beam model (present study) and 120 element plate model ............ ........ ............... .................102 51. Comparison of critical flutter speed and dynamic pressure with experiment and other num erical results ...................... .. .... ...........................................130 LIST OF FIGURES Figureage 11. Aeroelastic diagram of forces and associated phenomena............... ..................2 12. Flutter speed index prediction for AGARD 445.6 wing using several methods..........7 21. Sample MDICE environment for aeroelastic simulation .......................................17 22. Coupled fluidstructure flow diagram ............................................. ............... 27 23. Varying levels of complexity in modeling for fluids and structures ..........................27 31. Displacements Measured with respect to the Elastic Axis....................................46 32. Location of variables u, v and p on a 2D nonstaggered grid for the pressure based algorithm ................................. .... ...... ......... .. .............. .. 50 33. Overview of the SIMPLE algorithm ........................................... .................58 41. Schematic of the AGARD 445.6 wing used in the wind tunnel.............................67 42. Overview of the M ultiblock CFD grid............................. .................................... 69 43. CFD surface grid along with grid distributions at the leading and trailing edges......69 44. Schematic of the FEM grid on the AGARD wing............................... 70 45. Schematic to demonstrate interpolation technique..................................................71 46. Schematic of a super element: Portion of the entire structure............................... 72 47. Sample CFD mesh superimposed on the discretized beam structure.........................73 48. Schematic to demonstrate the extrapolation procedure............................................74 49. Top view of the CFD domain showing the type of boundary conditions specified at different surfaces ................................. ..... ..........................75 410. Steady state surface pressure contours on the AGARD wing ................................76 411. Steady state pressure coefficient distribution at different spanwise locations on the top surface ........................................................................76 412. Computational domain for flow past square cylinder along with imposed boundary conditions ............................................ ................. .. ...... 77 413. Periodic oscillation of the crossstream (v) component of velocity using PISO algorithm for square cylinder case at Re=215........ ....... ..... ......................... 79 414. Vordex shedding past a square cylinder using PISO algorithm for Re=215. A) A t=0.001, B ) At=0.0005 ............................................ .. ......... ........... 79 415. Periodic oscillation of the crossstream (v) component of velocity using SIMPLE algorithm for square cylinder case at Re=215 ..........................................80 416. Pressure residual history for unsteady flow over a square cylinder (Re=215).........81 417. Periodic oscillation of Crossstream velocity (v) using different number of stages for PISO algorithm ............................................... ............................. 82 418. Computational domain and boundary conditions imposed for flow over a circu lar cy lin d er............................. .................................................. ............... 83 419. Pressure residual history for unsteady flow over a circular cylinder (Re=100).......83 420. Periodic oscillation of crossstream velocity (v) for different number of corrector stag es ................................................................................... 84 421. Schematic of Cavity flow grid along with boundary conditions.............................85 422. Velocity and pressure contours for cavity flow at Re=100 using different momentum interpolation schemes for various time step sizes at y=0.5 location in th e cav ity ....................................................... ................ 86 423. Schematic of computational domain surrounding a cylinder ................................87 424. Velocity and pressure plot for flow around a cylinder at Re=40 using different momentum interpolation schemes for various time step sizes at the symmetry line dow stream of the cylinder.................................................................. 88 425. Computational grids for channel flow at different time instants............. ..............90 426. Velocity profile for channel flow with Re=100 at different time instants for coarse grid (151 x 11) using Backward Euler method .................... ..................91 427. Error norm versus grid velocity using various schemes for channel flow for 151x11 grid using Backward Euler method ................................. ................91 428. Velocity profile for channel flow with Re=100 at different time instants for fine grid (301 x21) using Backward Euler method ............. ........................................93 429. Error norm versus grid velocity using various schemes for channel flow for 301x21 grid using Backward Euler method ................................. ................93 430. Velocity profile for channel flow at different time instants for 151x11 grid using P ISO m ethod .........................................................................95 431. Plot depicting the arbitrary movement of the wing in the spanwise direction........97 432. Error norm versus grid velocity for the various schemes for AGARD wing using Backward Euler m ethod .................................... .... ................................. 97 433. Spanwise deflection of AGARD wing at four different time instants...................99 434. Schematic of multiblock grid used to validate moving mesh module ................100 435. Effect of the 2 parameters, FACMIN and 3, on the remeshing...........................101 436. Tip deflection of AGARD wing versus number of time steps for At=0.0001........104 437. Tip deflection of AGARD wing versus number of time steps for At=0.001..........104 51. Spanwise wing shapes at different time instants (Grid configuration I) ................106 52. Time varying displacement of wing at different spanwise locations (Grid con figu ration I) .................................................................... 10 7 53. Time history of lift coefficient for AGARD 445.6 wing subject to 1degree angle of attack for both grid configurations ......................................... ...............108 54. Time history of lift/drag ratio for AGARD 445.6 wing subject to 1degree angle of attack for both grid configurations .... ........... ......................................... 108 55. Pressure contour on the surface of the wing at steady state ................................... 109 56. Comparison of lift coefficient time history for AGARD wing subject to different angles of attack for grid configuration I............... ............................................. 110 57. Comparison of lift coefficient time history for AGARD wing subject to different angles of attack for grid configuration I............... ............................................. 110 58. Spanwise displacements at three different time instants to compare PISO and SIMPLE methods using incompressible flow around an AGARD wing example at a Re=3.64x105 based on unit root chord. ......... ................................ .......112 59. Diffusive and convective time scales near wing tip region for different time step sizes and grids .............. ..... .... ............................. ....... 114 510. Diffusive and convective nondimensional paramter at wing tip spanwise location for different grids and time step sizes ........ ........................................ 116 511. Generalized displacement versus time for three different dynamic pressures for At=5x 105 .................................... ........................... .... ..... .......... 118 512. Generalized displacement versus time for three different dynamic pressures for A t= lx 0 5 ........................................................................................ 1 19 513. Effect of grid resolution on generalized displacements using similar CFL n u m b e rs ................................ ............... ................................. .. 12 0 514. Flutter points for different choice of time step sizes (1x105 and 5x105 for grid I (350 K points) and grid configurations (350 K grid and 800 K grid) ..................121 515. Tip deflection of wing versus time for grid I (350 K points) using a time step size of 5x105 ................................... ......................... 122 516. Wing shapes at maximum and minimum tip deflection points ...........................123 517. Mach number contours representing supersonic region in the flow domain at m idspan plane. .................................................................... 123 518. Surface pressure contours indicating supercritical region or region of supersonic flow on the surface of the wing. ................................. ........................ ......... 124 519. Blending function plot at midspan plane................................... ............... 126 520. Comparison of filterbased model and standard ke model for q/qe=0.98 with At=5x105 and filter size, A=0.15. ........................................ ....... ............... 127 521. Flutter boundary comparison between filterbased turbulence model and standard kemodel using grid I and At=5x105 ................................. ............... 127 522. Spanwise wing shape at maximum and minimum tip deflection for (left) M =0.678 and (right) M =1.072 ........................................ ......................... 129 523. Mach number contours at maximum and minimum tip deflection points............129 524. Summary of flutter speed index prediction for AGARD 445.6 wing ....................131 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL AEROELASTICITY USING A PRESSUREBASED SOLVER By Ramji Kamakoti August 2004 Chair: Wei Shyy Major Department: Mechanical and Aerospace Engineering A computational methodology for performing fluidstructure interaction computations for threedimensional elastic wing geometries is presented. The flow solver used is based on an unsteady ReynoldsAveraged NavierStokes (RANS) model. A well validated ke turbulence model with wall function treatment for near wall region was used to perform turbulent flow calculations. Relative merits of alternative flow solvers were investigated. The predictorcorrectorbased Pressure Implicit Splitting of Operators (PISO) algorithm was found to be computationally economic for unsteady flow computations. Wing structure was modeled using BernoulliEuler beam theory. A fully implicit timemarching scheme (using the Newmark integration method) was used to integrate the equations of motion for structure. Bilinear interpolation and linear extrapolation techniques were used to transfer necessary information between fluid and structure solvers. Geometry deformation was accounted for by using a moving boundary module. The moving grid capability was based on a master/slave concept and transfinite interpolation techniques. Since computations were performed on a moving mesh system, the geometric conservation law must be preserved. This is achieved by appropriately evaluating the Jacobian values associated with each cell. Accurate computation of contravariant velocities for unsteady flows using the momentum interpolation method on collocated, curvilinear grids was also addressed. Flutter computations were performed for the AGARD 445.6 wing at subsonic, transonic and supersonic Mach numbers. Unsteady computations were performed at various dynamic pressures to predict the flutter boundary. Results showed favorable agreement of experiment and previous numerical results. The computational methodology exhibited capabilities to predict both qualitative and quantitative features of aeroelasticity. CHAPTER 1 INTRODUCTION The term computational aeroelasticity (CAE) generally refers to coupling high level computational fluid dynamic (CFD) methods with structural dynamic tools to perform aeroelastic analysis. Recently, CAE has gained interest as considerable progress has been made in CFD, computational structural dynamics (CSD), and in computer technologies. Extensive research on CFD and CSD has already been done. The aim of our study was to develop a closely coupled CAE model (comprising a detailed CFD model with a simplified CSD model) to perform fluidstructure interaction computations on threedimensional wing bodies. Aeroelasticity and the FluidStructure Interaction Problem Aeroelasticity can be defined as the phenomena associated with the interaction of aerodynamic forces and inertial forces within elastic structural systems. There are also aeroelastic phenomena associated with interaction between aerodynamic and elastic forces alone (Bisplingoff et al. 1955). Aeroelastic problems mainly arise from the flexible nature of the structure. In other words, rigid structures do not experience aeroelasticity of any sort. It is well known that external forces acting on a flexible structural system (such as a wing) lead to a deformation in the wing geometry, and this structural deformation thereby leads to additional aerodynamic loads. Figure 11. Aeroelastic diagram of forces and associated phenomena Generally, aeroelasticity has two major classes: dynamic and static. Dynamic aeroelasticity usually involves interactions among inertial, aerodynamic and elastic forces; whereas static aeroelasticity involves interaction between aerodynamic and elastic forces. Figure 1 shows different kinds of aeroelastic phenomena depending on how the different forces interact. Dynamic aeroelasticity: includes phenomena such as flutter, buffeting, and dynamic response. Flutter is an oscillatory dynamic instability (primarily caused by the elasticity of the structure) that occurs in an aircraft in flight at high speeds. The speed at which this occurs is called flutter speed or critical speed. Buffeting is a phenomenon that occurs because of transient vibrations of aircraft structural components (due to aerodynamic impulses such as wake behind wings). Dynamic response includes transient response associated with aircraft components; and is caused by rapidly applied loads (such as gusts, moving shock waves, or other dynamic loads). All of the abovementioned phenomena are all transient phenomena; hence the term dynamic aeroelasticity. Static aeroelasticiy: includes phenomena such as load distribution, divergence, and system reversal. Load distribution occurs because static deformation influences the Dynamic Aeroelasticity * Flutter * Buffeting * Dynamic response distribution of aerodynamic pressures over the structure. Divergence is another static instability (it occurs at a speed called divergence speed) in which the elasticity of the lifting surface plays a critical role in producing the instability. Another static aeroelastic phenomenon is control system reversal (it occurs at control reversal speed) in which the effects of structure displacements are cancelled by elastic deformations of the structure itself. Almost every flight vehicle (manned or unmanned) that flies through the atmosphere undergoes some degree of aeroelasticity. Catastrophic phenomena such as flutter must be avoided at all costs, and all vehicles must be cleared of such phenomena before they are put to use. Flight test and windtunnel testing are two ways to test for such phenomena, but they are both expensive and occur late in the design process. Hence, computational techniques are used first, to assess the aeroelastic characteristics of these flight vehicles. While computational methods that study different aspects of aeroelastic response have been studied for some time, numerous open research issues remain to be resolved. For example, many approaches in computational aeroelasticity seek to synthesize independent computational approaches for the aerodynamic and the structural dynamic subsystems. This strategy is known to be fraught with complications associated with the interaction between the two simulation modules. Some of the issues arise from the fact that CFD and CSD mesh systems are quite different. Frequently, the former uses a Eulerian or spatially fixedcoordinate system, while the latter uses a Lagrangian or material fixedcoordinate system. Hence, care must be taken to develop a suitable interfacing technique between the two modules. Also, the time scales can be very different for the two modules, hence one must be careful while performing unsteady calculations. There are three major classifications for CAE: fully coupled, closely coupled, and loosely coupled analyses. In loosely coupled analysis, the fluid and structure modules are treated as two separate modules, with only external interaction between them. This kind of methodology can be seen as a multidisciplinary problem. This method is limited to small perturbations with moderate linearity. In fully coupled analysis, the governing equations for fluids and structures are combined into one set of equations, and these equations are solved and integrated simultaneously. Since the matrices associated with structures are orders of magnitude stiffer than those associated with fluids, it is virtually impossible to solve the entire system using a single numerical scheme. Methods have been developed using fully coupled methods, but they are restricted to twodimensional problems and smallscale threedimensional problems. In the closely coupled approach, fluids and structures are modeled in separate domains, but are coupled into one module by an interface technique. The exchange of information between these modules takes place at the interface or the boundary. The coupling is integrated, thereby allowing the two modules to exchange information at the boundaries in an efficient manner. Our study emphasizes this kind of approach. Problem Statement The objective of our study was to develop a computational model that is capable of performing fluidstructure interaction computations on threedimensional geometries. Our model was based on a threedimensional, multiblock, structured CFD solver for the NavierStokes equations. Structural modal dynamic equations were solved simultaneously and were strongly coupled with the flow equations using fully implicit iterativee) and semiexplicit (noniterative) timemarching methods. Since the structure deformation is usually small, a linear structure model was found to be sufficient. To address the unsteady flow around deforming structures, since the flow can be complex because of compressibility, existence of shock waves, and effects of viscosity and turbulence, a more complex model was required. The flow solver addresses the full 3D Reynoldsaveraged NavierStokes (RANS) equations with wellvalidated turbulence models. The solver also has the capability to include effects for multiblock moving boundary treatment. Robust interfacing techniques were also embedded in the coupled solver to account for transfer of information between the two modules. Our study aimed to expand a wellvalidated CFD approach to coupled aeroelastic models and consider the complexity of coupling procedures in 3D wing models. A non iterative flow solver was used for flow computations, greatly reducing the overall cost of computations (as the fluids module is the most time consuming among all the modules). In developing this model, the following issues were addressed: * efficient moving boundary technique for multiblock structured grids * preservation of geometric conservation law * choice of time step of fluid and structure solvers * accurate computation of contravariant velocities using momentum interpolation method for collocated grids. The focus of this work was to study the fluidstructure interaction problem for 3D wing geometries. The flutter boundary was predicted for a transonic Mach number case. The AGARD 445.6 wing (Yates 1987) was used to demonstrate the methodology. This configuration was chosen because extensive research has been done in the field of aeroelasticity using this model (thus experimental and numerical results are readily available). Several flow solvers, ranging from transonic small disturbance models to full threedimensional NavierStokes solver and its thin layer approximations have been coupled to the normal modes of the structure to determine the flutter boundary for the AGARD wing geometry. Particularly, we are interested in the predicting the transonic dip observed at transonic mach numbers. This dip is important in determining the minimum velocity at which flutter can occur across the flight envelope of the vehicle; hence predicting this dip is critical. Figure 12 shows the measured and computed flutter boundary for the AGARD wing using several numerical methods. The solid line represents the measured flutter boundary originally published by Yates (1963). Both linear and nonlinear aerodynamic models have been employed to determine the flutter boundary. Linear analysis using transonic small disturbance model (CAPTSD; Bennett et al. 1989) was found to predict the flutter boundaries accurately at subsonic and supersonic speeds but failed to predict the dip, accurately, at transonic Mach numbers. Specifically, linear analysis was found to be unconservative in the transonic speed regime, where it predicted a significantly higher flutter speed. This is attributed to the highly nonlinear effects arising from formation and disappearance of shock waves at this Mach number regime as the aircraft undergoes unsteady, flexible motion. Inclusion of viscous effects to the transonic small disturbance model (CAPTSDV; Robinson et al. 1991) increased the predictive capability of the model at transonic Mach number, as seen from the figure. Within nonlinear models, one can use both inviscid as well as viscous analysis to determine the flutter boundary. The flutter boundary obtained by solving the unsteady Euler aerodynamics equation of motion coupled to the normal modes of the structure (CFL3DEuler; LeeRausch and Batina 1995) is shown in Figure 12. The result from this method was found to be overconservative, predicting a significantly lower flutter speed. Viscous effects such as boundary layer thickening and/or flow separation due to shock waves were found to be important factors in determining the transonic dip accurately (Schuster et al. 2003). The inclusion of viscous effects was found to improve the prediction of transonic dip (LeeRausch and Batina 1996; Gordnier and Melville 2001; Liu et al. 2003). 0.75 CFL3D (Euler), LeeRausch et al. (1995) A CFL3D (NS), LeeRausch et al. (1996) o CAPTSD; Bennett et al. (1989) CAPTSDV, Robinson et al. (1991) Liu et al. (2003) x C + Gordnier and Melville (2001) + $ 0.5 *Experiment 0.25 0.4 0.6 0.8 1 1.2 Mach number Figure 12. Flutter speed index prediction for AGARD 445.6 wing using several methods LeeRausch and Batina (1996) coupled an unsteady thinlayer approximation of the NavierStokes equations with the normal modes of structure. A moving mesh method based on spring analogy was incorporated to account for grid movement after each time step. Gordnier and Melville (2000) coupled an unsteady compressible NavierStokes model with normal modes of structure using a BeamWarming type implicit time marching scheme with subiterations. An overset grid approach with algebraic mesh deformation method was used to account for grid movement. The geometric conservation law, which takes care of certain geometric quantities associated with mesh movement, was invoked as well. Liu et al. (2003) coupled an unsteady RANS model with normal modes of structure to predict the flutter boundary for the AGARD wing. Spring analogy along with transfinite interpolation technique was used to move the multiblock mesh. An implicit time stepping scheme using subiterations was employed to march in time. Flutter boundary obtained using these methods are summarized in Figure 12. Significant improvement was also seen, while using nonlinear viscous models, at supersonic speed regimes. Not all of the abovementioned models incorporated all the features essential in producing a robust CAE model. Some of the limitations that the previous CAE models face can be listed as follows * Fixedgrid computations (Bennett et al. 1989; Robinson et al. 1991) this is often the case while using transonic small disturbance model * Use ofinviscid flow solvers (Bennett et al. 1989; LeeRausch and Batina 1995)  failed to predict viscous effects such as boundary layer growth and/or flow separation * Implicit timemarching schemes with subiterations (Gordnier and Melville 2001; Liu et al. 2003 being an iterative scheme, it can be computationally expensive. Although such an implicit scheme is unconditionally stable, the choice of time step is limited by the frequency of oscillations of structure. * Failure to include geometric conservation law (LeeRausch and Batina 1995, 1996; Liu et al. 2003) essential while solving problems on moving mesh systems We aim at addressing all of the abovementioned issues and to develop a robust CAE model capable of predicting the flutter boundary of threedimensional wing geometries accurately. 9 A review of the various existing methods in the field of computational aeroelasticity is given in Chapter 2. The governing equations used by the different models are addressed in Chapter 3. Chapter 4 discusses the computational procedure and setup involved with a CAE model. Individual modulevalidation results along with coupled simulation results and discussions are given in Chapter 5. Conclusions and thoughts on future directions are given in Chapter 6. CHAPTER 2 LITERATURE REVIEW Next we review various aspects and modules related to the field of computational aeroelasticity (Bennett and Edwards 1998; Friedmann 1999; Huttsell et al. 2001). First, the various models associated to unsteady aerodynamics are presented. Then we review various classes of CAE: fully coupled analysis (or unified fluidstructure interaction), closely coupled aeroelastic analysis, and loosely coupled analysis. We discuss advances in the field of moving mesh methods for remeshing purposes and interfacing techniques for exchanging information between different modules used in some coupled aeroelastic models. The various formulations of Geometric Conservation Law (GCL) are also reviewed. Aerodynamic Models To understand the fluidstructure interaction problem, we need to model both the structure and the fluid efficiently. However, since our emphasis was on the fluid (rather than structure) models, we first review some physical models from the fluids perspective undergoing timedependent motion. Different classes of coupled CAE models (explained earlier) and the issues associated therein are discussed later in this chapter. Physical Models Physical models used for treating fluidstructure interaction problems can vary enormously in their complexity, based on the applications. One of the simplest models is based on piston theory (Dowell and Hall, 2001), which expresses the pressure, p, at some point x, y and some time t on the oscillating body, as a simple function of the motion at the same point and instant. It can be expressed as follows PU aw aw p= +U /M S t ax where w is a function of x, y and t and it is the instantaneous deflection of the body. The symbols p, U and M represent freestream density, velocity, and Mach number, respectively. This simple method is only useful for a limited set of flow conditions, and is usually used to verify more complex models in the appropriate limit. An improved model to the piston theory is the fullpotential flow theory, which works under the assumption that the flow is inviscid and irrotational. The potential flow model solves the nonlinear wave equation for the velocity potential, from which the velocity (and thereby the pressure) can be obtained using Bernoulli's equation. If the body profile is assumed to be thin, the nonlinear equation can be cast into a linear convectedwave equation, which has found uses for many fluidstructure interaction problems such as flutter and gust response analysis (Bisplingoff et al. 1955; Fung 1955). The linear convectedwave equation has trouble satisfying the boundary conditions (Dowell and Hall 2001) because in the boundary condition, both the velocity potential and its gradient over different portions of the fluid domain are unknown (leading to a mixedboundary problem). This is resolved by reducing the convectedwave equation (partial differential equation) to an integral equation using Green's theorem or Fourier transform. This is also referred to as the boundary element approach. Another wellknown model is based on small perturbation theory (Bisplingoff et al. 1955; Fung 1955), but it was found to fail when the flow is transonic (when shock waves may appear and disappear). Another class of models is the timelinearized or dynamically linear model, in which a steadystate nonlinear solution is used as a starting point; then a small dynamic perturbation about this steady flow is considered, and all subsequent flow variables and shock motion are assumed to vary in a linear fashion. This model leads to an order of magnitude reduction in computer resources compared to the nonlinear model, and was found to be sufficient for many problems. However, this method was found to be less useful for turbomachinery problems. This approach can be extended to determine a full dynamically nonlinear solution, which involves solving a nonlinear convectedwave equation for potential flow or Euler or NavierStokes models. Either finitedifference or finitevolume schemes in spatial variables can be used to convert the system of partial difference equations to ordinary differential equations, which forms the basis for CFD. Additional models must be developed to account for turbulence flow features, and for transition from laminar to turbulent flows. Another class of models beginning to gain interest in the field of fluidstructure interaction is reducedorder modeling (ROM) techniques, discussed next. ReducedOrder Models For the past several decades, researchers have worked in the field of CFD to develop models for complex unsteady flows. The computational cost for high dimensionality model, especially for aeroelastic problems, has limited the use of full CFD models for such applications. Recently, advances are being made to develop a novel technique for unsteady flows based on the modal character of flows, which can be termed reducedorder models. In the structural dynamics world, over the years, finite element models for structural dynamics have been reduced in size by using the normal or eigenmodes of the structure, thereby reducing the model to a few degrees of freedom from thousands of degrees of freedom (Dowell and Hall 2001). This reduces the computation time for solving such problems, while maintaining the accuracy of the physical phenomena. This method has also gained interest in the field of fluid dynamics, because such an approach gives us great benefits (saving computational costs and giving insight into the dynamics of the fluid models by considering their different modal structures). This method involves constructing a computational aerodynamic model using the dominant eigenmodes of unsteady aerodynamic flows. Combining such a reducedorder aerodynamic model with a structural modal model is an efficient way to form an aeroelastic modal model with a modest number of degrees of freedom. Extracting the dominant eigenmodes for large dimensional systems can be potentially difficult. Hence another modal approach that seeks to include more information on the flow response to enhance the accuracy of the reduced model has been developed and it is called the proper orthogonal decomposition (POD) method (Ahlman et al. 2002; Zhang et al. 2003). It is a much simpler approach than the eigenmode approach, and it uses a methodology based on nonlinear dynamics and signal processing. One disadvantage of this method is that determining the POD modes can be computationally expensive compared to determining the eigenmodes. Extensive research is being done to construct nonlinear aerodynamic ROMs and to use the eigenmode ROM approach to develop better turbulence models. However, it is still unknown whether ROM or POD approach can accurately predict all the length scales associated with the turbulence models. Review of Coupled CAE Models Before looking at the various CAE models, the generalized equations of motion (Schuster et al. 2003) are given to explain CAE methodologies [M/]{q(t)}+ [C] {q(t)}+ [K] q(t)} ={F(t)} (21) N w(x, y, z, t)} = q,(t) {, (x, y, z)} (22) 1=1 where {w(x, y, z, t)} is the structural displacement at any time instant and position and {q(t)} is the generalized displacement vector. The matrices [M], [C], [K] are the generalized mass, damping, and stiffness matrices; respectively and 0 are the normal modes of the structure, with N being the total number of modes of the structure. The term on the righthand side of Eq. (21), {F(t)}, is the generalized force vector (which is responsible for linking the unsteady aerodynamics and inertial loads with the structural dynamics). Equations (21) and (22) show that the distinct terms representing the structures, aerodynamics, and dynamics disciplines give us the flexibility in choosing different methods for a given system. For example, linear structural models can be coupled with a 3D unsteady RANS model, to develop a CAE model without actually changing the overall formulation of the equations of motion. This example of a closely coupled model is the emphasis of our study. However, fully coupled models and loosely coupled or uncoupled models have been developed. Some of these models are discussed next. Fully coupled Analysis In this method, the governing equations are reformulated by combining fluid and structural equations of motion to obtain a unified set of equations, which are then solved and integrated in time simultaneously. While using a fully coupled procedure, one must deal with fluid equations in a Eulerian reference system, and structural equations in a Lagrangian system. This leads to the matrices being orders of magnitude stiffer for structure systems as compared to fluid systems, thereby making it virtually impossible to solve the equations using a monolithic computational scheme for largescale problems. Initially, Guruswamy and Byun (1993, 1994) combined Euler flow equations with plate finiteelement structures; and later combined the NavierStokes equations with shell finiteelement structure to perform fluidstructure calculations. They used a domain decomposition method, wherein fluids and structures are solved in separate modules. On the same note, Garcia and Guruswamy (1999) computed the transonic aeroelastic response of 3D wings by coupling a nonlinearbeam finiteelement model with Navier Stokes equations. This kind of fully coupled method has limitations on grid size, and is currently limited to 2D problems as they are computationally expensive. These models and the test cases used to study them are shown in Table 21. Table 21. Description and key results of a few fullycoupled analysis methods Author (s) Description of work Major Results (year) Guruswamy, Compute aeroelasticity by direct Validity of coupling plate Byun coupling using timeintegration method elements with Euler (1993, 1994) Fluid: Euler equations equation Structure: Plate finite elements Virtual surface method Aerodynamic loads are transferred by transfers loads more bilinear interpolation and by virtual accurately than bilinear surface methods interpolation technique CFD grid (151 x 30 x 35) FEM grid (36 plate elements) Fighter type wing with M=0.854 and c=l1 deg Garcia, Model for coupled nonlinear beam FEM FEM results are accurate Guruswamy model with NS solver for static except for deflections which (1999) aeroelastic analysis of high AR wings were smaller than modal Flow solver: ARC3D fluids module of results ENSAEROWING code Nonlinear and linear beam Structural code: Nonlinear beam FEM models predicted similar Aeroelastic research wing (ARW2) @ pressure coeff results M=0.85 and a=2 Geometrical nonlinearity was found to be negligible Loosely and Closely Coupled Analysis In this class of methodologies, unlike the fully coupled analysis, the structural and fluid equations are solved using two separate solvers. This can result in two different computational grids (structured or unstructured), which are not likely to coincide at the boundary. This calls for an interfacing technique to be developed, to exchange information back and forth between the two modules. This is true for both loosely and closely coupled approaches. We now review each of these methods separately. Loosely coupled analysis The loosely coupled approach has only external interaction between the fluid and structure modules; or the information is exchanged after partial or complete convergence (Smith et al. 1996a). This approach is like a multidisciplinary computing environment (MDICE) (Seigel et al. 1998), where one effectively controls the interaction between two commercial codes for each of the modules by means of interfacing techniques. This gives us the flexibility of choosing different solvers for each of the modules but the coupling procedure loses accuracy as the modules are updated only after partial or complete convergence. A typical block diagram of MDICE is shown in Figure 21. Here, the interface methodology has been divided into two categories: function matching interface and conservative interface. Function matching interfaces provide the closest match between data on the two computational grids. Conservative interfaces aim at conserving relevant properties (such as forces and momentum) during the transfer process. Fluids Module Interface Methodology Structural Module Figure 21. Sample MDICE environment for aeroelastic simulation. (Seigel et al., 1998) Closely coupled analysis This method will be the main focus of this thesis. Here, the fluid and structure equations are solved separately using different solvers but are coupled into one single module with exchange of information taking place at the interface or the boundary via an interface module. The information exchanged here are the surface loads, which are mapped from CFD grid onto CSD grid, and displacement field, which are mapped from CSD grid onto CFD grid. The transfer of surface displacement back to the CFD module implies deformation of the CFD boundary mesh and this calls for a moving boundary technique to enable remeshing the entire CFD domain for further computations as we march in time. This can cause potential problems for multiblock grids with complex geometries and will be looked at indepth shortly. Several models have been combined for individual modules to arrive at a coupled model. From the fluids perspective, models ranging from simple potential flow models to complex 3D RANS models have been used. On the other hand, models ranging from linear beam finite elements to nonlinear solid finite elements have been used for structure module. These models are interlinked via necessary interfacing techniques, the Panel methods Parabolized NavierStokes Euler equations Asymptotic expansion Boundary Layer Full NavierStokes other Function matching * Infinite plate spline * Thin plate spline * other Conservative interfaces * Infinite plate spline * Thin plate spline * Conservative/ consistent * other Modal analysis Influence coefficient Linear FEM Nonlinear EM other complexity of which depends on what two models are used for the individual modules. A brief summary of some of the models that have been developed in the past will be shown next. Cunningham et al. (1988) developed a computational scheme to perform transonic aeroelastic analysis by coupling transonic small disturbance (TSD) potential flow equations (CAPTSD) with the natural vibrational modes of the structure. Viscous effects were later incorporated into the flow solver by including an inverse integral boundary layer model. The equations of motion were solved on a sheared cartesian grid where the lifting surfaces were modeled as thin plates. This kind of approach simplified the task of generating grids and no moving boundary algorithm was required as the surface velocity boundary condition was applied at a mean plane. This technique of using TSD formulation failed in the presence of a strong shock or when viscous effects are dominant. To overcome this, Schuster et al. (1990) came up with a model that uses a 3D flow solver coupled with a linear structure model to study the aeroelastic analysis of a fighter aircraft (ENS3DAE). Thin layer approximations to the full threedimensional compressible RANS equations were used. A threedimensional implementation of the BeamWarming implicit scheme was employed for temporal integration. The equations were solved on multiblock curvilinear grids. The linear generalized mode shapes were used to model the structure. A grid motion algorithm that uses an algebraic shearing technique was used to account for the grid movement. A similar method (CFL3DAE), developed by LeeRausch and Batina (1995, 1996), couples a linear, normal mode structural dynamics model with the thinlayer three dimensional compressible RANS model. Time marching was accomplished by means of a second order accurate backward time differencing scheme. A pseudo time subiteration method was introduced to expedite the convergence at each time step. A moving mesh algorithm based on spring analogy was used here. This model was used to predict the wing flutter boundary. An overview of the abovementioned models, namely, CAPTSD, ENS3DAE and CFL3DAE, have been given by Bennett and Edwards (1998) and Huttsell et al. (2001). The main features and results of these methods are shown in Table 22. Table 22. Description of CAE simulations using CAPTSD, ENS3DAE and CFL3DAE Author (s) Description of work Major Results (year) Cunningham, Computational scheme for transonic Aerodynamic forces and Batina, Bennett aeroelastic analysis to perform flutter flutter characteristics (1988) analysis obtained using linear *Flow: Transonic small disturbance formulation compared well formulation with expt. Structure: Lagrange Equations of motion *Nonlinear flutter results based on the natural vibrational modes compared well with expt AGARD configuration with 45 deg but not so with linear sweep angle and M=0.3381.141 results Can treat configurations with arbitrary lifting surfaces Lewis and *External aeroelastic simulation for Results showed the engine Smith internal aerodynamics and shell liner to be dynamically (1998) structures stable SFlow: ENS3D Inner flow Mach no. had Predictorcorrector scheme for structural little effect on aeroelastic integration response Tested on an engine liner to study flutter Effect of pressure loadings with M=0.7 in inner region and M=0.4 on the shell structures were in the annular region not considered in this method Table 22. Continued Author (s) Description of work Major Results (year) Schuster, A 3D flow solver coupled with linear Aeroelastic analysis Vadyak, static structural model to study compared well with Atta aeroelastic response of aircraft experiment with respect to (1990) Grid deflection method is used to update pressure coefficient and the grid after each time step. twist SFlow solver: ENS3D Flexible wing/body Swept, tapered wing with constant cross configuration gave better section with M=0.9 and a=9 deg was results compared to rigid used body configuration oWing mesh: 92 x 32 x 32 points Separation on the upper surface was not predicted Lee Rausch *NavierStokes aerodynamics to compute *Difference in flutter speed and Batina AGARD 445.6 wing flutter index and frequency index (1993, 1995, *Flow: Implicit upwind Euler/NS solver between Euler and NS 1996) Structure: Modal analysis solver was pointed out SMoving mesh: Spring analogy Grid: 193 x 41 x 65 CH type M M=0.96, Re=364,600 per foot of chord Hartwich, Study LCO for a B1 configuration using *Predicted aerodynamic Dobbs, Arslan NS equations damping matched well with and Kim Flow: CFL3D a 3D NS solver experimental trends (2000) Structure: Lagrange's equations of Fell short of predicting a motion true LCO phenomenon Moving mesh: Spring analogy and TFI using master/slave concept *Grid: 281 x 137 x 65 CO type _M=0.975, a=7.38 deg and Re=5,900,000 Liu et al. (2000, 2003) presented an integrated CFDCSD code for flutter calculations based on a parallel, multiblock, multigrid flow solver for solving the full NavierStokes equations. The flow solver is strongly coupled with the structural modal dynamics equations. A dual timestepping scheme was introduced to enable simultaneous integration of flow and structural equations without a time delay. A moving mesh method based on transfinite interpolation (TFI) (Eriksson, 1981) and spring analogy (Hartwich and Agrawal, 1997) was also incorporated in the code. Message passing interface (MPI) was used to enable data transfer between the two modules. The method was tested to perform the static aeroelastic analysis and the wing flutter on the AGARD 445.6 wing. The key results from this model are shown in Table 23. Table 23. Summary of work with a moving mesh algorithm Author (s) Description of work Major Results (year) Liu, Cai, Zhu, AGARD 445.6 Wing flutter using a Flutter speed/frequency in Wong and Tsai coupled CFDCSD good agreement with (2000) Flow: Parallel multiblock Euler experiment Structure: Modal dynamic equations Transonic dip captured *Moving mesh: Arclength based TFI and spring analogy SInterface: Transformation spline matrix *Grid: 176,601 points (32 blocks) *M=0.3381.141 Cai, Liu and Static aeroelasticity of AGARD 445.6 Convergence was spedup Tsai (2001) wing using Euler/NS equations using relaxation technique. *Flow: Parallel multiblock NS *Difference in solutions Structure: Static elastic equations between rigid and flexible Moving mesh: Spring analogy and TFI wing were spotted _M=0.85 and a=5 deg A threefield formulation for solving transient nonlinear aeroelastic problems was suggested by Farhat et al. (2000) where they used an Arbitrary Lagrangian and Eulerian (ALE) method for solving the equations on a deforming mesh. In fact, most CAE problems can be formulated as a threefield problem: the fluid, the structure and the moving mesh. In the case of ALE formulation, separate set of equations are specified for grid movement that are directly coupled with the ALE flow equations. The fluid and structure equations are coupled by the interface conditions. Unstructured meshes were used for both fluid and structure solver. Farhat and Lesoinne (2000) improved upon the existing serial and parallel algorithms for nonlinear transient aeroelastic problems. A review of some of these methodologies is presented in Table 24. Table 24. Summary of work related to ALE formulation Author (s) Description of work Major Results (year) Farhat, Pierson Computational method to simulate Qualitative validation of and Degand transient aeroelastic response of flexible results was done (2000) aircraft during highG maneuvers Geometric conservation *Flow: Arbitrary LagrangianEuler law was incorporated equations are incorporated into the Viscous effects were unstructured flow solver (Euler) neglected Structure: Corotational formulation SM=0.901 and c=1 deg on Langley fighter Farhat and Serial and Parallel methodologies for *Partitioned algorithms were Lesoinne nonlinear transient aeroelastic problems found to be efficient than (2000) Flow: ALE formulation monolithic schemes SMoving mesh: Dynamic mesh equations coupled with the flow equations *M=0.901 on an AGARD 445.6 wing Geuzine, Threefield formulation for flutter *Energy conservative Brown and analysis of F16 configuration exchange of aerodynamic Farhat (2002) *Flow: ALE formulation and elastodynamic data was Structure: Elastodynamic equations shown Moving mesh: Dynamic mesh equations *Method was found to be combined with flow eqns. effective in the transonic *M=0.71.4 on F16 wing regime and not as effective Grid size: 403,919 (63,044 on wing in the subsonic and surface) supersonic regime Review of Moving Boundary Models Having reviewed the various developments in the field of computational aeroelasticity as far as coupling procedure, our focus shift towards one of the most key aspects of computational aeroelasticty, which is the deforming mesh method. Since the structure movement needs to be accounted for in the fluid domain, we need to ensure that the entire flow domain is remeshed appropriately. Also, an efficient moving boundary module is very important for performing unsteady flow calculations such as flutter simulation of wings and turbomachinery blades. Since the grid needs to be updated frequently in unsteady computations, a fast and automatic grid deformation procedure is an essential feature. Several models have been developed over the past decade and we will review some of the methods in this section and point out the advantages and disadvantages, if any. Initially, a spring analogy method, originally proposed by Batina (1989) for unstructured grids and later expanded by Robinson et al. (1991) to structured grids, was used to generate dynamic grids for structured and unstructured solvers. This method can handle large deformations but, being an iterative method resembling an elliptic grid generator, it was found to be computational expensive for larger grid sizes. Schuster et al. (1990) and Bhardwaj et al. (1998) used a simple algebraic shearing technique to deform the grid by redistributing the grid points along grid lines that are in the direction normal to the surface. This method can cause potential problems when the geometry becomes complex when it becomes difficult to locate the radial direction normal to the surface. Also, this method is limited to small deformations and large deformations may lead to poor grid quality and crossover of grid lines. A transfinite interpolation (TFI) method (Eriksson, 1982) is typically used for regenerating individual blocks in multiblock meshes. Hartwich and Agrawal (1997) combined the spring analogy method with the TFI method for regenerating multiblock grids. Spring analogy was used to move the boundary edges of the blocks whereas TFI was used to remesh the surface and interior volume of each block. A pointbypoint match was enforced between two abutting blocks. Potsdam and Guruswamy (2001) improved the above method and incorporated parallelization for mesh regeneration. Another class of methods for remeshing purposes is solving the moving mesh partial differential equations (Huang et al., 1994; Huang and Russell, 1999; Huang, 2001). In this method, a mesh equation is formulated and solved to move the nodes in a consistent fashion by accounting for clustering of nodes in regions of large solution variation. A monitor function was incorporated into the equation to enable mesh smoothing. This method can be computationally expensive for complex 3D problems. A comparison of some of the abovementioned methods is shown in Table 25. Table 25. Comparison of moving mesh algorithms Method Advantage Disadvantage Spring analogy Robust Needs more Memory and (Robinson et al., 1991) CPU Transfinite interpolation May not preserve original (Erikkson, 1982) grid quality Gordon's TFI based method (Wong et al., 2000) Eriksson's TFI based Er s TI b d Faster and Preserves grid May encounter crossover method (Hartwich and d (Harwich ad quality near the moving boundary Agrawal, 1997) Perturbation method (Reuther et al., 1996) M g mh p l Easy to implement and Moving mesh partial accounts for grid quality differential equation accounts for grid quality Computationally expensive near regions or large (MMPDE) (Huang, 2001) grains gradients Review of Geometric Conservation Law A key aspect of solving problems on a deforming grid is to ensure that the Geometric Conservation Law (GCL) is preserved. It takes care of certain geometric quantities associated with the deformed grid or the new grid. In the numerical perspective, it is called the discrete geometric conservation law (DGCL). The DGCL states that the computation of the geometric quantities associated with a moving grid should be computed in such a way that, independent of the mesh movement, the numerical scheme used for integrating the flow equations must preserve a uniform flow field (Guillard and Farhat 2000). This is in conjunction with the fact that preserving uniform field implies first order accuracy. In addition, Guillard and Farhat (2000) showed thatfor a porder timeaccurate scheme on afixed mesh, satisfying the corresponding p order DGCL is a sufficient condition for the scheme to be at least first order time accurate on a moving mesh. They established the requirement that preserving the uniform flow field on moving grids is related to a consistency condition. It has also been proven that not satisfying the DGCL introduces a weak instability in the numerical solution on moving grids (Lesoinne and Farhat, 1996). Substantial evidence exists showing that not satisfying the geometric conservation law leads to erroneous solutions or spurious oscillations in the solution (Guillard and Farhat 2000; Lesoinne and Farhat, 1996; Farhat et al., 2001 & 2003). For example, Shyy et al. (1996) demonstrated that without explicitly enforcing GCL, 0(1) error could be induced in the computation simply due to the grid movement effect. It has also been shown that satisfying the DGCL can improve the timeaccuracy of computations on moving grids (Koobus and Farhat, 1999). One of the widely used methods for fluid structure interaction problems is the ALE formulation. It formulates the NavierStokes equations in three coordinate systems namely, material or Lagrangian (for structure motion), spatial or Eulerian (for fluid motion) and referential (for grid movement). Farhat et al. (2001, 2003) showed that for ALE schemes, satisfying the DGCL leads to a necessary and sufficient condition for the numerical scheme to preserve nonlinear stability on a fixed grid. However, there have been a few cases where satisfying or not satisfying the GCL produced the same results (Morton et al., 1998). It should be noted that since GCL arises due to the numerical procedures devised based on grid movement, its implications are expected to be scheme dependent. Alternative forms of the GCL have been implemented over the years to study its impact on solution accuracy. Thomas and Lombard (1979) implemented the GCL for density based finite difference schemes on structured meshes by updating the value of the Jacobian at each time step. Shyy et al. (1996, 2001) implemented the GCL along the lines of Thomas and Lombard for pressurebasedfinite volume schemes by updating the Jacobian values after every time step using a first order backward Euler timeintegration scheme. Lesoinne and Farhat (1996) developed a first order, time accurate scheme preserving the GCL using the densitybased ALEfinite volume as well asfinite element schemes on unstructured grids. Koobus and Farhat (1999) proposed a GCL scheme for secondorder timeaccurate densitybased ALEfinite volume schemes. Farhat et al.(2001) summarized six different timeintegration schemes based on ALE formulation, some of them preserving the DGCL and some of them that did not, and showed the impact the different schemes have on solution accuracy. In this effort, we assess selected approaches for multiblock structured grids based onfinite volume formulation and do a comparative study on these methods. Most previously conducted studies employed the densitybased fluid flow solver; in the present effort, the pressurebased fluid flow solver (Shyy, 1994; Shyy et al., 1997 and Thakur et al., 2002) is utilized. The implications of different implementation of GCL and the fluid flow solver are of main interest. Together with the previously cited references, the present work offers a more complete assessment of the GCL. Review of Interfacing Techniques Having looked at the three major modules required for aeroelastic computations, namely, fluid, structure and moving mesh modules, we now take a look at the interfacing technique that links these individual modules in an efficient manner. For coupled analysis, the exchange of information between the fluid and structure models takes place at the common boundaries. A typical coupled fluid structure analysis diagram is shown in Figure 22. The interfacing module is highlighted here for convenience. As can be seen from the figure, for every time step, we need to map the surface loads, P, from the CFD grid system onto the structural grid to obtain the forces, F, on the CSD grid system, which are then used to obtain the displacements, w, on the CSD grid. These w's need to be interpolated onto the CFD grid to obtain the CFD surface grid. Map pressure CFD p F CSD > to FEM grid  Fluid/Structure Interface Move Interpolate to Grid CFD grid Figure 22. Coupled fluidstructure flow diagram. (Guruswamy, 2002) FLUID STRUCTURE I i iIi i Ii is S11111 iii I! i / .... .. J i i \ iy1 im in ,, o Figure 23. Varying levels of complexity in modeling for fluids and structures (Guruswamy, 2002) Since the fluid and structural module can be modeled at different levels of complexity, the fidelity of the interfacing technique depends on how the fluid and structure are modeled. This has been depicted in Figure 23. Maintaining accuracy in the data exchange process is very important in order to obtain correct aeroelastic results. Often times, the structural grid is unstructured or coarser than the CFD grid, thereby demanding accurate interpolation techniques to transfer surface loads from the CFD grid on to the structural grid. We will now review a few interpolation/extrapolation techniques employed in the recent years to accomplish this data exchange. Table 26. Summary of representative interface techniques Interface method Limitations *Infinite plate spline (IPS): based on superposition of the solutions for the PDE of equilibrium for an infinite plate * Multiquadraticbiharmonic (MQ): interpolation technique that represents an irregular surface makes use if quadratic basis functions * Thin plate spline (TPS): Characterizes an irregular surface by using functions that minimize an energy functional *Finite plate spline (FPS): Uses plate bending elements to represent a planform by a number of quadrilateral or triangular elements *Nonuniform Bsplines (NUBS): uses the fact that a 3D surface can be represented by a tensor product of 2 splines *Inverse isoparametric mapping (IIM): based on FEM scheme where an isoparametric element uses shape functions to perform interpolation *Minimum of 3 grid points required *Noncoincident points are required *Extrapolations are linear *No minimum number of grid points required but 3 are preferred for accuracy *No minimum number of grid points required but 3 are preferred for accuracy * Only 2D application was looked at * Four curves and four data points required *Points cannot be coincident * Valid for 2D interpolations only *No extrapolation possible Smith et al. (1996b, 2000) reviewed six interpolation methods: Infiniteplate splines (IPS), finiteplate splines (FPS), multiquadricbiharmonics (MQ), thinplate splines (TPS), NonUniform BSplines (NUBS) and Inverse Isoparametric Mapping (IIM). Moyroud et al. (2000) demonstrated a technique based on parent volume grid and child surface grid concept to perform interpolation on threedimensional unstructured triangular grids. A brief description along with the limitations of some of these methods is given in Table 26. Guruswamy (2002) reviewed interfacing techniques based on specific finite element techniques employed for the structural model. The flow solver used was the Euler/Navierstokes solver. The FE models considered were modal model, beam finite elements, plate/shell finite elements, wingbox FE model and the detailed FE model. For the modal analysis, where the structural modes are evaluated using the RaleighRitz approach, a simple bilinear interpolation method proved to be an accurate method for structured mesh systems. For the case when the structure mesh had irregular meshes, an area coordinate approach was used. When beam structures are employed, load vectors were used along with the shape functions to output transverse displacement, twist and bending along the elastic axis for different spanwise locations. When plate or shell elements are used as the finite element structures, a nodetoelement approach was used where shape functions were used to define the coordinates and planar displacements of the element. Another method found to be effective for plate/shell elements was the virtual surface method where a mapping matrix is used to exchange information between the two grids. More details of this approach can be found in Guruswamy (2002). When the wing is modeled as a wingbox, where only the components between the spars and ribs are considered for modeling purposes, a discrepancy might occur as there is a discontinuity in surface at the leading and trailing edges. In such cases, forces are lumped onto structural nodes and bending and twisting moment conservation is enforced. Deflection at the FEM nodes were obtained by using transformation functions by assuming that the wing is chordwise rigid. Brown (1997) proposed a method that combines the nodeto element approach used for plate/shell FE and the lumped method for wingbox structures. For detailed FE models, where the interior of the FE grid could be irregular and the surface elements could take both triangular and quadrilateral elements, the area coordinate method of the virtual surface method was found to be an efficient one. A different approach called the boundary element method was proposed by Chen and Gao (2001), Chen and Jadic (2000) and Chen and Hill (1999) to perform displacement interpolations between the two grid systems. In this method, a universal spline matrix is generated to transform the structural displacement, us, to aerodynamic displacement, ua. It is given by {ua} = [B]{us}, where [B] is the spline matrix. Brief description of this method is demonstrated in Table 27. Table 27. Summary of Boundary element methods Author's Name Description of work Major Results Chen, Jadic Direct boundary element method (BEM) Performs force transferal (1998) solver for CFD/CSD interfacing with good accuracy (2D case) Generation of universal spline matrix (a Performs accurate vector) to go back and forth between displacement extrapolation Chen, Hill CFD/CSD data CSD grid points should lie (1999) *Exterior BEM solver for CFD grid inside CFD surface grid. (3D case) regeneration Boundary element near Code used: ENS3DAE leading/trailing edge causes AGARD 445.6 at M=0.95 and ca=2 instability. CFD grid (63 x 26 points) Structure grid (121 points) Table 27. Continued Author's Name Description of work Major Results Chen, Gao, Indirect boundary element method Gives very good (2001) (IBEM) solver for CFD/CSD interfacing extrapolation results on the Multiblock BEM method to handle CFD grid discontinuous structures *Eliminates edge effects AGARD 445.6 found in the direct BEM CFD grid (145 X 37 points) solver Structure grid (121 points) Deals with complex configurations CHAPTER 3 GOVERNING EQUATIONS AND OVERVIEW OF ALGORITHM In this chapter, we discuss the formulation of the governing equations of various modules associated with aeroelasticity and then look at the numerical schemes associated with these modules. We will categorize them into different categories and describe each module in detail. Governing Equations We first take a look at the governing equations associated with the various modules used in our computations starting with the flow solver. Flow Module NavierStokes equations We use a full 3D compressible NavierStokes solver as our CFD model. The equations written in cartesian coordinates, using indicial notations, read as follows Continuity: S pu )= 0 (31) at axi Momentum: (pu)+ u,)= (32) at Ex n rx gx Energy: a (pH)+ (puH)=p q + (u ) (33) at ax at ax, axi where x, is the position vector, t is time, p is density, u, is velocity vector, p is pressure, C, is viscous stress tensor, q, is heat flux vector, obtained from Fourier's law, given by aT p a h qj = K = (34) axx PrL ax, where p is the molecular viscosity, Ais the thermal conductivity, and PrL is the laminar Prandtl number defined as: PrL = H is stagnation enthalpy given by 1 H =h+u,u (35) with h being the specific enthalpy. The constitutive relation between stress and strain rate for Newtonian fluid is used to relate the components of the stress tensor to velocity gradients: ( au, au'_ 2 au,, =/ pI ~ u (36) S xax ax 3 ax, Transformation to curvilinear coordinates For arbitraryshaped geometries, it is efficient to use bodyfitted curvilinear coordinates. We denote the curvilinear coordinates as (, r, g) where =J(x,y,z,t), ri= 7(x,y,z,t) and g=(x,y,z,t). The transformation of the physical domain (x,y,z) to the computational domain (j, r, ) is achieved via transformation metrics, which are related to the physical, coordinates as follows. ll jy J\\ J\2 JZ3 1T 17y 1h =7 f1 22 f23 (37) YX y A_ l A2 A33 where the metrics f's are defined as follows f I = zy f2 = Z7xl X'z, f3 = XIy y1 X f21 = Z y y z f22 = Z y y1z f23 = z4y y z (38) A31 = y5Z z y f32 = Z 4X X7 A 33 and J is the Jacobian given by J = xJy17Z + Xyqz+ + X'ycz XJycz xY7Z Xy zc (39) The governing equations can be rewritten in generalized bodyfitted coordinates as follows: Continuity equation (Jp) +a(pU)+ (pV)+ (pW)=0 (310) at aj a7 a umomentum equation a( J p u ) + a ++ a + a T r \ ap + fp + A ap 1 + u( au )+(pu)+(p1u)= + f21 31 a ( au + u u 13F( u 2u a3 l + (311) +q q a a31 _32 a 33 aC} vmomentum equation a(pV) \p p IpW +p +(pUv)+ (pVv)+ (pW f2 1 22 +32 av + +av+qv + qV+q+ +q av (312) S q 13 ] 2 23 ] + (312) SI r{ a3v a3v avYa 31 32 3 3 33 a v wmomentum equation a(pw) at + (pTUw)+ ;(plv)+ (pww)= f3 + f3 aj a77 g aj a q I +q +qT + a [q21 a [ +w a+w v aw 7 q 3+q32 +q33 aw (vw +qw} +q22 w+q23 +jJ+ aj aq a Energy equation 3(pH) [ I h h +q r, +(p)+( +()+(pWH)= \ +,++q + at 3j a77" g aJ aj a3 aj)] a rh ah 3 h ah3Y a r h h 3 h T3 h S1+q2+23 + + +3 + 3/7 J 85 3^ 3/ 31 3 8 837 3^ Sk Sk Sk k 3 Sk Sk Sk a()< J ak 3^ 3/g 7 3^J 3' L 3 3 3 3J a k ak ak ak + q7 3 a+q2 +q3 +3 3 a J I 37 3^ 3p+ a77 (313) (314) = f (q,,u+fU) +(q, u+fV) +(q, u+fW) +  (q2lu + fziU) +(q2u+ f2V) +(q2u+ fzW) + ag aj a7 a7  (q + f1U)+(q2u+ f31V) +(q33+ f1W) + a r f av av av ( JF + f2U) +(q12v+ f2V)+(q13v+ f2W)al 6r J d aj a7r a r av av a 3v1 (q1v f2U) +(q22v + fV) +(q23v+ f22W) + a7 J a go] SF (q + f3U) V +(q32V + f32V) + (q33V+ f3W) + F (q2,1( + f3U) +(q,2w+ f3) +(q13w+ f3W) + 3(qW + f23U) + (q22W + f23V) + (q23W + f23W) + a37 J [ a37 3^J (q wVV+ f3U) + (q32w + f33V) + (33W + f33W)  where F = u+,ut F= / + r, = + t (315) PrL Prt k q11 =f +fz+f2 =12 q21= Aflf2 + 1f2f22 +f3f23 q22= f2 + +2 q13 = 31 = lf31 + 2f32 + f3f31 (316) q33 = +32 +f323 q = q32 = 31f21 + f32f22 + ff23 Here, U, Vand Ware the components of the contravariant velocity, which are the scalar products of velocity vector and area vector at a control volume interface. They can be interpreted as the volume flux normal to control volume interfaces; specifically, Uis the local volume flux along the j coordinate, V along the q coordinate and W along the g coordinate. They are given by U= f(u x) + f(vy) +f (w 2) V= f21 (u +f22 v)f23 (w ) (317) W = f31 ( + 2(v + (w ) where x, y and z are the grid velocities which are approximated by a first order backward time difference scheme given by xx yy zz x= y= z= (318) At At At Here, At is the time step size of the flow solver and the superscript o refers to the previous time step. Since the flow equations are solved on a moving mesh system, care must be taken to preserve the geometric conservation law (GCL) as originally formulated by Thomas and Lombard (1972). Geometric Conservation Law The GCL is derived from the conservation of mass by setting p=1 and v=0. It can be written as follows: d dV = ws.dS (319) dtV S where w, is local velocity of cell boundary. The GCL can also be states as the change in volume of each control volume between two time instants, tn and t+1, must be equal to the volume swept by the cell boundary during that time At= t"+t. The above expression is referred to as the integral form of GCL. A differential statement of the GCL can be derived from the integral statement of GCL. We first perform a transformation from the cartesian coordinate system (x,y,z) to the bodyfitted coordinate system (j, r, ,), which leads to the following form of the integral statement: d Jdidrld = (V.w,)Jdrdid (320) dt V Here, Jrepresents the volume element in the transformed coordinate system hence each node is associated with a particular value of J. Therefore, the computed value of J must be consistent with the value of AV implied by the numerical scheme used for solving the flow equations. Earlier, arbitrary procedures were used to compute J, for e.g., instantaneous mesh distribution at a given time instant was used to evaluate J at that particular time, which lead to an erroneous solution. Expanding the right hand side of Eq. (320) and after performing necessary manipulations, we arrive at the following form for the differential statement of GCL. J, +( ) +(7 ) +( )+ =0 (321) where, J, are the metric terms given by , =[ x yz,y z' )+.y(izx zx,)+ z(X7y x y,)] l = [(yz yz, )+ y(zx zx,)+ (xy, xy)] (322) S= (yz yz )+y (zx zx)+ (xy Xy)] Equation (321) is solved numerically to update the Jacobian values at each time step. The numerical solution for Eq. (321) requires only an initial condition, which is obtained from the initial fixed grid and is given by Eq. (39). Turbulence Modeling The conservation equations (31) (33) hold good only for laminar flows. For turbulent flows, we need to reconstruct the equations based on an averaging process in the context of RANS methods. Two kinds of averaging can be used for the averaging process: (1) Reynolds averaging, or (2) Favre averaging. In solving for compressible flows, the use of Reynolds averaging alone introduces correlations involving density fluctuations and these can be difficult to model. Hence, a combination of both Reynolds and Favre averaging is used to overcome this. This also leads to a set of governing equations that resemble the incompressible counterpart. The averaging process of different flow properties is shown below. u,=u,+u " T=T+T" h=h+h" q, = q +q " p=p+p' Mass weighted Favre averaging is used for u1, j, T, h and qj (tilde and double prime denote Favreaveraged mean and fluctuating components respectively) whereas Reynolds averaging is used for p andp (bar and prime denote Reynoldsaveraged mean and fluctuating components respectively) to recast the conservation equation. The mean flow governing equations now take the form: ( pau)= 0 (323) at ax p + (T:Jpu ) (324) at ' ax )' I ax ax 3a~ 1 I a ~ 1 ,' ) ap a + 1 H+ Y +pU upuu, at 2 axx 2 at ax( (325) aI ' + a r[ (, p +u [2 U (u'u ) As can be seen from the equation, the averaging process leads to additional unknowns in the form of Reynolds stresses, pu'u which needs to be modeled. This leads to the wellknown closure problem for turbulent flow computations. In order to determine these quantities, turbulence models are required. The two major approaches to model these terms are eddy viscosity models and Reynoldsstress models. The former approach will be used in the present study. The eddy viscosity model makes use of the Boussinesq's eddy viscosity hypothesis, which is based on the assumption that the Reynolds stresses are a local property of the mean flow and are related to the mean flow gradients via a turbulent viscosity given by: c)U, cu,)U 2 aui 2 pu,u = 4,  + 4 8, pk 8, (326) axi a ) 3 ax, 3 where /u is the turbulent viscosity that needs to be modeled. ketransport equations The transport equations for the turbulent kinetic energy, k, and eddy viscosity, C, can be written by means of one equation given by: a(p)e + (pu,) = Ui+ oli +R +R, (327) where t k or with the following expressions for and R2 where 0 = k or with the following expressions for R, and R2 tuR for the k equation R, = C ,uR (328) RR for the c equation k p C p2 k k for the k equation = at (329) C2P k C2P k for the e equation with iau f2av2 2 au' av ] u aw 12 av aw (330) R= +r\+ + + + + + + (330) [ax) ay a z ay ax /az ax _az ay where C Cs, C2, Gk and a are constants, the values of which, as suggested by Launder and Spalding (1972) are as follows C =0.09, C, =1.44, C12 =1.92, Gk=1.0 and a, =1.3 (331) The abovementioned governing equations for k and care transformed to body fitted coordinates, integrated over the control volume and discretized in a manner similar to the momentum equation, which will be described later. The turbulent viscosity, ltt, is the medium through which the time and length scale effects of turbulent flows are introduced into the equations. Thus, modeling ltt requires specification of local length and time scales. The kc models provide these scales via the modeled turbulent kinetic energy (k) and the rate of viscous dissipation of turbulent kinetic energy (c). Dimensional analysis yields the following expression for turbulent viscosity: , pCk2 (332) ju= E where C, is a dimensionless constant to be defined later. Filterbased turbulence model for unsteady RANS computations One of the fundamental problems with unsteady RANS model is that it is based on steady state mean flow data whereas the main turbulent scales returned by the model are the macro or the most energetic length scales. Hence, these models are found to have difficulties in resolving flow structures resulting from multiple length scales. The model is also known to overpredict turbulent production and hence effective viscosity in stagnation flow regions. A filterbased model (Johansen et al., 2003) is investigated here to improve upon the predictive capability of the standard ketwoequation model. In RANS computations, the true resolution is not only dictated by the mesh size but also by the magnitude of eddy viscosity. These in turn affect the local Reynolds number, which needs to be 0(1) magnitude in order to resolve the flow structure satisfactorily. Also, one should note that an excessive eddy viscosity can smear out the flow structures within the reach of a grid resolution. In such cases, the effective viscosity in the model should be reduced to resolve the structures satisfactorily. To achieve this, a filter is imposed on the turbulence model via eddy viscosity, which does not resolve structures smaller than the filter size. The filter size is chosen based on the maximum grid size in the domain. The filterbased viscosity model can be summarized as follows ,t = pC, Min 1, ki2 (333) This choice of the function, fMin 1, k assures that in near wall nodes the function will always return a value of 1.0 and hence wall functions can still be applied at near wall regions. More details about the filterbased model can be found in Johansen et al. (2003). Boundary conditions Boundary condition needs to be specified due to the elliptic nature of the kc equations. All the variables at the inlet are either specified or estimated. If the inlet velocity is specified, k is estimated based on a fraction of the square of the inlet velocity and cis estimated from k and a characteristic length scale representing the size of the turbulent eddies (usually a fraction of the inlet dimension). The estimates for k and care given by 32 k= (U )2 (334) C3/4k3/2 0.07L where T is the turbulence intensity usually set in the range of 0.02 and 0.05 and L is a typical inlet dimension length scale. The outlet boundary conditions are determined via a zero gradient condition based on the assumption that the flow is convectiondominated. Wall treatment Near wall boundaries, the local Reynolds number is of the order of one and hence viscous effects are more dominant, hence the ke model cannot be used as it is formulated based on the assumption of high Reynolds number. Two approaches have been proposed to handle nearwall effects, one being the low Reynolds number model and the other being wallfunction method. The former method requires a very fine grid resolution near the wall and hence makes the computations expensive. The latter method is based on the assumption that there exists local equilibrium between production and dissipation of turbulent kinetic energy. It was proved to be an accurate and robust approximation and it does not require a fine grid near the wall or boundary. For a fully developed turbulent flow near a noslip wall, the normalized nearwall tangential velocity, assuming a twolayer structure (viscous sublayer followed by log layer), can be written as follows Sy+ for laminar sublayer (y+ <11.63) u In(Ey) for law of the wall layer (y+ > 11.63) where K is the von Karmen constant which has a value of 0.42 and E is a factor that depends on wall roughness, streamwise pressure gradient, etc. It is assigned a value of 9.793 for a smooth wall. Here, y' is the local Reynolds number given by y= p (336) and u+ is the nondimensional velocity given by u = (337) where u, is the friction velocity at the wall given by u, = (338) The expression for the shear stress at the wall, Tw,,, expressed in terms of turbulent kinetic energy, in the log layer, is given by pC14k' /2Au T.11 = (339) ln(Ey ) Structural Dynamics Model For the current structural modeling of the wing, only linear effects will be considered (Kamakoti et al., 2002). This simplification allows for a good description of the motion of the wing, without being computationally hampered by complex nonlinear effects. Since the wing is to be modeled as a linear structure, it is possible to model the deformations as a summation of different modes of deformation without looking at the complex interaction of the modes. To this end, we choose to model the wing as a linear finite element structure that can undergo bending and torsion. It is also important to note that we choose to keep the cross sections of the wing rigid, further uncoupling the bending and torsional displacements. The linear finite element that we choose to model the wing is a beam that has mass, stiffness, and damping matrices of the actual wing. Thus, the deformations become that of a BernoulliEuler beam bending and torsion, the equations for which reads as Ed2 EI dW2 = P (340) dx dx2 where p is the distributed loading (force per unit length) acting in the same direction as the outofplane displacement (w), E is the Young's modulus of the beam, and I is the area moment of inertia of the beam's cross section. Since, the cross sections of the wing are assumed to be rigid, there is a point at which the vertical displacement of the beam is a result of only the bending of the wing. This point is where the elastic axis of the wing intersects the cross section. The generalized displacements from bending and torsion are measured from this point. This has been depicted in Figure 31. elastic axis c 4y) Figure 31. Displacements Measured with respect to the Elastic Axis To find the equations of motion, Lagrange's equation was used. The equations take the form given by dt Dq ) q d q aq where q represents the generalized displacements, vertical and torsional displacements, F represents the Rayleigh dissipation function, and Q represents the generalized forces. The kinetic energy and the potential energy of the wing are given by T and V, respectively. The generalized coordinates for the wing are functions of the position of the cross section along the span of the wing and time. Here, the generalized coordinates are referred to as w, representing the classical generalized coordinates of bending, and 6, representing the classical generalized coordinates of torsion. The equations of motion that govern the structural dynamics of the wing take the wellknown form given by Mq + Cq + Kq= R(t) (342) where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, R(t) is a vector containing the generalized forces associated with aerodynamic loads, and q, q, and q are the displacement, velocity and acceleration vectors of the finite element assembly, respectively. The above equation can be solved using modal approach by composing the solution with the eigenvectors of the vibration problem, i.e., the displacement vector can be written as q = OZ; q = OZ; q = OZ (343) where a( is the modal matrix containing the eigenvectors, orthonormalized with the mass matrix, and Z is the generalized displacement. The eigenvectors are orthogonal to both mass and stiffness matrices and if we assume Rayleigh damping, it is also orthogonal to the damping matrix. Premultiplying Eq. (342) by (T, we get Z(t)+ TCOZ(t)+n2Z(t) = TR(t) (344) where (TK = 2, DTM =1 (345) The initial conditions on Z(t) are obtained using Eq. (343) at time 0, as follows Z0 = T Mq, Z0 = TM4q0 Equation (344) can be written as n individual equations, one for each mode, as follows z(t)+ 2 ,(t)+ 2z (t) = r(t) 12 n (346) r,(t)0= R(t) with initial conditions X It=0 (DMq, j^ It=0 Here 0 is the natural frequency for the ith mode and J is the corresponding damping parameter for that mode. The solution to the above equation can be obtained for each mode using direct integration algorithms. The structural solver integrates these equations of motion in time for one time step given the time step size, the pressures on structural nodes at the initial time, and the initial geometry of the wing. It uses an implicit Newmark scheme for temporal calculations and hence is not limited by a time step size unlike other explicit methods. Moving Grid Module For fluid/structure interaction problems, we must account for grid movement along the deformed surface. Since the structure moves after every time step, we need to accommodate this movement in the CFD domain. This is usually done with some type of dynamics related mesh algorithm. In our approach, we use the perturbation method for single block grid, which is a TFIlike method as proposed by Eriksson (1982). A main feature of the perturbation method is that it regenerates the computational grid based on the original grid, hence it has the capability to preserve the quality of the original grid. For a onedimensional problem, it uses a simple formula based on spring analogy, originally proposed by Batina (1989) for tetrahedral cells and later extended by Robison et al. (1991) for hexahedral cells, given by new old old new old(347) x, ~ +x 3 x xs ) (347) where x, represents the location of interior grid point, x, the location of grid point on a boundary, and S the normalized arc length along the radial mesh line measured from the outer domain. Specifically, 1(1  x_)2 + (Y1+1 Yi)2 + (Z1+1 Z)2 1 (x1+ x)2 + (1 2) + (z z (348) To use this onedimensional perturbation method, we only need to know the displacements of the end points. The positions of the interim points are solely determined based on the displacements of these end points. Here, we consider such displacements as the source of perturbation. By design, a perturbation will spread throughout the whole domain, but exerts more effect on points near the boundary perturbed. Based on Eq. (347), a more general 3stage perturbation method is proposed for threedimensional problems treated by the singleblock approach. As far as the multiblock grid movement is concerned, we use the master/slave strategy to establish a relationship between the moving surface points (master points) and vertices located at the other blocks (slave points). The movement of the master points is based on the displacements obtained from the structure solver. The movement of the slave points follows its master point. A slave point, which has the coordinate ofxs, moves when its master point moves from xm to xm. A simple but effective formula suggested by Hartwich and Agrawal (1997), based on spring analogy, is expressed as follows: ,=x, + 0(mxm) (349) where the subscripts m and s represent master and slave, respectively, and tilde () indicates the new position. 0 is the decay function; here, we use the Gaussian distribution function, as suggested by Hartwich and Agrawal (1997), given by S= exp {,8min[FACMIN,dv/(c +dm)]} (350) where dv =(x, Xm)2 +(yv m)2 +(Z Zm)2 (351) dm = J(i m)2 +( m )2 +(Zm Zm)2 (352) and eis an arbitrary small number to eliminate division by zero. In Eq. (350), the coefficient 3 affects the stiffness. A larger 3 causes the block to behave more like a rigid body and a smaller value makes the body behave like a softball. The factor, FACMIN, in Eq. (350) plays an important role in the remeshing part when the displacement of the master nodes is small. This will be explained with an example later on. Overview of Algorithm The numerical procedure employed by the different modules will be discussed in this section. Discretized Form of Equations The governing equations presented in Eq. (310)(314) are discretized on a structured grid. The velocity components and the scalar variables (pressure, density, kinetic energy, etc.) are located at the cellcenter of the control volume whereas the contravariant components of velocity (U, V, W) are located at the center of the control volume surfaces as shown in Figure 32. V u,v,w ,p.p 0 U Figure 32. Location of variables u, v and p on a 2D nonstaggered grid for the pressure based algorithm. The discretized form of the continuity equation is written as: Jp +[pU]e +[pV] + [pW] =0 (353) The discretized form of the momentum equations can be obtained in a similar manner. Here, we present the details for the umomentum equation based on first order backward Euler method. The v and wmomentum equations follow similar derivation: At + ( +T + T3 a )+f P7 + JpuJJopu i F w Jw u + fpUul(q +qiz+q ~)i)+ fp + S ajv a7w 3w pVu(q21 +q22 +q23 )+ f2p + (354) pWu (q3 +q32 +q33 )+f p =0 J 6J a77 aC ]b The above equation is comprised of three flux terms namely, convection fluxes, pressure fluxes and diffusion fluxes. A second order difference operator is used for the flux terms. The final form of the discretized umomentum equation takes the form: aUp =EUE +au u +au +aus +auu +auBu +(fJP + f2P + f3P ). +Su (355) which can also be written as au = +rubr S (356) where 'nbr' represents the six neighboring points to the point P. The various coefficients in the Eq. (355) are given by pU F pU F aE 11 IT + q11 2 J 2 J a + as +q22 2 J 2 J (357) a = +P q33 aB +q a 2 J 2 J = a + a + a, +as +a; +a +t At and the source term is given by r, u 3 e r F au au Y r ( au au\ Se e n t S. = q12 +q1 3 q2w +q 23 + 313 2+ JJJ T5 87 b (358) +p (f i+f21 +, 1 3b) Evaluation of Contravariant velocities on Nonstaggered Grid While performing computations on nonstaggered grids, the velocity field is evaluated at the cell centers whereas the contravariant components of the velocity are evaluated at the cell faces, hence necessary interpolation technique need to be employed to evaluate the contravariant components of velocity in a consistent fashion. We use the momentum interpolation technique originally proposed by Rhie and Chow (1983) to accomplish this. We rewrite the umomentum equation (Eq. (356)) as follows: appu, = I anurubr + S* I ,p (359) = H" fpi The source term S* contains the same terms as S except the component of the pressure gradient term in the Jdirection. Using indicial notation, the above equation can be written as: u 1 1 {(l)l+l 2 (fl 2 (360) A similar equation can be written for the east neighbor //H+1 7 1+1 a1 U+1 = (f11P)1+3 2 11+1 2 (361) ap +1 aP+1 We use momentum interpolation method to evaluate interface velocity based on the nodal pressure at the interface. The expression for interface velocity is given by 91=+1/2 (P+1 P) (362) aP )+1 2 P 1+1 2 The first time on the right is evaluated using Eq. (360) and (361) as follows S =I H + (1I (363) ap l/2 a a zl = [I~u_+(1 1)uH+,] + {((fp)1+1(;2p)_(/2 ) (1){(f;+p) +/2 (_fp)1+1/2} aP+1 aP+1 Here, I, is the interpolation coefficient along the xcoordinate direction. Substituting the above expression into Eq. (362), we get 1 I (1I,) +1/2 = 1/2 1, 2+1 (PI+1 P,) (T+1/2 1/2) (p+3/2 P+1/2) (364) aP )+1/2 aP a+ where the over bars represent interpolation using the nodal values on either side of the interface. Similar expressions can be written for v,+1 2 and w,+1/2. Finally, the expression for contravariant velocities at the cell face for a stationary grid (grid velocities are zero) can be evaluated using Eq. (317) as follows U,1/2 = [fil,+1/2 + f2+1v/2 + f13w1/2] (fl2 2 f2 + 2 2)(1 12 (365) 1 (12 13 +1/2 (PIP') (P+3/2 /2)I (I+/2PA1/2) a, +1/2 a+ aP, For unsteady computations, situations may arise where we might be forced to use a small time step size based on the CFL type stability condition. It has been shown previously (Choi, 1999; Yu et al., 2002a and Yu et al., 2002b) that a small time step leads to minor oscillations in pressure and velocity field while using the original RhieChow momentum interpolation method. They proposed modified momentum interpolation schemes to calculate cellface velocity to eliminate the effect of time step size on the solution. One such proposed scheme to calculate the cartesian uvelocity component at the cell face is shown here. Similar expressions follow for v and wcomponents. 1 I, (1 /I( ) u1+1/2 = U+1/2 + 1 1 (P,+1 ) (P,3/2 12) ( +1/2 1/2) l P ) +1/2 a aP + /1 /u, a+1/ 2 ul (l) 1 P ,+1/2 "+1/2 0 i (366) where a' = J'. The superscript 1 denotes previous time step. In the above equation, the At terms inside the first box bracket constitute the original RhieChow momentum interpolation scheme and the term inside the second bracket is the additional term that was proposed to eliminate the time dependency effect. The above formulation was developed for cartesian grids. We extend the above scheme to make it suitable for curvilinear grids, which means calculating contravariant velocities in an appropriate manner from the cartesian velocity components at the cell faces. The expression for u component of contravariant velocity reads as ,+1/2 = [fil +1/2 + 12V+1/2 + f13W+l/2 (f211 +/2 +(1 I,(_^P, (367) (1 f +13 (+1 P, p P (1/ (P,+3/2 P +1/2) ( 13+l/2 (1/267) aP +1/2 a1+ a, + 1 U'+12 1 (1I) (fllu;1 l+fl3wl) (f 1 + fl2V + f13W J) a\ P ,+1/2 P,+ P, Similar expressions can be derived for V and W components of contravariant velocities. The terms in the last brackets are the additional terms stemming from the unsteady terms added to the original RhieChow momentum interpolation method in Eq. (366). We will carry out test cases on a uniform cartesian grid (cavity flow) as well as a curvilinear nonuniform grid (cylinder flow) to study the impact of these terms for different time step sizes and the results will be presented at a later stage. PressureBased Flow Solver (SIMPLE) We use a pressurebased flow solver for computing both compressible and incompressible flows as well as laminar and turbulent flows. The flow solver is based on the SIMPLE (SemiImplicit Method for PressureLinked Equations) algorithm originally developed by Patankar and Spalding (1972). This method was originally developed for incompressible flows and it has been extended to solve compressible flows by modifying the pressure correction equation to include the effects of density on pressure (Shyy & Braaten 1988). Additionally, the algorithm has been extended to bodyfitted curvilinear coordinates in order to handle arbitrarilyshaped flow boundaries. This approach can handle flows at all speeds without any fundamentally different treatment for any particular flow regime. The SIMPLE algorithm can be divided into two steps, the prediction and correction steps. In the prediction step, we start with a guessed or intermediate solution of the pressure field, p* and the velocity fields, u*, v* andw* and the velocity and pressure fields are obtained by adding a correction term as follows: u=u +u v=v +v (368) w=w +w p=p +p where the primed quantities denote the correction terms. To obtain the corrected quantities, we start with the umomentum equation for the intermediate solution i.e., by substituting the predicted quantities into the momentum equations. The expression for u is given here. Expressions for v* and w* follow similar procedures. au* =au+ +a"Wu +a"Nu +a"u +a"u, +au +(f p + fP +3p ) +S, (369) Subtracting Eq. (369) from Eq. (355) and neglecting all the terms on the right hand side related to velocity (this is the main approximation of the SIMPLE algorithm), we get the following expressions for the velocity corrections. u = fLp ++f2p +f 3P' S= flP21 + 22P + f23P~ (370) w' = f31p + f32P + f33 The corresponding corrections for the contravariant velocities can be written as U' = f1' + f1' + f3w V = f21 +f22 +f23 (371) ' = f3' + f3' + f33w Substituting Eq. (370) into Eq. (371) and after necessary approximations, we obtain the following expressions for U, Vand W U =U*+ + 2 + p a" a" a" V=V*+ 2L+ff22+f23 p (372) a" a; aw "p W= W* + + A2+A3 p' au a; a< where U = f,(u X)+ f(v* y)+ (W* z) V = f2 ( )+f(v )+f23(w ) (373) w* = f3 ( X)+f32( )+f33( ) We now obtain the pressure correction equation by substituting the above expressions for the contravariant velocities into the continuity equation. The pressure correction equation can be written as app = aEPE + aww + aNPN + asps + T + + p +Sp (374) where the coefficients are given by [f2 /2 /2 /2 /2 /2 au av a a a aw aN f2 + f22 23 f22 +f23 (375) f =E + fa +f f a + af+ + af o _Jp (376) S' = p(U U +U* Un +Ub U )+ P At In this method, pressure correction is used to update contravariant velocities and the cartesian components of the velocity are used to update cartesian velocity and pressure fields. This procedure is iterated until the convergence criterion is met. An overview of the SIMPLE algorithm is shown in Figure 33. Solve discretized momentum equations u*, v* Solve pressure Iterate correction equation Correct pressure and velocity u, v, p No o Converged? Yes STOP Figure 33. Overview of the SIMPLE algorithm PISO Algorithm for unsteady computations For unsteady flow computations we use a predictorcorrector type method to avoid the iterative procedure associated with the SIMPLE family of methods. One such method is the PISO algorithm, originally developed by Issa (1985). It is modified to suit multi block, curvilinear, structured grids. First, we will look at the discretized equation for fully implicit scheme based on SIMPLE algorithm. ap + = H" (u n+l pu J p+w) f2l(P pA)t f3(Pt b,)n (377) At At In the above equation, ap is the nodal coefficient containing the special discretization terms and H" contains all the convective and diffusive flux discretization terms as well as any source term that may be present excluding the pressure gradient terms. The superscript n denotes the old time step and n+ the new time step. The implicit solution for above equation requires several iterations at each time step, which can be computationally expensive. In order to eliminate such an iterative procedure, we use the PISO algorithm, which uses a series of predictor and corrector stages to arrive at the converged solution for a given time step. The intermediate solution at the predictor stage is denoted by and the solution at the corrector stages are denoted by **, ***, etc. The solution at the last corrector stage is considered the final solution at new time step, n+1. Predictor step: In the predictor step, the u, v, wmomentum equations are solved implicitly but the pressure field is provided by the previous time step. The umomentum equation is shown below. The v and wmomentum equations follow a similar pattern a + u = H'(u)+ f( p" w)n fl(PnP ) f31(t Pb)n (378) SAt At Corrector step: In the first corrector stage, a new velocity (u**, v**, w**) is sought along with a corresponding pressure field (p*). The velocity field obtained from the predictor step (u*, v*, w *) is used to arrive at the desired quantities. The u** momentum equation is shown here. a, + u H =H(u + f(PePw)* f2i(PP)* f3l(P Pb)* (379) At At To determine* in the above equation, we combine the above equation with the continuity equation to formulate an equation for pressure correction, p', which is defined by P =P*P" Once we solve forp', we correct the velocities u*, v*, w* to obtain u**, v**, w**. Subsequent corrector steps are carried out in a similar fashion by assuming u**, v**, w**, p* as the intermediate velocity field to obtain u***, v***, w***, p** and so on. A total of four corrector stages were used to arrive at the final solution at the new time level, n+1. The PISO algorithm, being a semiimplicit method, is limited by the choice of time step based of CFL number, which is defined as follows CFL =uAt Ax where u is the free stream velocity, Ax is the grid spacing. Time step size is chosen based on free stream velocity and grid spacing so as to obtain a CFL number with an 0(1) magnitude. An overview of the PISO algorithm is shown in Figure 34. A detailed procedure of the algorithm can be found in Thakur et al. (2002). Updating Jacobian values for moving boundary treatment While formulating the abovementioned solver for the moving boundary problem, we need to make sure that the geometric conservation law is satisfied. This is ensured by updating the Jacobians in a consistent fashion. The equation that needs to be satisfied is rewritten here for convenience J,+(~ +(r), +({7)1 =0 (321) We will now look at how to calculate the Jacobian values in a consistent fashion. In this regard, we will implement four different timeintegration schemes for evaluating the Jacobian values: first order implicit scheme, first order timeaveraged scheme, second order implicit scheme and a second order timeaveraged scheme. The formulation for the different GCL schemes is as follows. START Solve discretized momentum equations with pressure field from previous time step u*, v* Solve pressure correction equation Sp Correct pressure and velocities Predictor step I Corrector step 1 u**, v**, p* Solve second pressure correction equation vr P' Corrector step 2 Correct pressure and velocities : STOP Figure 34. Overview of PISO algorithm for unsteady flow computations First order Implicit Scheme: This scheme was initially suggested by Thomas and Lombard (1972) for density based finite difference schemes and later was implemented using a pressurebased finite volume method by Shyy et al. (1996 & 2001). The difference equation for the fully implicit scheme is given by: J" =" At [(~)"+1 +(i):1 +( 1 (380) where the metric coefficients are calculated from Eq. (322) at the (n+l)th time step and the grid velocities are calculated as follows Sn+1 n n+1 n z n+1 n S x Y n+l (381) At At At Firstorder timeaveraged scheme: It was suggested (Farhat et al., 2001) that a timeaveraged evaluation of the metric terms leads to a more consistent evaluation of the Jacobian especially when the time step is large. Hence, we reformulate the governing equation (Eq. (321)) by timeaveraging the evaluation of metrics over more than one mesh configuration as opposed to evaluation of metrics on a single mesh configuration (Eq. (380)). It can be written as follows: jn+l jn (n+ )l nl (382) 111 (t) 2 +() 2 + t 2 =0 (382) At L< 1 The metric terms are computed as follows: / h\n+ 2 1 m n+1 \n ()n+ = n+ +(j)n (, )n+1 1 [(,7 )n+1 + (7 )n ( n+n ) +2 + t. t 2 2 As can be seen from the above expression that the metrics are evaluated at a mesh configuration, which is in between the mesh configurations at tn and t+1. The grid velocities for the metric term are evaluated as follows: 1 Xn+l n 1n yn+l yn n+1 n+1 n n+ x n+ y n+ z z x2 2 2 2 =(383) At At At Second order implicit scheme The difference equation for a second order implicit scheme (3point backward), as suggested by Koobus and Farhat (1999) is used here. For structured meshes, it can be written as 3 Jn+l 2Jn + IJI 2 A 2 + [( +(71 +( +1 = I 0 (384) Here, the metrics are evaluated in a similar fashion to the one used for first order implicit scheme at the (n+1)th grid configuration. The grid velocity is also calculated in a similar fashion as given by Eq. (381). Second order timeaveraged evaluation of Jacobian Here we use a timeaveraged procedure to evaluate the fluxes. Koobus and Farhat (1999) employed such a procedure for unstructured meshes and we extend such a scheme to a structured mesh here. The scheme can be written as follows. 3 j~n+ 2J + I jnI 2 2 1 jl, + ,7 Jn+ + + jn At 2. t 7 ; (385) +1 [() +() +n ( + = o n 2 77 where Sn+1 i n+l + (1 lnn n 22 2 ,2 2)+() ) (386) t) n+ 1 2 )n+l "+ nn n + ) ( )_11 Here, the averaging is done in between configurations [tn, t n+] and [tn, t"] and the mesh velocities for each of these metric terms are evaluated as follows: 1 n+l n 1 n+l n 1 n+l n n+ nx x y y n+ z z x 2 ,y = At At At 21 _" 1 yn 1 Yzz (387) xA A ,' At At At As can be seen from the above equations that we need three different mesh configurations to arrive at this scheme. All of the abovementioned schemes will be implemented using several test cases to arrive at the appropriate scheme for aeroelastic computations. Newmark Integration Method for Structure Solver In this section, the numerical procedure used to solve the structure equations of motion is presented. The equations of motion of the structure, written for the individual modes, can be written in a generalized form for any particular mode as follows S+ 2(co + o2z = r (388) Equation (388), written at a time t+At, reads as follows ;+At + 2(coz+At + o2t+At rt+At (389) In the above expression, c, o, and r,+A, are known before hand. The task is to evaluate the displacement, velocity and acceleration at the time t+At. The following expressions for velocity and displacement are formulated at time t+At first as a function of acceleration at t+At and displacement, velocity and acceleration from previous time level t. t+t = + [(P1g)zt + 'tt]At (390) zt+At =z,+ ,At+ [( /a)Y,+at, ] At2 (391) where, aand 5are parameters that are chosen based on desired stability and accuracy. For the Newmark scheme to be unconditionally stable, values of 0.25 and 0.5 are chosen for aand 6j respectively. Substituting Eqs. (390) and (391) into Eq. (389) and solving for 2,+, and then using Eqs. (390) and (391) to calculate ztA, and z,,,+, the following relation can be established for displacement, velocity and acceleration. t+At t t+At =A z, + L r,+A (392) Zt+At Zt where 1 1 ( a) 2(1 8)A (3 2K) A( At At2 A= At[1 ( a)8 2(1 6)6] 1 8 26K (8) (393) At At2[ a(Ia)af 2(1 )aAK] At(1aI 2aAc) (1a a) 02At2 L= # (394) (0) 1 2; #2 + +a ; = A7= mw At cAt m At CHAPTER 4 COMPUTATIONAL PROCEDURE AND CODE VALIDATION In this chapter, the computational procedure for performing fluidstructure interaction computations will be discussed in depth. The computational setup for the test case that is used to carry out computations will be shown as well. We also present validation results associated with several modules of the computational aeroelasticity model. Computational Procedure The overall procedure for carrying out computational aeroelastic computations can be divided into the following major steps. 1. Constructing the geometry for aeroelastic computations and also to supply appropriate boundary conditions and initial conditions. 2. Perform steady state CFD computation to obtain initial guess for starting coupled computations. 3. Perform unsteady CFD computations using steady state result as initial guess and obtain necessary aerodynamic forces on the surface of the wing. 4. Interpolate aerodynamic forces onto the structural mesh. 5. Perform CSD computation to obtain the deformation of the geometry 6. Extrapolate the displacements onto CFD surface grid. 7. Remesh CFD grid based on the deformation obtained from the CSD calculations using the moving boundary module. 8. Repeat steps 37 using current solution as the initial guess for the subsequent steps. Now we will take a closer look at the abovementioned steps along with the grid generation details, interpolation/extrapolation schemes used. Geometry definition and Computational Grids Geometry Definition The geometry used for all of our test cases for aeroelastic computations is the well documented AGARD (Advisory Group for Aerospace Research and Development) 445.6 wing (Yates, 1987). This is the first AGARD standard aeroelastic configuration. It was first tested in the Transonic Dynamics Tunnel (TDT) at the NASA Langley Research Center (Yates et al. 1963). The AGARD 445.6 wing is a swept back wing with a quarter chord sweep angle of 450 with a NACA 65A004 airfoil (4% thickness) crosssection. It has a panel aspect ratio of 1.65 and a taper ratio of 0.66. The root chord of this model is 1.833 feet with the semispan of 2.5 feet. The wing tested at NASA Langley was a semi span, wallmounted model made with laminated mahogany. A schematic of the AGARD wing is shown in Figure 41. The computational grids used by the fluid and structure solvers are described next. Figure 41. Schematic of the AGARD 445.6 wing used in the wind tunnel (Yates, 1963) Computational Grids CFD grid A CFD mesh is generated around the AGARD 445.6 wing by placing the wing in the middle of the computational domain, which has dimensions of 18x9x9 units. The geometry could be generated by using the CAD module of any commercial mesh generating software such as ICEMCFD or PATRAN, the latter being used for this case. ICEMCFD was used to construct the CFD mesh around the wing. The entire computational domain with all the blocks is shown in Figure 42. As can be seen from the figure, it is a multiblock domain comprising of 10 blocks. An Ogrid was employed around the wing to preserve orthogonality of grid near the wing. Since it is a very thin wing, care must be taken while generating mesh around the wing tip and trailing edge to avoid any negative cell volumes. The initial coarse gird had a total of 4838 points distributed over the wing surface (118 points in the chordwise direction and 41 points in the spanwise direction). The entire CFD domain had a total of 322,622 points. The CFD surface grid along with the meshing system at the leading and trailing edges are shown in Figure 43 for clarity. Figure 42. Overview of the Multiblock CFD grid Figure 43. CFD surface grid along with grid distributions at the leading and trailing edges CSD grid For the structure solver, since we use beam elements, a tenelement beam mesh spanning the semispan of the wing was used. However, in order to make the interpolation and extrapolation procedures between the CFD and structure mesh efficient k, and easier, an intermediate surface mesh was generated with QUAD4type elements. Equal width spanwise elements, four per beam element, along the spanwise direction were used to generate this intermediate mesh. The QUAD4 elements, however, were non uniform along the chordwise direction to comply with the geometry. This intermediate surface mesh on the AGARD wing was generated using PATRAN. It is comprised of 2400 elements surrounding the wing (60 along the chordwise direction and 40 equal width elements along the spanwise direction). The intermediate or temporary mesh is shown in Figure 44. Figure 44. Schematic of the FEM grid on the AGARD wing Coupling and Interfacing Procedure Grid generation is only the preliminary step towards performing coupled aeroelastic analysis. The next step is to formulate an efficient and robust interface technique that links the two different grid systems (CFD and CSD) mentioned in the previous section. This is essential while formulating a closelycoupled CAE model. The interface technique as well as the coupling procedure will be discussed at length in this section. Once the aerodynamic pressures are obtained on the CFD surface grid of the wing, they need to be interpolated onto the CSD grid to act as the external force for the structural equations of motion. A bilinear interpolation procedure was used to transfer the scalar pressures onto the structure grid. To do this, the top and bottom surface of the wing are treated as twodimensional surfaces. The pressures are then mapped from the CFD surface grid onto the corresponding intermediate CSD surface grid. This is done by locating the four CFD grid points engulfing a given CSD grid point and employing a straightforward bilinear interpolation procedure. This is demonstrated in Figure 45. Once we locate the CFD grid points engulfing the FE grid point (p, q), we calculate intermediate pressures at points (1) and (2) by a linear interpolation procedure and these intermediate pressures are then used to evaluate pressure at point (p, q) by performing another linear interpolation. Such a scheme gives an order of accuracy between one and two. (ij 1) 0+ (i++) /' ( ) CFD grid lines ( pq)  FEM grid lines (2) (i, j) : CFD grid locations (ij),' (i lj) (p,q) : FEM grid location Figure 45. Schematic to demonstrate interpolation technique Now that the surface pressure distribution on the CSD grid is obtained, they have to be converted into pressure forces. This is accomplished by computing the unit normal and the surface area for each QUAD4 element. The scalar pressures are then multiplied to the area of the element and the unit normal to obtain the pressure forces for each element. As mentioned before, the intermediate CSD grid is divided into 40 equalwidth elements along the spanwise direction. A set of 4 such elements along with all the corresponding chordwise elements are treated as a "super element", which constitutes a single beam element. Each super element has 4 equal width surface elements in the spanwise direction and 60 surface elements in the chordwise direction as shown in Figure 46. The net pressure force for each super element is then calculated by summing up all the pressure forces for the elements contributing to the super element. This is evaluated for each of the super elements and the resulting net force calculated acts as the external force for the beam element. The mass and stiffness matrices for the structure are obtained from full NASTRAN models. Figure 46. Schematic of a super element: Portion of the entire structure The structural equations of motion are now solved using the interpolated forces and the matrices to obtain the displacements. The displacements are output in terms of a vertical deflection, w and torsional angle, 0. These are used to obtain the new shape of the wing and hence the new location of the CFD surface grid location. Since BernoulliEuler beam theory is assumed, the crosssection of the wing remains unchanged. The new location of the CFD surface grid points are determined by assuming a rigid link connecting each CFD grid point to the beam element. The links are assumed to be perpendicular to the elastic axis as shown in Figure 47. Finite element nodes Rigid links  .t *t! CFD element nodes Figure 47. Sample CFD mesh superimposed on the discretized beam structure The state of the beam at any instant along the spanwise direction is given by w={w, w2 w301 62 03} , where w represents the vertical deflection and Othe twist at each spanwise section and subscripts 1, 2, 3 denote displacements in the x, y and z directions respectively. In the present study, there is no deflection in the x andy directions and no rotation about the x and z axes. In other words, the deflection of a CFD grid point P can be written as Wp = WpT + Wp' , where wp is the translation component and wR is the rotational component. The translation and rotation component corresponding to CFD spanwise grid locations are obtained by performing a linear interpolation using CSD spanwise grid points. Once the translation and rotation component for each CFD spanwise grid location is obtained, an extrapolation technique is used to obtained the new locations of the CFD surface grid. This is demonstrated through Figure 48. In the schematic, quantities with subscript o denote the original location of the CFD surface grid points and subscript 1 denotes the new location. By knowing w and Oat each CFD grid spanwise location, the new location is obtained using equations (395) and (396). X =XoCos0+ZoSin0 (395) z1 ZoCosOXoSin0+w (396) S (xI',zI) r G(y,t) (xo,zo) w(y,t) Elastic Axis Figure 48. Schematic to demonstrate the extrapolation procedure Once the new CFD surface grid is obtained, the CFD domain needs to be re meshed. This is accomplished using the moving boundary module explained earlier. The new CFD surface grid acts as the source of perturbation needed to enable remeshing as explained in chapter 3. The next step is to ensure that the grid velocities are updated in a consistent fashion on the boundary and also, care should be taken to preserve the geometric conservation law mentioned in the previous chapter. This entire procedure is then repeated for subsequent time steps to arrive at desired aeroelastic results. Code Validation One of the most important aspects of developing a computational tool is to validate the model with theory or experiment or prior computational results. Since there are numerous issues associated with a CAE model, this section aims at looking each of the modules individually and to ensure consistency with previously published results for these modules. First, steady state CFD results will be shown, which will be followed by various issues that one may experience with unsteady flow computations, for e.g., accurate computation of contravariant velocities (for small time step sizes) and Jacobians (for moving mesh problems). Finally, results from the simplified structure model are then compared with experimental results as well. Steadystate CFD Computations First, we validate our CFD code using the AGARD 445.6 wing defined in the previous chapter and present some steady state results. The top view of the computational domain along with the boundary conditions is shown in Figure 49. We perform a viscous compressible flow calculation based on a Reynolds number of 366,400 (based on a unit chord length at the root) and a transonic Mach number of 0.96. The steady state pressure contours on the top surface of the wing is shown in Figure 410. Pressure coefficients at different spanwise locations at the top surface of the wing are shown in Figure 411. Results from both figures are in excellent agreement with prior numerical calculations performed by Lee Rausch and Batina (1993). NO SLIP WALL O U N T SL L E E T T NO SI IP \V AI I Figure 49. Top view of the CFD domain showing the type of boundary conditions specified at different surfaces /7~7 05 Figure 410. Steady state rl=0.25 X surface pressure contours on the AGARD wing il=0.50 xlc 02 04 06 08 1 I I I I I 04 04 il=0.75 xlc 02 04 06 08 rl=0.96 01 02 04 06 08 1 I I I I I xlc 02 04 06 08 1 , , 1 1 1 1 1 1 1, ,,, ,, . Figure 411. Steady state pressure coefficient distribution at different spanwise locations on the top surface 0 1  Unsteady Computations using PISO Algorithm To validate the pressurebased algorithm for timedependent computations, the unsteady problem of flow around a square cylinder is studied. Computations are performed at a Reynolds number of 215, for which experimental results are readily available. The flow is laminar and exhibits the wellknown vortex shedding phenomenon at this Reynolds number. A perturbation is given to the initial flow field to quickly reach the periodic oscillatory flow pattern. The computational domain and the boundary conditions employed for performing the unsteady computations are shown in Figure 412. NOSLIP BOUNDARY Figure 412. Computational domain for flow past square cylinder along with imposed boundary conditions On the left boundary, the steadystate inlet condition with constant horizontal velocity is used. On the right boundary, the outlet condition employing zero gradient extrapolation of the velocity and specified pressure is invoked. The top and bottom boundaries are treated as noslip boundaries, which sets all the velocity components to zero. The plot of the periodic behavior of the crossstream component of the velocity (v) obtained by probing at a single point downstream of the cylinder is shown in Figure 413. The vcomponent of the velocity alternates due to the vortex shedding behind the cylinder. Three different time step sizes are used to perform unsteady computations using the PISO algorithm. The time step sizes chosen were 0.002, 0.001 and 0.0005 which correspond to CFL numbers of 0.27, 0.13 and 0.07 respectively, which are all well within 0(1) magnitude ensuring stability and accuracy of PISO algorithm As can be seen from the figure, minute differences are seen in the amplitude of oscillations for different time step sizes. Specifically, the amplitude of the oscillations was found to increase with decrease in time step size. The wellknown Karmen vortex street is plotted for two choices of time step sizes to visualize the differences, if any, in the vortex structure. This is shown in Figure 414. It can be seen from the plot that the frequency of vortex shedding is identical for the two choices of time steps; however, the vortex structure itself looks a little different in terms of the rotation angle of some of the structures. For verification purposes, identical computations were run using the SIMPLE family of methods using the first order backward Euler time marching scheme. The plot of v velocity is shown in Figure 415. It can be seen from the figure that similar amplitude increase pattern is seen using SIMPLE method as well. Such a difference in amplitude and vortex structure associated with different time step sizes can be attributed to the first order time marching scheme used by the PISO and SIMPLE family of algorithms. 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 03 04 0.5 0.6 0.7 0.8 0.9 1 Figure 413. Periodic oscillation of the crossstream (v) component of velocity using PISO algorithm for square cylinder case at Re=215. A 02  015 01 005 > 0 005 0 0 01 02 03 04 X 05 06 07 08 . B 02 r 01 oL 01 0 1  0 01 02 03 04 05 06 07 08 Figure 414. Vordex shedding past a square cylinder using PISO algorithm for Re=215. A)At=0.001, B) At=0.0005 02 I I I I I I I I ~1 I) I * 0.5 . dt=0.001 Sdt=0.002 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.4  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 415. Periodic oscillation of the crossstream (v) component of velocity using SIMPLE algorithm for square cylinder case at Re=215. Effect of number of stages on accuracy and stability of PISO algorithm The number of momentum and energy corrector stages can have a significant impact on the stability of the PISO algorithm. The original PISO algorithm, proposed by Issa (1985), included two momentum corrector steps to go with one energy corrector step. This method was found to be unstable for skewed meshes due to the omission of crossderivative p' terms, which is necessary when solving problems on nonorthogonal grids (Thakur and Wright 2004). A method proposed by Thakur and Wright (2004), which includes tow additional corrector stages, is used for performing unsteady computations here. As a first step to this, we study the impact of the number of corrector stages on stability and accuracy using two test cases. First, we perform unsteady computations over a square cylinder, identical to the one mentioned in the previous section. All computations were run using a time step size of 0.002, corresponding to a CFL number of 0.27. Number of stages of 2, 4 and 8 were looked at to study the accuracy and stability. The residual (pressure) history for the different number of stages is shown in Figure 416. It can be seen from the figure that using 4 or 8 corrector stages produced much better convergence, compared to using 2 corrector stages. Plot of vvelocity at a point downstream of the cylinder, using different number of stages, is shown in Figure 417. It can be inferred from the plot that the number of stages for this case did not have a significant impact on solution accuracy. Since a Cartesian type gird is used for these computations, the number of stages did not have a significant impact on the stability or the solution accuracy of the PISO algorithm. 2 stages 0.5 4 stages 8 stages 1 1.5 2.5 0 3  3.  4 4.5 5. i ^"  0 20 40 60 80 100 120 140 160 180 200 time steps Figure 416. Pressure residual history for unsteady flow over a square cylinder (Re=215) 0 1 2 3 4 5 6 7 8 9 10 01 3456789r Figure 417. Periodic oscillation of Crossstream velocity (v) using different number of stages for PISO algorithm As a second test case, unsteady, laminar flow computations, over a circular cylinder is carried out. Computations were run at a Reynolds number of 100 with a time step size of 0.002 seconds, corresponding to a CFL number of 0.27. This case is chosen to test the impact of number of stages on a curvilinear grid. The grid, along with the boundary condition, used for these computations is shown in Figure 418. Using 2 corrector stages for this case was found to be unstable due to the curvature in the geometry. The residual (pressure) plot using 4 and 8 corrector stages are shown in Minor differences were observed in the convergence tolerance obtained using the different number of stages; however, the convergence rate is in comparable range for both cases. Similar vvelocity plot at a point downstream of the cylinder is shown in Figure 420. Both cases produced identical results and the Strouhal number evaluated using either 4 or 8 stages were found to be 0.164, which is in excellent agreement with experiment (Williamson 1989). Similar results were observed for compressible, turbulent flow over a threedimensional wing as 83 well. Based on the above results, it was found that using 4 stages was necessary for performing computations on curvilinear mesh systems, as well as skewed geometries. Using 8 corrector stages increased the computational cost by twofold, compared to using 4 stages, and did not affect the accuracy significantly. Hence, we chose 4 corrector stages for all further computations. Some cases may require additional stages depending on the complexity of the problem. OUTLET I OUTLET circular cylinder r 4 2tawg F.g5 8 stage c 1.S 2 ^ 2.5 z D 4 5 40 B 100 13fl140 1M 1 0 2)0 Figure 419. Pressure residual history for unsteady flow over a circular cylinder (Re=100) 21 S3 Figure 419. Pressure residual history for unsteady flow over a circular cylinder (Re=100) S25 4 stages 0.25 J \. ^ 8 stages 0.2 f 0.15 I 0.1 0.05 0.05 0.1 0.16 0.2 o 0.25 " 0 1 2 3 4 5 6 7 8 9 10 t/T Figure 420. Periodic oscillation of crossstream velocity (v) for different number of corrector stages Momentum Interpolation Techniques for Computing Contravariant Velocities The impact of different momentum interpolation methods to compute contravariant velocities, specifically the effect of time step size, on solution accuracy will be demonstrated here. Prior work in this field showed indepth analysis for various Reynolds numbers and grid resolutions using the 2D cavity flow case (Choi, 1999). We chose a similar case to validate the proposed momentum interpolation scheme and then extend this method to study the flow around a cylinder employing curvilinear mesh. Schematic of the cavity is shown in Figure 421 along with the boundary conditions. The cavity flow case was run at a Reynolds number of 100 based on cavity height of 1 unit on a 21x21 uniform grid. Three different time step sizes of 1, 0.1 and 0.01, corresponding to CFL numbers of 40, 4 and 0.4, were chosen to study the impact of time step size on solution accuracy. The plot of uvelocity and pressure using both the original RhieChow momentum interpolation scheme as well as the modified momentum interpolation scheme is shown in Figure 422 for different time step sizes along with steady state solution for comparison. The velocity plots using both the original and modified schemes produced almost identical results, whereas noticeable deviation from steady state solution can be seen in the pressure plot for the original momentum interpolation scheme for 0.01 time step size. This is in agreement (Choi, 1999) with the finding that smaller time step size leads to pressure fluctuations for the original momentum interpolation scheme. This deviation is not observed with the modified scheme. CAVITY LID (Re=100) 1  0 0 0 NO SLIP WALL x Figure 421. Schematic of Cavity flow grid along with boundary conditions Modified RhieChow scheme Steady  dt=1  dt=0.1 S dt=0.01 x Original RhieChow scheme x Modified RhieChow scheme x x Figure 422. Velocity and pressure contours for cavity flow at Re=100 using different momentum interpolation schemes for various time step sizes at y=0.5 location in the cavity. Next, the results for the various momentum interpolation schemes on a curvilinear mesh using unsteady flow around a cylinder example are shown. The overall computational domain along with the grid is shown in Figure 423. The cylinder is placed in the middle of a circular domain as can be seen in the figure. An orthogonal Ogrid with increasing grid density towards the cylinder is used for this test case. There are 201 points in the circumferential direction and 121 points in the radial direction. Unsteady computations were performed until steady state is achieved, which occurs at a non dimensional time of about 250. Two different time step sizes were chosen (0.1 and 0.01) Original RhieChowscheme 