COMPUTATIONAL STUDY OF THE EFFECT OF REYNOLDS NUMBER AND MOTION TRAJECTORY ASYMMETRY ON THE AERODYNAMICS OF A PITCHING AIRFOIL AT LOW REYNOLDS NUMBER By Patrick R. Hammer          A DISSERTATION  Submitted to Michigan State University in partial fulfillment of the requirements for the degree of  Mechanical EngineeringÑDoctor of Philosophy 2016   ABSTRACT COMPUTATIONAL STUDY OF THE EFFECT OF REYNOLDS NUMBER AND MOTION TRAJECTORY ASYMMETRY ON THE AERODYNAMICS OF A PITCHING AIRFOIL AT LOW REYNOLDS NUMBER By Patrick R. Hammer  It is well established that natural flyers flap their wings to sustain flight due to poor performance of steady wing aerodynamics at low Reynolds number.  Natural flyers also benefit from the propulsive force generated by flapping.  Unsteady airfoils allow for simplified study of flapping wing aerodynamics.  Limited previous work has suggested that both the Reynolds number and motion trajectory asymmetry play a non-negligible role in the resulting forces and wake structure of an oscillating airfoil.  In this work, computations are performed to on this topic for a NACA 0012 airfoil purely pitching about its quarter-chord point.   Two-dimensional computations are undertaken using the high-order, extensively validated FDL3DI Navier-Strokes solver developed at Wright-Patterson Air Force Base.  The Reynolds number range of this study is 2,000-22,000, reduced frequencies as high as 16 are considered, and the pitching amplitude varies from 2¡ to 10¡.  In order to simulate the incompressible limit with the current compressible solver, freestream Mach numbers as low as 0.005 are used.  The wake structure is accurately resolved using an overset grid approach.   The results show that the streamwise force depends on Reynolds number such that the drag-to-thrust crossover reduced frequency decreases with increasing Reynolds number at a given amplitude.  As the amplitude increases, the crossover reduced frequency decreases at a given Reynolds number.  The crossover frequency data show good collapse for all pitching amplitudes considered when expressed as the Strouhal number based on trailing edge-amplitude  for different Reynolds numbers.  Appropriate scaling causes the thrust data to become nearly independent of Reynolds number and amplitude.  An increase in propulsive efficiency is observed as the Reynolds number increases while less dependence is seen in the peak-to-peak lift and drag amplitudes.   Reynolds number dependence is also seen for the wake structure.  The crossover reduced frequency to produce a switch in the wake vortex configuration from von K⁄rm⁄n (drag) to reverse von K⁄rm⁄n (thrust) patterns decreases as the Reynolds number increases.  As the pitching amplitude increases, more complex structures form in the wake, particularly at the higher Reynolds numbers considered.  Although both the transverse and streamwise spacing depend on amplitude, the vortex array aspect ratio is nearly amplitude independent for each Reynolds number.  Motion trajectory asymmetry produces a non-zero average lift and a decrease in average drag.  Decomposition of the lift demonstrates that the majority of the average lift is a result of the component from average vortex (circulatory) lift.  The average lift is positive at low reduced frequency, but as the reduced frequency increases at a given motion asymmetry, an increasing amount of negative lift occurs over a greater portion of the oscillation cycle, and eventually causes a switch in the sign of the lift.  The maximum value, minimum value, and peak-to-peak amplitude of the lift and drag increase with increasing reduced frequency and asymmetry.   The wake structure becomes complex with an asymmetric motion trajectory.  A faster pitch-up produces a single positive vortex and one or more negative vortices, the number of which depends on the reduced frequency and asymmetry.  When the airfoil motion trajectory is asymmetric, the vortex trajectories and properties in the wake exhibit asymmetric behavior.                          Copyright by PATRICK R. HAMMER 2016  v                      I dedicate this document to my loving wife, Katherine,  and precious daughter, Anastasia,  who helped me navigate this journey  and for all the happiness they have brought to my life.    vi ACKNOWLEDGMENTS   I first, and most, wanted to thank my wife, Katherine, and daughter, Anastasia, for their love and support during this process.  Thank you to my dad, Chris, my mom, Debra, brother, Chris, and sister, Laura, for their continued love and support during my doctoral studies, as well as life in general.  Thank you also to my family-in-laws, extended family, and friends.  My cat, Sparta, was also a great source of moral support. I would like to thank my advisors, Dr. Manoochehr Koochesfahani and Dr. Ahmed Naguib for their encouragement and patience as I tried to navigate this computational work in a world of experiments.  Thank you to my AFRL advisor, Dr. Miguel Visbal, for allowing me to use FDL3DI in my research and acting as a computational mediator with my advisors.  I also appreciate the feedback from my additional PhD. committee members, Dr. Farhad Jaberi and Dr. Chang Wang. I owe a special thanks to my colleagues at WPAFB, in particular Dr. Daniel Garmann, Dr. Caleb Barnes, and Dr. Scott Sherer, for all their assistance helping me understand the FDL3DI and overset grids, extending their tools for my use, and general brain-picking.   I want to thank the past and current members of TMUAL who acted as technical and moral support in the trenches.  In particular, thank you to Dr. Shahram Pouya, David Olson, Rohit Nehe, Fariborz Daneshvar, Nic Roberts, and Alireza Safaripour for their moral and technical support.  I extend a special gratitude in particular to David for our enjoyable conversations and debates, daily trips to SpartyÕs, and technical assistance.    vii I would like to acknowledge the MSU High Performance Computing Center (HPCC) for providing the computational resourced necessary to do this work.  I would finally like to acknowledge the support of the U.S. Air Force, which provided funding through AFOSR grant FA9550-10-1-0342.  viii TABLE OF CONTENTS   LIST OF TABLES ........................................................................................................................ x!LIST OF FIGURES ..................................................................................................................... xi!KEY TO ABBREVIATIONS .................................................................................................. xxv!CHAPTER 1. INTRODUCTION ................................................................................................ 1!1.1.!Background .................................................................................................................................1!1.2.!Reynolds Number Dependency .................................................................................................9!1.3.!Kinematic Trajectory Effects ..................................................................................................14!1.4.!Present Research ......................................................................................................................17!CHAPTER 2.  COMPUTATIONAL CONSIDERATIONS ................................................... 20!2.1.  Geometry and Trajectory Definition ...........................................................................................20!2.2.  Computational Method .................................................................................................................24!2.2.1.  Governing Equations ...............................................................................................................24!2.2.2.  The Time Marching Scheme .....................................................................................................27!2.2.3.  Spatial Discretization and Filter ..............................................................................................29!2.3.  Computational Domain .................................................................................................................31!2.4.  Initial Conditions and Boundary Conditions ..............................................................................35!2.5.  Convergence Studies .....................................................................................................................36!2.5.1.  Convergence Study Summary ...................................................................................................36!2.5.2.  Grid, Time-step, Subiteration Studies ......................................................................................36!2.5.3.  Effect of Freestream Mach Number .........................................................................................38!2.6.  Validation .......................................................................................................................................43!2.6.1.  Validation of Sinusoidal Pitching ............................................................................................44!2.6.2.  Validation of Asymmetric Pitching ..........................................................................................51!2.7.  Final Computational Considerations ...........................................................................................54!2.7.1.  Time-Periodic Convergence .....................................................................................................54!2.7.2.  Grid Remapping .......................................................................................................................54!CHAPTER 3.  EFFECT OF REYNOLDS NUMBER ON AERODYNAMIC LOADS ...... 55!3.1.!Average Thrust Coefficient .....................................................................................................55!3.2.!Thrust Crossover Frequency ...................................................................................................63!3.3.!Propulsive Efficiency ................................................................................................................69!3.4.!Force Fluctuation .....................................................................................................................71!3.5.!Comparison with Inviscid Linear Theory ..............................................................................78!CHAPTER 4.  EFFECT OF REYNOLDS NUMBER ON WAKE STRUCTURE .............. 81!4.1.!Flow Near Trailing Edge .........................................................................................................81!4.2.!Wake Structure Configuration for 2¡ Pitching .....................................................................87!4.3.!Transverse Spacing Crossover Frequency ...........................................................................100! ix 4.4.!Effect of Pitching Amplitude .................................................................................................102!CHAPTER 5.  EFFECT OF TRAJECTORY ASYMMETRY ON AERODYNAMIC LOADS ....................................................................................................................................... 109!5.1.!Average Forces ........................................................................................................................111!5.2.!Lift Decomposition .................................................................................................................117!5.3.!Force Fluctuation ...................................................................................................................130!CHAPTER 6.  EFFECT OF TRAJECTORY ASYMMETRY ON WAKE STRUCTURE..................................................................................................................................................... 133!6.1.!Wake Structure .......................................................................................................................133!6.2.!Vortex Formation at Trailing Edge ......................................................................................138!6.3.!Vortex Properties ....................................................................................................................148!CHAPTER 7.  CONCLUSIONS AND FUTURE RECOMMENDATIONS ....................... 152!7.1.!Reynolds Number Dependence .............................................................................................153!7.2.!Trajectory Asymmetry ...........................................................................................................155!7.3.!Future Recommendations ......................................................................................................155!APPENDICES ........................................................................................................................... 157!APPENDIX A.  SMOOTHED ASYMMETRIC MOTION TRAJECTORY ................................158!APPENDIX B.  COMPUTATIONAL METHOD ............................................................................165!APPENDIX C.  CONVERGENCE STUDIES ..................................................................................198!APPENDIX D.  ADDITIONAL VALIDATION ..............................................................................210!APPENDIX E.  SURFACE PRESSURE AND SHEAR STRESS DISTRIBUTIONS ..................217!APPENDIX F.  LINEAR THEORY ..................................................................................................222!APPENDIX G.  ADDITIONAL VORTEX PROPERTY RESULTS .............................................226!APPENDIX H.  ADDITIONAL PITCHING AMPLITUDE FLOW VISUALIZATION RESULTS .............................................................................................................................................232!APPENDIX I.  ADDITIONAL ASYMMTRIC MOTION TRAJECTORY FLOW VISUALIZATION RESULTS ...........................................................................................................235!BIBLIOGRAPHY ..................................................................................................................... 246! x LIST OF TABLES   Table 2.1.  Final simulation parameters. ...................................................................................... 36!Table 2.2.  Spacing information for three grids of varying resolution. ........................................ 37!Table 2.3.  Vortex transverse spacing crossover reduced frequency for Re ! 12,000 and !o = 2¡. ............................................................................................................................................... 50 Table B.1.  Coefficients for the interior finite difference schemes. ........................................... 180!Table B.2.  Coefficients for boundary point 1. ........................................................................... 183!Table B.3.  Coefficients for boundary point 2. ........................................................................... 183!Table B.4.  Interior point filter coefficients. .............................................................................. 187!Table B.5.  Filter scheme utilized in current work. .................................................................... 191 Table C.1.  Mach number effect on average thrust coefficient for different cases. ................... 203!Table C.2.  Average lift coefficient for different asymmetries. ................................................. 209!Table C.3.  Average thrust coefficient for different asymmetries. ............................................ 209 Table D.1.  Vortex properties for P1 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%. .................................................................................. 212!Table D.2.  Vortex properties for N1 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%. .................................................................................. 212!Table D.3.  Vortex properties for N2 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%. .................................................................................. 212!    xi LIST OF FIGURES   Figure 1.1.  a) Operating Reynolds numbers for different flyers (both natural and man-made) b) maximum lift-to-drag ratio versus Reynolds number (Lissaman, 1983).  The broken line corresponds to the Reynolds number of 105, where there is a sharp drop in performance.  Figure courtesy of Annual Reviews, Inc. ................................................................................ 2!Figure 1.2. Schematic of flat plate undergoing a) pitching and b) plunging. ................................ 3!Figure 1.3.  Illustration of vorticity shed into the wake and its relationship to the airfoil circulation (von K⁄rm⁄n and Burgers, 1935), with resulting average velocity behavior in the wake based rotation of vorticity sheet.  Figure courtesy of Springer-Verlag. ........................ 4!Figure 1.4. Flow visualization of traditional von K⁄rm⁄n street, aligned vortices, and reverse von K⁄rm⁄n street with accompanying velocity profile. ........................................................ 5!Figure 1.5.  Thrust coefficient dependency on reduced frequency and pitching amplitude (Koochesfahani, 1989).  Flow visualization of vortical patterns from current computations are included as reference.  Figure courtesy of AIAA. ............................................................ 6!Figure 1.6.  Schematic of vortex streamwise (a) and transverse (b) spacing. ............................... 8!Figure 1.7.  Thrust coefficient vs. reduced frequency (left), crossover reduced frequency vs. Reynolds number for NACA 0012 pitching about its quarter-chord with a 2¡ amplitude at a Reynolds number of approximately 12,000 (right). .............................................................. 10!Figure 1.8.  Propulsive efficiency versus Strouhal number (Liu and Kawachi, 1999).  Figure courtesy of Elsevier. .............................................................................................................. 13!Figure 1.9.  a) Schematic of airfoil and kinematics considered and b) oscillation amplitude vs. Strouhal number.  The chord Reynolds number is approximately 1,200.  The blue line in the right image represents the switch from traditional von K⁄rm⁄n street to reverse von K⁄rm⁄n vortex patterns.  Figure courtesy of American Physical Society. ......................................... 14!Figure 1.10.  Schematic of asymmetric trajectory definition (Koochesfahani, 1989).  Figure courtesy of AIAA. ................................................................................................................. 15!Figure 1.11.  Flow visualization for three asymmetries: S = 38%, 50%, and 61% (Koochesfahani, 1989).  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  Figure courtesy of AIAA. ..................................................................................................... 16!Figure 1.12.  Schematic of investigated flow with important parameters labeled. ...................... 18  xii Figure 2.1.  NACA 0012 geometry with rounded trailing edge. ................................................. 20!Figure 2.2.  Angle of attack time-histories for S = 38%, 50%, and 62%, normalized by the pitching amplitude, !o. .......................................................................................................... 22!Figure 2.3.  Pitching acceleration and lift time-history for S = 38%, illustrating influence of ". 23!Figure 2.4.  Instantaneous spanwise vorticity fields, #zc/U$, using only O-grid.  The Reynolds number is 12,000, reduced frequency is 6.68, and 2¡ pitching amplitude. ........................... 31!Figure 2.5.  Overset computational domain shown at three different scales. .............................. 33!Figure 2.6.  Comparison between the instantaneous spanwise vorticity fields, #zc/U$, using O-grid (top) and overset grid (bottom).  The Reynolds number is 12,000, reduced frequency is 6.68, and 2¡ pitching amplitude. ........................................................................................... 34!Figure 2.7.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at "/T = 0.50 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.005, k = 6.68, !o = 2¡, and S = 30%. ................................................................................. 38!Figure 2.8.  Mach number effects on the peak-to-peak amplitude of lift and drag coefficients, and average thrust coefficient for Reynolds number of 12,000 and 2¡ pitching amplitude. 40!Figure 2.9.  Lift history for S = 38% illustrating effect of decreasing freestream Mach number for Re = 12,000, k = 6.68, and !o = 2¡. ................................................................................. 41!Figure 2.10.  Spanwise vorticity fields, #zc/U$, S = 38% at %/T = 0.0, illustrating influence of M$.  Dye-visualization from the experiment by Koochesfahani (1989) is included for reference at the bottom.  Bottom figure courtesy of AIAA. ................................................. 42!Figure 2.11.  Density variation towards the incompressible limit with decreasing freestream Mach number. ....................................................................................................................... 43!Figure 2.12.  Peak-to-peak lift coefficient and peak-to-peak drag coefficient for 2¡ and 4¡ pitching amplitudes at a Reynolds number of approximately 12,000.  Data from the literature is included for comparison. ................................................................................... 44!Figure 2.13.  Average thrust coefficient for 2¡ and 4¡ pitching amplitudes at a Reynolds number of approximately 12,000.  Data from the literature is included for comparison. .................. 46!Figure 2.14.  Instantaneous spanwise vorticity field (#z/cU$) comparison at a reduced frequency of approximately 5.2 (left column) and 12 (right column) between current compressible computations with experiments by Bohl and Koochesfahani (2009) and incompressible computations by Jee and Moser (2012).  The Reynolds number is approximately 12,000 and the pitching amplitude is 2¡ amplitude.  The airfoil is at zero angle of attack and is pitching down.  Top and middle figures courtesy of Cambridge University Press and Elsevier. ...... 47! xiii Figure 2.15.  Vortex properties versus reduced frequency for Reynolds number of 12,000 and 2¡ amplitude.  Data from available literature is included for comparison.  Properties calculated at x/c = 0.5; except for Naguib, et al., (2011), which are done at x/c = 1.0. ......................... 50!Figure 2.16.  Instantaneous spanwise vorticity field #z/cU$ comparison between incompressible dye-visualization experiment (top) and compressible computation (bottom) at a reduced frequency of 6.68 and S = 38%.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  Top image courtesy of AIAA. .................................................................................... 52!Figure 2.17.  Average and r.m.s. velocity profiles measured at x/c = 1.0 using analytical pitching trajectory with S = 38% and k = 6.68. Trajectory comparison between analytical formula and experimental data is included on top for reference. ....................................................... 53 Figure 3.1.  Total, pressure-induced, and friction-induced average drag versus Reynolds number for static NACA 0012 at ! = 0¡. The broken lines represent a least squares power-law fit to the data. ................................................................................................................................. 56!Figure 3.2.  Average thrust coefficient versus reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡. ............................................................................... 57!Figure 3.3.  Average thrust coefficient, with static drag removed, versus reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡. ................................................ 58!Figure 3.4.  Average thrust coefficient components, with static drag components removed, versus reduced frequency for different Reynolds numbers (color).  The pitching amplitude is 2¡. ...................................................................................................................................... 60!Figure 3.5.  Average thrust coefficient versus reduced frequency for different pitching amplitudes (symbol).  Data for Re = 2,000 on left and for Re = 22,000 on right. ................ 61!Figure 3.6.  Average thrust coefficient versus Strouhal number for different pitching amplitudes (symbols).  Data for Re = 2,000 on left and for Re = 22,000 on right. ................................. 62!Figure 3.7.  Instantaneous spanwise vorticity field (#zc/U$) for two reduced frequencies (1.0 and 4.0) and two Reynolds numbers (2,000 and 22,000).  The pitching amplitude is 10¡ and the airfoil is at its maximum angle of attack position and is about to pitch down (%/T = 0.25). . 63!Figure 3.8.  Thrust crossover reduced frequency versus Reynolds number for different pitching amplitudes (symbols). ........................................................................................................... 64!Figure 3.9.  Thrust crossover Strouhal number versus Reynolds number for different pitching amplitudes (symbols). ........................................................................................................... 65!Figure 3.10.  Thrust crossover Strouhal number versus Reynolds number data comparison with literature.  Only the 2¡, 6¡, and 10¡ from the current data set are included.  Note that Godoy-Diana, et al., (2008) used a tear-drop (TD) shaped airfoil. ................................................... 66! xiv Figure 3.11.  Average thrust coefficient versus Strouhal number for different Reynolds numbers (color) and pitching amplitudes (symbols).  Only the 2¡, 6¡, and 10¡ are included. ............ 68!Figure 3.12.  Average thrust coefficient versus the Strouhal number scaled by the crossover Strouhal number as a function of Reynolds number (color) and pitching amplitude (symbol).  The thrust coefficient is scaled by the static airfoil drag. .................................... 69!Figure 3.13.  Average propulsive efficiency versus Strouhal number as a function of Reynolds number (color) and pitching amplitude (symbol). ................................................................ 71!Figure 3.14.  Reynolds number effect of the fluctuating force characteristics for static NACA 0012 airfoil at ! = 0¡. ........................................................................................................... 72!Figure 3.15.  Lift and drag fluctuations for different reduced frequencies.  The Reynolds numbers are Re = 2,000 (blue), 7,000 (red), and 22,000 (yellow).  The pitching amplitude is 2¡ amplitude.  Note that the scale is different for the different plots. ................................... 73!Figure 3.16.  Magnitude of lift fluctuation FFT for the Re = 22,000 case.  Instantaneous spanwise vorticity field at %/T = 0.25 is included for visual reference. ................................. 74!Figure 3.17.  Schematic of the phase angle definition. ................................................................ 75!Figure 3.18.  Peak-to-peak lift amplitude (left), peak-to-peak drag amplitude (middle), and lift phase lag (right).   Different Reynolds numbers are shown with the symbol color and different pitching amplitudes are shown with symbol type. ................................................. 76!Figure 3.19.  Peak-to-peak lift amplitude normalized by !" (left) and peak-to-peak drag amplitude normalized by !"# (right).   Different Reynolds numbers are shown with the color and different pitching amplitudes are shown with symbols. ....................................... 77!Figure 3.20.  Lift (top) and drag (bottom) histories, showing contribution from pressure and shear stress distribution compared to the total force for Re = 2,000 (left) and Re = 22,000 (right) at k = 6.0 and !o = 2¡. ................................................................................................ 77!Figure 3.21.  Average thrust coefficient with static drag removed versus reduced frequency for different Reynolds numbers (color) and pitching amplitudes (symbol).  Data shown is for 2¡ and 10¡ amplitudes.  The linear theory prediction is included for comparison. ................... 78!Figure 3.22.  Average thrust coefficient from pressure with static pressure drag removed versus reduced frequency for different Reynolds numbers (color) and pitching amplitudes (symbol).  Data shown is for 2¡ and 10¡ amplitudes.  The linear theory prediction is included for comparison. ...................................................................................................... 79!Figure 3.23.  Peak-to-peak lift amplitude (left), peak-to-peak drag amplitude (middle), and lift phase lag (right).   Different Reynolds numbers are shown with the symbol color and  xv different pitching amplitudes are shown with symbol type.   The linear theory prediction of Theodorsen (1935) and Garrick (1936) is included for reference. ....................................... 80 Figure 4.1.  Upper surface u velocity and z-vorticity profiles measured at trailing edge for static airfoil at ! = 0¡ for different Reynolds numbers. ................................................................. 82!Figure 4.2.  Boundary layer thickness as a function of Reynolds number for static NACA 0012 at ! = 0¡. ............................................................................................................................... 83!Figure 4.3.  Streamwise velocity boundary profiles, where Y is scaled by the static airfoil BLT, at eight phases of the oscillation cycle for different of Reynolds numbers.  The reduced frequency 6.0 and pitching amplitude is 2¡.  Measurements are made at the trailing edge on the upper surface.  The coordinate system starts at the airfoil surface.  The airfoil position at every other phase is included for reference.  The cyan line denotes u = 0 and the magenta line represents u = U$. .......................................................................................................... 85!Figure 4.4.  Instantaneous spanwise vorticity field (#zc/U$) at 8 phases of the oscillation for k = 6.0, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The green line is the measurement location and the white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075...................................................................................................................................... 86!Figure 4.5.  Instantaneous spanwise vorticity fields (#zc/U$) for static NACA 0012 at ! = 0¡ and different Reynolds numbers (increasing from the top down). ....................................... 88!Figure 4.6.  Wake vortex transverse (left) and streamwise (right) spacing dependence on Reynolds number for static NACA 0012 at ! = 0¡. ............................................................. 89!Figure 4.7.  Instantaneous spanwise vorticity field (#zc/U$) for !o = 2¡ and different Reynolds numbers (increasing left to right) and reduced frequencies (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.5 " y/c " 0.5. .................................................... 90!Figure 4.8.  Instantaneous spanwise vorticity field (#zc/U$) for k = 6.0, !o = 2¡, and different Reynolds numbers (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ...................................................................................................... 92!Figure 4.9.  Vorticity profiles normalized by peak vorticity at x/c ! 1.0 for different Reynolds numbers at k = 6.0 and !o = 2¡. ............................................................................................. 92!Figure 4.10.  Schematic of $# variables for a Gaussian vortex. .................................................. 93! xvi Figure 4.11.  Instantaneous spanwise vorticity field (#zc/U$) for k = 6.0, !o = 2¡, and different Reynolds numbers, showing the spatial radius (increasing from the top down). ................. 95!Figure 4.12.  Positive and negative vortex trajectories for different Reynolds numbers at k = 6.0 and !o = 2¡. ........................................................................................................................... 96!Figure 4.13.  Wake vortex transverse spacing (top) and streamwise spacing (bottom) at k = 6.0 and !o = 2¡ for different Reynolds numbers. ........................................................................ 98!Figure 4.14.  Wake vortex transverse spacing (top) and streamwise spacing (bottom) dependency on k and Re for !o = 2¡.  Measurements were done at x/c = 1.0. ...................... 99!Figure 4.15.  Crossover reduced frequency for b = 0 dependence on Reynolds number for !o = 2¡.  Data from literature are included for reference.  Note that Monnier, et al., (2015) used a modified NACA series airfoil. ............................................................................................ 101!Figure 4.16.  Instantaneous spanwise vorticity field (#zc/U$) for visually aligned vortex array and !o = 2¡ for different Reynolds numbers (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ............................................................. 101!Figure 4.17.  Instantaneous spanwise vorticity field (#zc/U$) at k = 3.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4 ................................................... 103!Figure 4.18.  Wake vortex transverse (top) and streamwise (bottom) spacing dependence on k, for Re = 2,000 (left) and 22,000 (right) and !o = 2¡-10¡.  Data is presented for only cases with well-formed single vortex pairs. ................................................................................. 104!Figure 4.19.  Wake vortex transverse (top) and streamwise (bottom) spacing dependence on k, for Re = 2,000 (left) and 22,000 (right) and !o = 2¡-10¡.  Data is presented for only cases with well-formed single vortex pairs. ................................................................................. 105!Figure 4.20.  Wake vortex transverse spacing switch crossover reduced frequency (left) and Strouhal number (right) dependence on Reynolds number for different amplitudes.  Measurements from computations done at x/c = 1.0.  Data from literature are included for reference.  Note that Godoy-Diana, et al., (2008) used for a tear-drop shaped airfoil and Monnier, et al., (2015) used a modified NACA series airfoil. ........................................... 106!Figure 4.21.  Wake vortex configuration aspect ratio (b/a) dependence on St for different Reynolds number and pitching amplitudes.  Measurements from computations done at x/c = 1.0........................................................................................................................................ 107! xvii Figure 4.22.  Wake vortex configuration aspect ratio (b/a) dependence on St for different Reynolds number for !o = 2¡ (circles) and 10¡ (squares).  Measurements from computations done at x/c = 1.0. ................................................................................................................. 108 Figure 5.1.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) for different values of S. ...................................................................................................... 110!Figure 5.2.  CL (left) and CT (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. ............................................................ 112!Figure 5.3.  Thrust differential vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. ............................................................................ 113!Figure 5.4.  Pressure (top) and shear stress (bottom) components of average CL (left) and CT (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. ..................................................................................................... 114!Figure 5.5.  Schematic illustrating lift averaging definitions for the example of k = 6.68 and S = 38%.  The angle of attack and pitching acceleration are included for reference. ............... 115!Figure 5.6.  Ratio of cycle time-average of positive to negative lift (left), ratio of the duration of positive to that of negative lift period (middle), and ratio of duration time-average of positive to negative lift (right). ........................................................................................... 117!Figure 5.7.  Ratio of lift circulatory and non-circulatory components to total lift at the phase of maximum lift.  The lift coefficient time history for k = 6.68 (indicated by the blue line) using TheodorsenÕs formula is included in the figure for reference, where the phase of maximum lift is indicated by the green dashed line. .......................................................... 118!Figure 5.8.  Schematic of fluid domain, Vf, with control surface # about solid body B. ........... 119!Figure 5.9.  Schematic of integration domain for computing the vortex lift. ............................ 121!Figure 5.10.  Time history of angle of attack (top left), acceleration (top right), vortex lift coefficient (bottom left), and total lift coefficient (bottom right) for k = 6.68 and S = 38%.............................................................................................................................................. 122!Figure 5.11.  CL vs. CVL for each reduced frequency and asymmetry.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. ............................................................................ 123!Figure 5.12.  The time histories of lift, integrated Lamb vector contribution, and acceleration from computations by Wang, et al., (2013) for a flat plate at non-zero angle of attack undergoing combined pitching and plunging at Re = 300.  Instantaneous fields of spanwise vorticity (top), Lamb vector projected in the vertical direction (middle), and acceleration projected in the vertical direction (bottom) at % = 18.8 (indicated by the red dashed line) are included on the right for reference. Figure courtesy of AIP Publishing, LLC. .................. 124! xviii Figure 5.13.  Angle of attack (top) and vortex lift coefficient (bottom) time histories for k = 4 (left) and 8 (right). ............................................................................................................... 125!Figure 5.14.  Schematic illustrating vortex lift averaging definitions for the example of k = 6.68 and S = 38%.  The angle of attack is included for reference. ............................................. 126!Figure 5.15.  Ratio of the duration of positive to that of negative lift period (left), and ratio of duration time-average of positive to negative lift (right). ................................................... 128!Figure 5.16.  Vortex lift coefficient time histories for k = 4 (left) and 8 (right), with S = 38%. 129!Figure 5.17.  Lift and drag fluctuation for different asymmetries at k = 6.68.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  The angle of attack and pitching acceleration history is included for reference. .................................................................... 131!Figure 5.18.  Maximum (top), minimum (middle), and peak-to-peak (bottom) values of CL (left) and CD (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. ............................................................................................... 132 Figure 6.1.  Instantaneous spanwise vorticity field (#zc/U$) for Re = 12,000 and !o = 2¡ for different reduced frequencies (increasing to the right) and asymmetries (increasing from top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The range covered in the figure is 0.5 " X/c " 3 and -0.3 " Y/c " 0.3. ................................................ 135!Figure 6.2.  Schematic of vortex labeling for both S = 50% and S = 38%. ............................... 136!Figure 6.3.  Positive and negative vortex centroid trajectories at three reduced frequencies numbers for S = 50% and 38%.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. .................................................................................................................................... 137!Figure 6.4.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 6.68, S = 50%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 140!Figure 6.5.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 6.68, S = 50%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 141!Figure 6.6.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 6.68, S = 38%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which  xix the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 143!Figure 6.7.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 6.68, S = 38%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 144!Figure 6.8.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 6.68, S = 30%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 146!Figure 6.9.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 6.68, S = 30%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 147!Figure 6.10.  Schematic of vortex labeling for both S = 50% and S = 38% for vortex properties............................................................................................................................................. 148!Figure 6.11.  Vortex properties of positive vortex P (circles) and negative vortex N (triangles) for different values of S at different reduced frequencies.  The Reynolds numbers is 12,000 and the pitching amplitude is 2¡. ........................................................................................ 150!Figure 6.12.  Schematic of vortex spacing definitions for asymmetry trajectory. ..................... 151!Figure 6.13.  Streamwise spacing of centroid between opposite sign vortices (top) and same-sign vortices (bottom) for different values of S at different reduced frequencies.  The Reynolds numbers is 12,000 and the pitching amplitude is 2¡. .......................................................... 151 Figure A.1.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) coefficients for k = 6.68 and S = 38%, illustrating number of Fourier series terms. .......... 159!Figure A.2.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) coefficients for k = 6.68 and S = 38%, illustrating influence of ". ..................................... 162!Figure A.3.  Time history of lift lift (left) and drag (right) coefficients for k = 6.68 and S = 38%, illustrating influence of ". ................................................................................................... 163!Figure A.4.  Average lift (left) and thrust (right) coefficients for k = 6.68 and different asymmetries, illustrating influence of ". ............................................................................. 163! xx Figure A.5.  Average u (left), r.m.s. u (middle), and r.m.s. v (right) for k = 6.68 and S = 38%. 164 Figure B.1.  Notation for the interior and boundary points of a 1-D mesh. ............................... 177!Figure B.2.  Notation for the interior point stencil of a 1-D mesh. ............................................ 178!Figure B.3.  Modified wavenumber versus scaled wavenumber for different finite difference schemes. .............................................................................................................................. 182!Figure B.4.  Effect of filter order (left) and &f (right) on filter characteristics. ......................... 188!Figure B.5.  Resolution comparison between the compact scheme C6 and filter F8 (%&'()*().............................................................................................................................................. 189!Figure B.6.  Spectral filter response for the one-sided boundary point 2 formula for various filter orders and %&. ...................................................................................................................... 190!Figure B.7.  Spectral filter response for boundary point 2 for a one-sided filter (F4) and centered filter (F8) with similar spectral characteristics. .................................................................. 191!Figure B.8.  Schematic of donors points (from Mesh 1) and receiver points (from Mesh 2). ... 192!Figure B.9.  Illustration of block decomposition with five-point overlap. ................................ 195 Figure C.1.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at "/T = 0.625 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.015, k = 11.9, and !o = 2¡, and S = 50%. ......................................................................... 199!Figure C.2.  Instantaneous vorticity measured at x/c = 1.0 at "/T = 0.25 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.015, k = 11.9, and !o = 2¡, and S = 50%. 199!Figure C.3.  Lift time histories (top) and r.m.s vorticity profiles (bottom) measured at x/c = 1.0 for different symmetric pitching cases illustrating solution convergence by Mach number.............................................................................................................................................. 201!Figure C.4.  Error percentage for lift and drag histories versus Mach number for sinusoidal pitching cases. ..................................................................................................................... 202!Figure C.5.  Error percentage for lift and drag maximum and minimum values versus Mach number for sinusoidal pitching cases. ................................................................................. 202!Figure C.6.  Error percentage for r.m.s. vorticity profiles versus Mach number for sinusoidal pitching cases. ..................................................................................................................... 203! xxi Figure C.7.  Instantaneous surface pressure coefficient distributions for k ! 12 at four phases illustrating effect of Mach number. .................................................................................... 204!Figure C.8.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at "/T = 0.50 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.005, k = 6.68, and !o = 2¡, and S = 30%. ......................................................................... 205!Figure C.9.  Instantaneous vorticity profiles measured at x/c = 1.0 at "/T = 0.25 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.005, k = 6.68, and !o = 2¡, and S = 30%. .................................................................................................................................... 206!Figure C.10.  Lift time histories (top) and fluctuation vorticity profiles (bottom) measured at x/c = 1.0 for different asymmetric pitching cases illustrating solution convergence by Mach number. ............................................................................................................................... 207!Figure C.11.  Error percentage for lift and drag histories versus Mach number for sinusoidal pitching cases. ..................................................................................................................... 208!Figure C.12.  Error percentage for lift and drag maximum and minimum values by Mach number for asymmetric pitching cases. .............................................................................. 208!Figure C.13.  Error percentage for r.m.s. vorticity profiles versus Mach number for sinusoidal pitching cases. ..................................................................................................................... 209 Figure D.1.  Average u (left), r.m.s. u (middle), and r.m.s. v (right) profiles vs. y/c for Reynolds number of approximately 12,000 and 2¡ amplitude.  The reduced frequency is approximately 12.  Data from Bohl and Koochesfahani (2009) is included for comparison.  Profiles are measured at x/c = 1.0 for both the computational and experimental results. .. 210!Figure D.2.  Comparison of the downstream decay of peak vorticity between computation and theory.  The initial condition location for the theory is highlighted by the yellow circle.  Flow conditions: Re = 12,000, M$ = 0.015, k = 11.9, and !o = 2¡, and S = 50%. .............. 211!Figure D.3.  Schematic of vortex labeling for both S = 50% and S = 38% for vortex properties.............................................................................................................................................. 213!Figure D.4.  Comparison between the experimental asymmetric trajectory and smoothed trajectory for S = 38%. ........................................................................................................ 213!Figure D.5.  Average and r.m.s. velocity profile measured at x/c = 1.0 using actual experimental pitching trajectory for Re = 12,000 at S ! 38%., illustrating effect of k. ............................ 215! xxii Figure D.6.  Wake vortex trajectories for k = 6.20 (left), 6.40 (middle), and 6.68 (right) using actual experimental pitching trajectory for Re = 12,000 at S ! 38%. ................................. 215 Figure E.1.  Average pressure coefficient (top) and skin friction coefficient (bottom) for static NACA 0012 at ! = 0¡ for different Reynolds numbers. ..................................................... 218!Figure E.2.  Average pressure and skin friction coefficient for a Reynolds number of 2,000 (left) and 22,000 (right) oscillating at different k.  The pitching amplitude is 2¡.  The airfoil is shown at the bottom for spatial reference. .......................................................................... 219!Figure E.3.  Average projected pressure and skin friction coefficient for a Reynolds number of 2,000 (left) and 22,000 (right) oscillating at different k.  The pitching amplitude is 2¡.  The airfoil is shown at the bottom for spatial reference. ............................................................ 220 Figure F.1.  Schematic of oscillating plate (Theodorsen, 1935; Garrick, 1936). ...................... 223!Figure F.2.  F and G vs. 1/k (Theodorsen, 1935; Garrick, 1936). ............................................. 224 Figure G.1.  Vortex properties as a function of reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡. ............................................................................. 227!Figure G.2.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 2,000. .................................................................................................. 228!Figure G.3.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 7,000. .................................................................................................. 228!Figure G.4.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 12,000. ................................................................................................ 229!Figure G.5.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 22,000. ................................................................................................ 229!Figure G.6.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 2,000. .................................................................................................. 230!Figure G.7.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 7,000. .................................................................................................. 230!Figure G.8.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 12,000. ................................................................................................ 231!Figure G.9.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 22,000. ................................................................................................ 231 Figure H.1.  Instantaneous spanwise vorticity field (#zc/U$) at k = 2.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).   xxiii The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 232!Figure H.2.  Instantaneous spanwise vorticity field (#zc/U$) at k = 2.67 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 232!Figure H.3.  Instantaneous spanwise vorticity field (#zc/U$) at k = 3.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 233!Figure H.4.  Instantaneous spanwise vorticity field (#zc/U$) at k = 3.3 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 233!Figure H.5.  Instantaneous spanwise vorticity field (#zc/U$) at k = 4.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 234!Figure H.6.  Instantaneous spanwise vorticity field (#zc/U$) at k = 5.2 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4. .................................................. 234 Figure I.1.  Instantaneous spanwise vorticity field (#zc/U$) for k = 4.0 at Re = 12,000 and !o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ................................................................................... 235!Figure I.2.  Instantaneous spanwise vorticity field (#zc/U$) for k = 5.20 at Re = 12,000 and !o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ................................................................................... 236!Figure I.3.  Instantaneous spanwise vorticity field (#zc/U$) for k = 5.80 at Re = 12,000 and !o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ................................................................................... 237!Figure I.4.  Instantaneous spanwise vorticity field (#zc/U$) for k = 6.68 at Re = 12,000 and !o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ................................................................................... 238!Figure I.5.  Instantaneous spanwise vorticity field (#zc/U$) for k = 8.00 at Re = 12,000 and !o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0). ................................................................................... 239! xxiv Figure I.6.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 5.20, S = 50%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 240!Figure I.7.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 5.20, S = 50%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 241!Figure I.8.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 5.20, S = 38%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 242!Figure I.9.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 5.20, S = 38%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 243!Figure I.10.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-up for k = 5.20, S = 30%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 244!Figure I.11.  Instantaneous spanwise vorticity field (#zc/U$) at 12 phases of the pitch-down for k = 5.20, S = 30%, Re = 12,000, and !o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 " X/c " 1.2 and -0.075 " Y/c " 0.075. ............................................................. 245!   xxv KEY TO ABBREVIATIONS   Abbreviations AFRL  Air Force Research Laboratory AoA  Angle of attack BLT  Boundary layer thickness CFD  Computational fluid dynamics FC  Fractional change FFT  Fast Fourier Transform GCL  Geometric Conservation Law LE  Leading edge r.m.s.  Root mean square RMSD  Root mean square difference TE  Trailing edge UnPb  Unpublished results WPAFB Wright-Pattern Air Force Base Author Abbreviations B & K  Bohl and Koochesfahani  H & L  Huang and Lin R & S  Ramamurti and Sandberg J & M  Jee and Moser M & W Mackowski and Williamson X & L  Xiao and Liao Y & L  Young and Lai    1 CHAPTER 1. INTRODUCTION  1.1.!Background Since the early 20th century, oscillating airfoils have been the focus of a large number of investigations.  Although early work focused on understanding the mechanisms related to flutter (Theodorsen, 1935), more recent studies of oscillating airfoils have been driven by developing Micro Air Vehicles (MAVs).  Fixed wings perform well to generate lift at high Reynolds numbers (Re # 106) where aircrafts operate.  As a reference, the operating Reynolds number for different natural and man-made flyers, and the corresponding maximum lift-to-drag ratio, is illustrated in Figure 1.1 (Lissaman, 1983).  However, the performance of fixed wing aerodynamics deteriorates rapidly at the lower Reynolds numbers in which MAVs typically to operate (Re = 102-105).  At Reynolds numbers of order 105, the maximum lift-to-drag drops significantly due to higher drag and lower lift experienced at the lower Reynolds numbers, where the boundary layer is still laminar.  Here, the Reynolds number is based on the airfoil chord length and is defined as the ratio of inertia to viscous forces (i.e., Re = '$U$c/µ$, where '$ is the freestream fluid density, U$ the freestream flow speed, c the airfoil chord length, and µ$ the freestream fluid dynamic viscosity).  To overcome the loss in steady aerodynamic performance at these low Reynolds numbers, natural flyers flap their wings.    2  Figure 1.1.  a) Operating Reynolds numbers for different flyers (both natural and man-made) b) maximum lift-to-drag ratio versus Reynolds number (Lissaman, 1983).  The broken line corresponds to the Reynolds number of 105, where there is a sharp drop in performance.  Figure courtesy of Annual Reviews, Inc. The flows generated by flapping motions are predominantly vortical and have been studied extensively, primarily due to their connection with the propulsion of the majority of natureÕs flyers and swimmers (e.g. Wu, 1971; Lighthill, 1975) through analytical approaches (e.g. Theodorsen, 1935; Garrick, 1936; von K⁄rm⁄n and Burgers, 1935; von K⁄rm⁄n and Sears, 1938), experiments (e.g. Katzmayr, 1922; DeLaurier and Harris, 1982; Freymuth, 1988; Koochesfahani, 1989; Jones and Platzer, 1998; Lai and Platzer, 1999), and computations (e.g. Katz and Weihs, 1978; Stanek and Visbal, 1989; Liu and Kawachi, 1999; Ramamurti and Sandberg, 2001).  The idealized form of two-dimensional flapping typically studied consists of an oscillatory action comprised of two distinct motions: cross-stream translation (plunging or heaving) and rotation about a specific axis along the chord length (pitching).  Although true flapping is more complex, the study of the simplified oscillatory motions elucidates the basic aerodynamic mechanisms that natural flyers take advantage of.  One mechanism of flapping comes from pitching causing a geometric angle of attack, +, and plunging generating an effective angle of attack, +,, respectively, between the chord line and the direction of the relative approach flow.  This is shown schematically in Figure 1.2 for a flat plate.  The combination of  3 the pitching and plunging represents the effective angle of attack experienced by flyers in idealized, two-dimensional kinematic conditions.    Figure 1.2. Schematic of flat plate undergoing a) pitching and b) plunging. It is well known from inviscid theory that an airfoil at a fixed angle of attack experiences a lift force because of the circulation around the airfoil (Anderson, Jr., 1984).  When the angle of attack continuously changes, as in the case of oscillatory motions, there is an accompanying change in circulation around the airfoil (von K⁄rm⁄n and Burgers, 1935), which causes circulation of the opposite sign as the circulation bound to the airfoil to shed into the wake due to KelvinÕs circulation theorem (Thomson, 1869) in the form of a wavy vorticity sheet (see Figure 1.3).  In the figure, the circular arrows indicate the direction of rotation in the vorticity sheet.  The resulting average velocity distribution in the wake induced by the rotation of the vorticity sheet causes a jet-like velocity profile, also shown in the figure.  Early flow visualization (Bratt, 1953) and inviscid numerical studies (Gieseng, 1968; Katz and Weihs, 1978) have shown that this vorticity sheet rolls up into distinct vortices at high enough frequencies.   4  Figure 1.3.  Illustration of vorticity shed into the wake and its relationship to the airfoil circulation (von K⁄rm⁄n and Burgers, 1935), with resulting average velocity behavior in the wake based rotation of vorticity sheet.  Figure courtesy of Springer-Verlag. At the Reynolds numbers in which natural flyers operate, it is known (Koochesfahani, 1989; Jones and Platzer, 1998) that there exist oscillation frequency requirements to generate the thrust-producing vortical pattern.  At lower oscillation frequencies, the wake structure oftentimes takes the form of a traditional von K⁄rm⁄n street (see Figure 1.4 for an example), which is known to produce drag based on work by von K⁄rm⁄n (1911, 1912) for cylinders in a cross-flow.  As the oscillation frequency increases, the vortices become aligned and lead a uniform velocity profile that has no net momentum deficit or surplus.  The vortices then switch in orientation as the oscillation frequency increases further and produce the reverse von K⁄rm⁄n street in the airfoilÕs wake, which leads to thrust due to the jet-like velocity profile.    5  Figure 1.4. Flow visualization of traditional von K⁄rm⁄n street, aligned vortices, and reverse von K⁄rm⁄n street with accompanying velocity profile. Although force measurements (DeLaurier and Harris, 1982) and qualitative deduction based on the flow patterns (Freymuth, 1989) have shown that thrust could be generated at low Reynolds numbers through pitching and plunging, Koochesfahani (1989) quantitatively connected the wake structure with the thrust generation.  Koochesfahani (1989) performed dye flow visualization and Laser Doppler Velocimetry (LDV) measurements of the wake structure behind a pitching airfoil and demonstrated the role of varying reduced frequency and pitching amplitude on the evolving wake structure and thrust force.  In these experiments, a NACA 0012 airfoil pitched about the quarter-chord with two amplitudes (2¼ and 4¼).  With the qualitative and quantitative information of the wake, the thrust coefficient, CT, was estimated based on the transverse profile of the time-averaged streamwise component of the velocity measured one chord length downstream of the trailing edge.  This allowed for a direct correlation between the  6 vortical structure and the thrust force.  At low reduced frequencies (k = f $ c/U$, where f is the pitching frequency) and amplitude, the wake structure was drag producing (see Figure 1.5).  As the reduced frequency increased, the drag decreased until it reached a critical reduced frequency at which the rows of clockwise and counter-clockwise vortices became aligned on the wake centerline, leading to a uniform average velocity profile.  This critical reduced frequency is called the crossover reduced frequency, kcr.  Further increasing the reduced frequency generated the jet-producing reverse von K⁄rm⁄n vortex street, leading to thrust.    Figure 1.5.  Thrust coefficient dependency on reduced frequency and pitching amplitude (Koochesfahani, 1989).  Flow visualization of vortical patterns from current computations are included as reference.  Figure courtesy of AIAA. By increasing the pitching amplitude, Koochesfahani (1989) observed an increasing number of vortices shed per half-cycle at lower reduced frequencies.  At !o = 4¼, two same-sign vortices could be shed during each half-cycle (referred to as a bifurcating wake), producing a double-wake time-averaged velocity profile.  By increasing the amplitude further, three vortices of the same sign could be shed during each half-cycle.  Multiple-same sign vortical patterns, however, were not generated at lower amplitudes.  Similar complex structures were also  7 observed by Schnipper, et. al., (2009) for pitching airfoils and Jones and Platzer (1998) and Lai and Platzer (1999) for plunging airfoils.  Increasing the pitching amplitude produced higher thrust and, thus, a lower crossover reduced frequency.  It was shown by Triantafyllou, et al., (1993) that the thrust becomes amplitude independent if scaled by a Strouhal number based on trailing edge peak-to-peak amplitude (St = fATE/U$; where ATE is the trailing edge peak-to-peak amplitude) as opposed to the reduced frequency.   However, simulations by Streitlien and Triantafyllou (1998) demonstrated that the time-averaged thrust would be over predicted if only the time-averaged streamwise velocity profile in the wake of an oscillating airfoil is used to estimate the thrust, as done by Koochesfahani (1989).  This overestimation was confirmed in computations by Liu and Kawachi (1999) and Ramamurti and Sandberg (2001) and experiments by Bohl and Koochesfahani (2009).  By including terms involving velocity fluctuations in the momentum integral estimation previously neglected, Bohl and Koochesfahani (2009) quantitatively showed that the thrust was lower, and the crossover reduced frequency was higher, than that obtained using only the time-averaged streamwise velocity.  This work also showed that the switch in vortex arrangement (von K⁄rm⁄n to reverse von K⁄rm⁄n) occurs before the switch from drag-to-thrust.  Thus, two crossover reduced frequencies are defined: one based on the switch in vortex arrangement and the other based on the drag-to-thrust switch.  Godoy-Diana, et al. (2008) also observed that the vortex arrangement switch occurs before the switch in streamwise force using a different airfoil and Reynolds number.  The reason the vortex arrangement switches before the drag-to-thrust switch is due to the additional momentum flux from the streamwise velocity fluctuations required to overcome the pressure reduction in the wake of the airfoil produced by the vortical structures (Bohl and Koochesfahani, 2009).    8 In order to quantitatively relate the behavior of the thrust to the vortical structure in the wake, Bohl (2002) and Bohl and Koochesfahani (2009) characterized the effect of reduced frequency on a variety of properties related to the strength, size, speed, and the pattern of the vortices in the wake.  The vortex strength was characterized by the peak vorticity (#z,peak) and circulation ((), the speed by the vortex convection speed (UC), the size by the core radius (rc), and the pattern by the streamwise and transverse spacing (a and b, respectively) between the vortices (see Figure 1.6 for a schematic).  Measurements were done at a location 0.5c downstream of the trailing edge.  The combination of these parameters would then impact the resulting thrust force on the airfoil, as demonstrated by Monnier, et al., (2015).  It was found that as the reduced frequency increased, the peak vorticity, circulation, and convection speed increased.  Conversely, the streamwise spacing decreased while the transverse spacing increased until it reached an observable asymptotic value.  The core radius, however, did not vary significantly with increasing reduced frequency at the measurement location.    Figure 1.6.  Schematic of vortex streamwise (a) and transverse (b) spacing.  9 1.2.! Reynolds Number Dependency A comparison between the thrust coefficient for a NACA 0012 purely pitching about its quarter-chord with 2¡ amplitude from experiments (Bohl and Koochesfahani, 2009, unpublished; Naguib, et al., 2011; Mackowski and Williamson, 2015) and viscous computations (Ramamurti and Sandberg, 2001; Young and Lai, 2004; Xiao and Liao, 2009; Jee and Moser, 2012; Zhang, et al., 2012) with inviscid flat plate theory (Garrick, 1957) and inviscid computations for a NACA 0012 (Ramamurti and Sandberg, 2001) is provided in Figure 1.7.  The Reynolds number for the experiments and viscous computations is approximately 12,000.  In general, the low Reynolds number data show lower CT values than the inviscid predictions.  This produces a drag-to-thrust crossover reduced frequency, kcr, that is generally higher than the inviscid predictions of approximately 1 for the flat plate (Garrick, 1957) and 3 from inviscid flow around the NACA 0012 (Ramamurti and Sandberg, 2001).  Koochesfahani (1989) commented that the discrepancy between the higher crossover reduced frequency in the viscous flow compared to the inviscid flow Òmay be expected because here is there substantial viscous drag that must be overcome that does not of course exist in the inviscid caseÓ.  The higher kcr of the low Reynolds number data and remarks by Koochesfahani (1989) are suggestive that there may be a Reynolds number dependence on the crossover reduced frequency.    10  Figure 1.7.  Thrust coefficient vs. reduced frequency (left), crossover reduced frequency vs. Reynolds number for NACA 0012 pitching about its quarter-chord with a 2¡ amplitude at a Reynolds number of approximately 12,000 (right). Another significant point from Figure 1.7 is that it is troubling to see the great deal of scatter amongst the computational data for this ÒcanonicalÓ case commonly used for low Reynolds number unsteady aerodynamic validation, where the crossover reduced frequency varies from approximately 4 to 10.  For the specific instances where the thrust curves from the literature were close to the inviscid predictions, the authors noted only that there was agreement in the trends when making comparisons with other data without making comment on the quantitative agreement.  However, there is good quantitative agreement between the experiments by Naguib, et al., (2011) and Bohl and Koochesfahani (2009; unpublished) with compressible flow computations by Young and Lai (2004) and Yu, et al., (2010).   Liu and Kawachi (1999) were the first to quantitatively observe a Reynolds number dependency on thrust, where the thrust increased as the Reynolds number increased.  This work used an incompressible Navier-Stokes solver.  However, only two Reynolds numbers within a very narrow range were used in this study: 7,200 and 12,000.  Other authors have also remarked  11 on a Reynolds number dependency on the thrust.  Dong, et al., (2006) performed computations to study the effect of aspect ratio on the wake structure and propulsive characteristics for ellipsoid wings undergoing combined pitch and plunge at insect-range Reynolds numbers (Re = 102).  The authors noted that at low Reynolds numbers, airfoils produce thrust only at high Strouhal numbers due to the high profile drag that must be overcome at the low Reynolds numbers.  Experiments on a low aspect ratio pitching plate by Buchholz and Smits (2008) over a range of Reynolds numbers between 3,000 and 15,000 showed that the thrust generally increased with increasing Reynolds number.  However, the size of the error bars (especially at high Strouhal number) made this trend difficult to ascertain.    High-order computations were also undertaken by Visbal (2009) to study transitional flow around an SD7003 undergoing pure plunging with a non-zero mean angle of attack.  For a specific reduced frequency of 3.93, an order of magnitude reduction in Reynolds number saw a switch in sign of the average thrust coefficient.  Ashraf, et. al., (2011) studied airfoil thickness effects at different Reynolds numbers (Re = 102-106) for an airfoil undergoing plunging and combined pitching and plunging.  The plunging amplitudes were 0.25 and 0.50 while the pitching amplitudes were 15¼-30¼.  The authors found that at low Reynolds number (Re = 102), the airfoils that were 12% or less produced thrust.  Airfoils thicker than this produced drag, which increased as the thickness increased.  As the Reynolds number increased to 2,000+, there was a peak in thrust and efficiency when the thickness was in the range 20%-30%.  Yu, et al., (2013) also investigated the effect of thickness on purely pitching, purely plunging, and combined pitch/plunge on two-dimensional airfoils at Re = 103, also using a high-order spectral difference method.  Yu, et al., (2013) noted that as the thickness-based Reynolds number increases, the wake structure gradually changes from drag-indicative to thrust-indicative.    12  Directly tied to thrust is the propulsive efficiency.  The propulsive efficiency, -., is the ratio between work caused by the propulsive force versus the input power (Garrick, 1936; Triantafyllou, et al., 1993).  Written in terms of standard non-dimensional quantities, the propulsive efficiency is defined as  -.'/0/1. (1.1) where the time-averaged /0 is   /0'23/04546706 (1.2) the time-averaged input power coefficient, /1., is   /1.'23/849:/;4+546706 (1.3) In the above expression, /84 is the time-varying lift coefficient, /;4 is the time-varying pitching moment coefficient and T is the oscillation period.  The plunging velocity, 9, and pitching velocity, +, are normalized by U$ and U$/c, respectively.   Triantafyllou, et. al., (1993) theorized that there exists an optimal Strouhal number range where the peak propulsive efficiency occurs.  This optimal Strouhal number range occurs between St = 0.25-0.35, even though parameters such as the trailing edge amplitude, pitch-axis location will affect the maximum efficiency.  By comparing the operating range of various fish and cetaceans (Re = 104-106), Triantafyllou, et. al., (1993) noted that the optimal Strouhal number range was not significantly dependent on Reynolds number.  Computational results by  13 Liu and Kawachi (1999) on a pitching airfoil at Re = 7,200 and 12,000 also showed a propulsive efficiency within this range (see Figure 1.8).  It is possible that since a small Reynolds number range was used by Liu and Kawachi (1999), they could not fully ascertain the dependency of Reynolds number on the range of optimal Strouhal number.    Figure 1.8.  Propulsive efficiency versus Strouhal number (Liu and Kawachi, 1999).  Figure courtesy of Elsevier. With the limited work found in the literature on the Reynolds number effect on the propulsive characteristics, the effect of Reynolds number on the vortex arrangement and its accompanying configuration switch has received even less attention in the literature.  Bohl and Koochesfahani (2009) characterized the effect of reduced frequency on the switch in orientation from a traditional von K⁄rm⁄n street to reverse von K⁄rm⁄n street for a sinusoidally pitching NACA 0012 at Re = 12,600 and 2¡ amplitude.  The authors observed that a visually aligned vortex array was produced at a reduced frequency of approximately 5.7 (St = 0.095).  Godoy-Diana, et al., (2008) investigated wake transitions for the tear-drop shaped airfoil at chord Reynolds number of approximately 1,200 and found a switch in vortex configuration at St > 0.1  14 for AD < 1, where AD is the trailing edge amplitude normalized by the airfoil thickness (see Figure 1.9).  Monnier, et al., (2015) examined the effect of flexibility on the wake properties of an airfoil with a flexible trailing edge.  The airfoil had a ÒheadÓ with a NACA 0036 profile and a ÒtailÓ with uniform thickness and a sharp trailing edge.  The chord Reynolds number was 2,100.  It was observed that the switch in transverse spacing, b, occurred at St ! 0.135-0.20 for the rigid airfoil, which was higher than the conditions for the pattern switch observed by Bohl and Koochesfahani (2009).  It is unclear whether or not this higher crossover Strouhal number is a result of shape effects or Reynolds number.  Figure 1.9.  a) Schematic of airfoil and kinematics considered and b) oscillation amplitude vs. Strouhal number.  The chord Reynolds number is approximately 1,200.  The blue line in the right image represents the switch from traditional von K⁄rm⁄n street to reverse von K⁄rm⁄n vortex patterns.  Figure courtesy of American Physical Society. 1.3.!Kinematic Trajectory Effects  The vast majority of investigations found in the literature have focused on airfoils oscillating sinusoidally while less attention has been given to motions that are non-sinusoidal in nature.  The first work involving non-sinusoidal oscillations was Koochesfahani (1989).   In this work, the pitching trajectory was made asymmetric, where the degree of asymmetry was denoted  15 by a symmetry parameter, S, defined as the duration of the pitch-up (or downstroke for a plunging airfoil) divided by the total oscillation period.  This is shown schematically in Figure 1.10.  Based on this nomenclature, a sinusoidal trajectory becomes a special case of the asymmetric trajectory when S = 50%.  Figure 1.10.  Schematic of asymmetric trajectory definition (Koochesfahani, 1989).  Figure courtesy of AIAA. Koochesfahani (1989) found that asymmetric pitching had a very pronounced influence on the evolving vortical structure (see Figure 1.11).  For faster pitch-up (S = 38%), a single counter-clockwise vortex and two clockwise vortices were shed during each pitching cycle.  A slower pitch-up and a faster pitch-down produced a mirrored structure, where two counter-clockwise vortices shed during the slow pitch-up and a single clockwise vortex shed during the fast pitch-down.  It was hypothesized that the force history would be affected due to the breaking of symmetry, though this effect was not characterized.  By increasing the pitching amplitude, more complex vortical structures were observed.  Interestingly, Naguib, et. al., (2011) reproduced the very intricate mean and r.m.s. profiles of Koochesfahani (1989) using a simple finite-core vortex array model whose only inputs were the relative position, core radius, and circulation of the vortices.    16  Figure 1.11.  Flow visualization for three asymmetries: S = 38%, 50%, and 61% (Koochesfahani, 1989).  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  Figure courtesy of AIAA. More recently, Sarkar and Venkatraman (2006) studied an asymmetric plunging trajectory using an incompressible, viscous discrete vortex method.  Flow visualization of the discrete vortex particles showed complex flow patterns in the wake, reminiscent of those seen in Koochesfahani (1989).  This work further illustrated that the fast part of the motion, a single positive-sign vortex sheds while multiple vortices shed during the slower part.  Sarkar and Venkatraman (2006) also were the first to directly show the impact of asymmetry on the forces and efficiency.  They found that the airfoil experienced an increase in thrust as the plunging motion became more asymmetric.  The increase in thrust was also symmetric about S = 50%, such that S = 38% would produce the same thrust as S = 62%, for example.  In general, an asymmetric trajectory produced an increase in propulsive efficiency.  The authors also commented that the time-averaged lift was non-zero when using an asymmetric motion, and increased with asymmetry and reduced frequency, though no quantitative data was presented.  Asymmetric flapping using oscillatory rolling kinematics was studied by Devranjan, et. al., (2013) using an inviscid discrete vortex method and experiments.  The discrete vortex  17 method incorporated flexibility using a cantilever simplification.  This work showed that asymmetric flapping can be an effective means of attaining lift, with an optimal flexibility existing for a given amplitude.  Yu and Tong (2005) also studied a complicated planar motion representative of insect flight.  In this work, stroke asymmetry was studied using an incompressible, potential flow solver.  Yu and Tong (2005) showed that asymmetry enhances time-averaged thrust while there was little impact on the time-averaged lift.  Asymmetric pitching was also investigated by Lu, et al., (2013) in the context of higher Reynolds number (Re = O(105)) and low reduced frequencies (k = O(10-1)).  They found that the force peak values, force time-histories, and hysteresis loop widths are noticeably affected by asymmetry.  Asymmetry also delayed both the formation and development of the leading edge vortex. 1.4.! Present Research It is clear from the literature that a dependency of the vortical flow on the Reynolds number exists.  However, it is even more apparent that an extensive, quantitative analysis of the effect of Reynolds number on both the forces, efficiency, and wake structure is lacking.  One of the main purposes of the present research is to characterize the effect of Reynolds number on the forces, efficiency, and wake structure for a pitching airfoil.  It is also of interest to further explore the role of asymmetry on the forces and flow structure of a pitching airfoil.  The underlying interest in asymmetric airfoil motion trajectory comes from the fact that natural flyers utilize asymmetric wing beats (Park, et. al., 2001; Ros”n, et al., 2004).  Thus, it is beneficial to study the underlying physics of the flow created by asymmetric kinematics which natural flyers use and MAVs may exploit.   In this work, high-fidelity, computations are undertaken to study both the effect of Reynolds number and trajectory asymmetry on the forces and flow structure of an oscillation  18 airfoil at low Reynolds number (Re = 2,000-22,000).  A schematic of the considered flow is given in Figure 1.12.  The flow is assumed to be two-dimensional, limiting the flow to be laminar.  This assumption is valid for the low Reynolds numbers of interest and confirmed by experiments (Koochesfahani, 1989).  This assumption is also the first step in studying these flows given the large parameter space.  The oscillatory motion selected is pitching performed about the quarter-chord.  The choice of kinematics is to allow for direct qualitative and quantitative comparison with the experiments by Koochesfahani (1989), Bohl and Koochesfahani (2009), and Naguib, et. al., (2011).  Two coordinate systems are used.  A global coordinate system is denoted by X and Y and originates at the leading edge.  A wake coordinate system is also used and is denoted by x and y, which originates at the trailing edge.  All forces were obtained by integrating the pressure and shear stress distributions around the surface of the airfoil.  The lift, L, is positive in the positive Y direction while drag, D, is positive in the positive X direction.  Note that positive thrust is in the negative X direction.  For this work, the high-order, extensively-validated, FDL3DI Navier-Stokes solver is used (Gaitonde and Visbal, 1998; Visbal and Gaitonde, 1999).   Figure 1.12.  Schematic of investigated flow with important parameters labeled. What sets this work apart from other high-order computations for oscillating airfoil studies (Visbal, 2009; Liang, et al., 2011; Ou, et al., 2011, Yu, et al., 2011; Yu, et al., 2013) are as follows.  The first, and most distinctive, difference is that this is the first attempt, to the  19 authorÕs knowledge, to use high-order computational tools to obtain quantitative data of the wake flow structure (velocity profiles, vortex properties, etc.).  This is made possible by implementing an overset grid approach.  The second distinctive difference is that this work utilizes freestream Mach numbers lower than those typically found in the literature employing a compressible flow code.  In this work, freestream Mach numbers as low as 0.005 are used.  The appropriateness of the freestream Mach number selections for simulating incompressible flow is examined systematically.  The remainder of this document is divided as follows.  Chapter 2 discusses computational considerations, such as the asymmetric motion trajectory definition, the numerical method, validation, Mach number effects, and grid remapping/vortex tracking.  Chapters 3 and 4 present the effect of Reynolds number on the forces and vortex properties of a sinusoidally pitching airfoil.  Chapters 5 and 6 demonstrate the role asymmetry has on the forces and flow structure for the pitching airfoil.  Chapter 7 summarizes the outcomes of the present research and provided recommendations for future research.     20 CHAPTER 2.  COMPUTATIONAL CONSIDERATIONS Before the results of this work can be presented, several computational considerations must be addressed.  This chapter is organized as follows: Section 2.1 provides a discussion of the geometry and the motion trajectory of the pitching airfoil.  Sections 2.2 and 2.3 detail the computational method and domain, respectively.  Section 2.4 contains a description of the initial and boundary conditions. Section 2.5 is focused on issues related to the solution convergence and freestream Mach number.  Sections 2.6 presents validation for the symmetric and asymmetric pitching cases.  Finally, Section 2.7 provides additional details related to time-periodic convergence and data remapping. 2.1.  Geometry and Trajectory Definition The geometry considered in this work was the NACA 0012 airfoil, which is 12% thick and symmetric about the centerline.  The airfoilÕs trailing edge was rounded with a radius of 2.2x10-3c to facilitate the use of an O-grid topology (see Figure 2.1).  This radius is approximately 3% of the boundary layer thickness at the trailing edge for the static airfoil at the highest Reynolds number of 22,000 studied here, which had the lowest boundary layer thickness of the investigated Reynolds numbers.    Figure 2.1.  NACA 0012 geometry with rounded trailing edge.  21 The generalized pitching motion used an asymmetric pitching formulation, in which sinusoidal pitching was a special case.  The asymmetric pitching motion is defined in equation (2.1) and accurately describes the trajectory created by the function generator waveform used in experiments by Koochesfahani (1989).  Part 1: 0 " t < TPU/2    +<4'+=>?@A4B3  Part 2: TPU /2 " t < T-TPU /2    +C4'D+=>?@A2DB43D2# (2.1) Part 3: T Ð TPU /2 " t < T    +E4'+=>?@AB43D2  where +<4 - +E4 are the angle of attack time-history definitions for parts 1-3, T is the oscillation period, TPU = TS is the pitch-up period, and S is the symmetry parameter defined by the pitch-up period divided by the entire oscillation period (Koochesfahani, 1989).  Based on this definition of S, a sinusoidal trajectory occurs when S = 50%, while a faster pitch-up and slower pitch-down occurs at S < 50% and vice-versa for S > 50%, as illustrated in Figure 2.2, where f is the oscillation frequency.   22  Figure 2.2.  Angle of attack time-histories for S = 38%, 50%, and 62%, normalized by the pitching amplitude, !o. Since the trajectory described by equation (2.1) is piecewise continuous, there is a physically unrealistic discontinuity in the acceleration at t = TPU/2 and t = T Ð TPU/2.  In order to temporally resolve this acceleration jump, the angle of attack was smoothed using a Gaussian filter.  This resulted in a smoothed trajectory +4, defined by equation (2.2), where the size of the filter window, ", was based on a desired percentage of the pitch-up period, TPU.    +4'F4D4+454GHG'2#AICJHK6H6LCML+454GHG (2.2) For completeness, an analytical derivation of the smoothed trajectory using a Fourier series expansion is given in Appendix A.  The impact of the smoothing on the acceleration, and resulting time-varying lift, for S = 38% can be seen in Figure 2.3 for the case of reduced frequency of 6.68 and pitching amplitude of 2¡.  A small filter size of 0.5%TPU follows the unsmoothed trajectory (0.0%TPU) closely.  For the example shown in the figure, there are approximately 240 points capturing the sharp acceleration jump for 0.5%TPU. As the filter size increased, there is an increasingly dampened response in the forces around the phase where the angle of attack changes direction (see Figure 2.3).  This is also shown in Figure A.3 in Appendix.   23 The effect on the average forces is a slight decrease in thrust and a slight alteration in shape of CL vs. S curves (see Figure A.4 in Appendix A).  Trajectory smooth also slightly affected the wake structure, as shown in Figure A.5 in Appendix A.  The smoothing did not alter the sinusoidal trajectories.  For each asymmetry, the 0.5%TPU data collapsed on the unsmoothed (0.0%TPU) trajectory data and thus is used for all asymmetric motion trajectory computations in this work.  Only values of S " 50% were considered, with the assumption that S > 50% would produce mirrored results.    Figure 2.3.  Pitching acceleration and lift time-history for S = 38%, illustrating influence of ". In order to ensure a smooth startup of the airfoil motion, the filtered asymmetric trajectory, +4, was multiplied by a ramping function (Visbal, 2009), given by   +4'+42DJHN)O6P6Q (2.3) In equation (2.3), to is the time-scale for the motion to achieve 99% of its periodic state.  The value of to in this work is selected so that the airfoil reached 99% of its periodic state after one pitching cycle.  Remaining artifacts from the initial transients are minimized by advancing simulations twenty or more pitching cycles.    24 2.2.  Computational Method A brief summary is now given of the numerical method method based on two-dimensional flow, with a more complete description based on three-dimensional flow given in Appendix B.  This work used the high-order, extensively validated FDL3DI Navier-Stokes solver, developed at WPAFB (Gaitonde and Visbal, 1998; Visbal and Gaitonde, 1999).   2.2.1.  Governing Equations  The governing equations being numerically solved are the unsteady, compressible, two-dimensional Navier-Stokes equations.  Using an appropriate transformation from Cartesian coordinates (x, y, t) to curvilinear coordinates (), *, %), the Navier-Stokes equations are solved in strong conservative form, given by  RRSTU:RV1RW:RX1R-'2YJRVZRW:RXZR- (2.4) The solution vector of conservative variables, T, is defined by  T'[\[]\[^\[_0 (2.5) In the above expression, ' is the density, u and v are the Cartesian velocity components, and E is the specific total energy defined by  _'3``D2aGC:2#]C:^C (2.6) where, T is the temperature, + is the ratio of specific heats, and M$ is the freestream Mach number (M$ = U$/a$, where U$ and a$ are the freestream velocity and speed of sound, respectively).  The transformation Jacobian between Cartesian and curvilinear coordinate is U, where U'RW\-\SPRb\c\4.    25 The inviscid fluxes, V1 and X1 are defined by   V1'2U[d[]d:Wef[^d:Wgf[_:fdDW6fhX1'2U[i[]i:-ef[^i:-gf[_:fiD-6f (2.7) The subscripts x, y, and z in equation (2.7) denote a partial derivative.  The directional contravariant velocities, U and V, in equation (2.7) are  d'W6:We]:Wg^i'-6:-e]:-g^ (2.8) The viscous fluxes, VZ and XZ, in are given by  VZ'2U(WeSee:WgSgeWeSeg:WgSggWe]See:^SegDje:Wg]Sge:^SggDjgXZ'2U(-eSee:-gSge-eSeg:-gSgg-e]See:^SegDje:-g]Sge:^SggDjg (2.9) The shear stress components See, Seg, and Sgg are defined by equation (2.10) upon implementing StokeÕs hypothesis for bulk viscosity (i.e., , = -2/3µ).  See'#kl#We]m:-e]nDWg^m:-g^nSeg'Sge'lWg]m:-g]n:We^m:-e^nSgg'#kl#We^m:-e^nDWg]m:-g]n (2.10)   26  The heat flux, jo and jp, is  jo'D2`D2aGClqrWe3m:-e3njg'D2`D2aGClqrWg3m:-g3n (2.11) where, µ is the dynamic viscosity while Pr is a constant Prandtl number (Pr = 0.72 for air).   In order to close the set of equations, the Ideal Gas Law (equation (2.12)) and SutherlandÕs law for viscosity (equation (2.13), where Sc = 0.38 for air) are used.  All flow variables are normalized by their freestream value, except pressure which is normalized by twice the dynamic pressure (i.e., '$U$2).  The length scale is the chord length, c, and the physical time has been normalized by U$/c.  f'[3`aGC (2.12)  l'3EPC2:Bs3:Bs (2.13) As stated before, these equations correspond to the two-dimensional Navier-Stokes equations since three-dimensional computations are very prohibitive given the very large parameter space considered in this work.  However, the assumption of two-dimensional flow is supported by experiments (Koochesfahani, 1989).      27 2.2.2.  The Time Marching Scheme The time marching algorithm (equation (2.14)) is the linearized, approximately factored method of Beam and Warming (1978), which has been diagonalized (Pulliam and Chaussee, 1981) to improve computational efficiency.  The time-marching algorithm is supplemented by Newton-like subiterations to maintain second-order accuracy by reducing errors caused by the linearization, approximation factorization, and diagonalization, and explicit implementation of physical boundary conditions (Visbal, et al., 2003).  Although, three subiterations are typically used (Visbal and Gaitonde, 2002), seven subiterations are used in the current work.  In ÒdeltaÓ form, the time marching scheme is given by  2Ut7<:uS2:vwmCRV1tRTD2YJRVZtRTUt7<2Ut7<:uS2:vwnCRX1tRTD2YJRXZtRTuT'DuS2:v2Ut7<:2:vTtD2:#vTxDvTxH<uS:Tt2Uyt:wmOV1tD2YJVZt:wnOX1tD2YJXZt (2.14) In equation (2.14), uS is the time-step size, p denotes subiteration step, and uT'Tt7<DTt, where Tt7< corresponds to the p + 1 approximation of the n + 1 time level and p denotes the subiteration level within the current time-step.  In the first subiteration, p = 1 corresponds to Tx.  As p $ %, Tt$Tx7<.  The inviscid flux Jacobians are RV1PRT and RX1PRT while RVZPRT and RXZPRT are the viscous flux Jacobians.  First-order Euler-implicit and second-order three-point backward schemes result when - = 0.0 and - = 0.5, respectively.    28 The . operator on the left- and right-hand-side of the equation represents spatial derivatives computed using finite difference schemes in the direction of the subscript (i.e., ) and *) of the indicated order of accuracy, denoted by the number in the superscript.  For example, wmC denotes second-order central difference scheme in the ) direction while wmO is the sixth-order compact scheme in the same direction.  Since the solution is converged with subiterations, the lower order of accuracy on the implicit side does not impact the final solution (Garmann, 2013).    Not shown are second- and fourth-order, nonlinear artificial dissipation terms (Jameson, et al., 1981; Pulliam, 1986) that are appended to the implicit, left-hand-side for enhancing stability.  This results in a scalar, pentadiagonal system of equations that is more computationally expensive than solving a scalar tridiagonal system of equations but dramatically less expensive than solving the original block tridiagonal system.  The use of subiterations eliminates the impact of the dissipation terms on the final solution, and thus can be chosen specifically for stability (Visbal and Gaitonde, 2002).  The implementation of this time-marching method has been successfully applied to a variety of flows (e.g. Visbal and Rizzetta, 2002; Visbal and Gaitonde, 2002; Visbal, et al., 2003; Rizzetta, et al., 2008).  The second-order time marching scheme is used for this research.      29 2.2.3.  Spatial Discretization and Filter High-order spatial derivatives are computed with the compact finite difference scheme, described by Lele (1992) to obtain spectral-like resolution.  For the interior grid points, the finite difference approximation of variable / at point i can be expressed as a linear combination of the functional values of / at point i and the surrounding points.  For a five-point stencil, this combination is represented by  z{|H<}:{|}:z{|7<}'~{|7<D{|H<#:{|7CD{|HC* (2.15) The solution of the equation (2.16) requires the inversion of a tridiagonal matrix and results in both explicit and compact schemes of varying orders of accuracy with appropriate choice of the coefficients a, b, and z.  If z = 0, the derivative at point i is decoupled from its neighbor points and the traditional second-order and fourth-order central difference schemes are recovered based on a and b.  If z & 0, the derivative at point i is coupled to its neighbor points and results in fourth-order and sixth-order compact difference schemes based on z, a, and b. Although eighth- and tenth-order accuracy can be achieved using a seven-point stencil, it is much more computationally expensive due to the coupling of additional points.  Garmann (2010) commented that the sixth-order scheme is an effective balance between computational accuracy and efficiency.  For this work, the sixth-order accurate compact scheme was used for computing spatial derivatives at the interior points, with a = 14/9, b = 1/9, and z = 1/3.  In order to maintain the tridiagonal nature of the compact scheme, special treatment is required at the boundary points 1, 2, N, and N-1.   Since the high-order compact finite difference scheme is based on a centered stencil and thus non-dissipative, it is susceptible to numerical instabilities due to the growth of unstable  30 frequency modes that originate from mesh non-uniformity, boundary condition approximations, and non-linear flow features (Gaitonde and Visbal, 2000).  In order to maintain solution stability, a highly discriminating filter dampens the poorly evolved high-frequency content of the solution not properly resolved by the compact difference scheme (Gaitonde and Visbal, 1999; Gaitonde and Visbal, 2000).   The high-order, low-pass filter (Gaitonde, et al., 1999) for interior points is derived in a similar way as the high-order compact scheme and is based on templates presented by Lele (1992) and Alpert (1981).  The filtered values { of the variable { are obtained by solving the following tridiagonal system of equations   z•{|H<:{|:z•{|7<'~†#{|7†D{|H†‡†…— (2.16) Up to a tenth-order accurate filter can be obtained with an eleven-point stencil (N = 5).  As the boundary is approached, special treatment is needed due to the large centered stencil.  Two approaches have been used in the literature.  One approach reduces the filter order while maintaining a centered stencil while a one-sided formula is used at point 2 (Visbal and Gaitonde, 1999).  The second approach uses high-order one-sided formulas (Gaitonde and Visbal 2000).  In this work, the first approach is used.   The filter is applied to the conservative variables during each time-step or subiteration.  Since the filter is one-dimensional, it is applied sequentially one direction at a time.  Before filtering a subsequent direction, the solution is updated with the filtered values.  To minimize potential bias, this sequence is alternated between various permutations (Visbal and Gaitonde,  31 1999).  The filtering sequence alternates once per time-step; thus the same sequence is applied after each subiteration of a specific time-step.   For the current work, the filter was defined as follows.  The interior point used an eighth-order filter in conjunction with z– = 0.4.  The order of the centered filter was reduced from eighth-order at point 5 (and N-4) to sixth-order at point 4 (and N-3) and fourth-order at point 3 (and N-2)  while maintaining z– = 0.4.  Point 2 and N-1 used a fourth-order one-sided filter with z– = 0.45.  Similar filtering schemes were adopted by Galbraith (2009), Garmann (2010), and Garmann (2013).   2.3.  Computational Domain One of the goals of the present work is to study the wake evolution far downstream from the airfoilÕs trailing edge.  In computational aerodynamic research, it is common to use an O-grid to discretize the fluid domain due to its simplicity in implementation, lower computational cost, and effectiveness in accurately capturing the forces and near wake region.  However, stretching of the grid causes the vortices to rapidly diffuse (see Figure 2.4).  Resolving the wake with the O-grid requires the clustering of grid points, which creates the disadvantage of having high resolution in locations where it is not necessary.  Figure 2.4.  Instantaneous spanwise vorticity fields, #zc/U$, using only O-grid.  The Reynolds number is 12,000, reduced frequency is 6.68, and 2¡ pitching amplitude.  32 An alternative is to use an overset grid approach, which allows for grid resolution flexibility while maintaining a structured grid capable of using the high-order compact difference and filter schemes.  Connectivity between the grids was established during pre-processing using the NASA-developed software Pegasus 5 (Suhs, et al., 2002).  This was done by first identifying donor points and receiver points, then calculating their interpolation offset, w, in computational space.  For reference, donor points are points in which data is interpolated from and receiver points are points in which data is interpolated to.  The low-order interpolation coefficients are then extended to high-order using AFRL-developed software BELLERO (Sherer, et al., 2006).  It computes the high-order interpolation coefficient and offsets with the following procedure (Sherer and Scott 2005; Sherer and Visbal, 2007).  First, the donor point stencil obtained from Pegasus is expanded such that the number of points in each coordinate direction equals the desired interpolation order of accuracy (N).  The receiver point is made as centered as possible within the donor point stencil during the expansion process.  Next, the high-order offsets in each coordinate direction are obtained by solving a set of isoparametric mapping equations based on the expanded stencil.  The final step is the calculation of the interpolation coefficients.  The interpolated flow variable at the receiver point r is calculated from the donor point d using the following formula   {ƒ'q⁄nq|m‡H<|…—K{1‹7|\›‹7⁄‡H<⁄…— (2.17) where {ƒ is the interpolated scalar at the receiver point, {− is the value from the donor point, and the indices Id and Jd represent the location of the base donor point in computational space (Sherer and Visbal 2007).  The interpolation coefficients q⁄n, and q|m are calculated using  33  qx‰'D2‡7xH<„D“D2”“”w‰D‘‡H<|…—|’⁄ (2.18) where n = i, j and ‚ = ™\Kfi)KK Sixth-order interpolation was used in this work. The overset computational domain in this work consisted of a body-fitting O-grid around the airfoil and four Cartesian grids around the O-grid (see Figure 2.5).  Although the forces were not affected using the overset grid, its implementation allowed for the wake to be captured (see Figure 2.6).  As verification of the diffusion rate of the vortices, the decay of peak vorticity was compared to the theoretical decay of a two-dimensional vortex with a Gaussian vorticity distribution.  There was excellent agreement between the computational result and theory, which is shown in Appendix D.  Figure 2.5.  Overset computational domain shown at three different scales.  34  Figure 2.6.  Comparison between the instantaneous spanwise vorticity fields, #zc/U$, using O-grid (top) and overset grid (bottom).  The Reynolds number is 12,000, reduced frequency is 6.68, and 2¡ pitching amplitude. The O-grid around the airfoil was generated using the commercial software Pointwise (Gridgen User Manual, 1997), which uses a marching method to solve the system of two-dimensional partial differential equations, given by equation (2.19), where subscripts denote partial derivatives and iW\- is the cell volume distribution.   bmbn:cmcn'(bmcn:cmbn'iW\- (2.19) The O-grid was composed of 651 grid points in the tangential, ), direction and 180 grid points in the normal, *, direction.  The grid spacing in the )-direction at the leading edge and trailing edge was 5.0x10-4c and 2.5x10-4 c, respectively.  The initial grid spacing in the *-direction away from the airfoil surface was 5.0x10-5 c.  The Cartesian grid behind the airfoil had a uniform spacing of 2.5x10-3c over four chord lengths downstream of the trailing edge in the streamwise direction and approximately one and a half chord lengths in the transverse direction.  The cell size in the overlap regions was approximately uniform between all grids.  The four grids were then stretched rapidly to 140c away from the airfoil.  Two additional planes were added ±0.2c from  35 the center plane in the spanwise direction, 0, in order to simulate two-dimensional flow using the current three-dimensional Navier-Strokes solver.  The total number of grid points was 5 million, which was decomposed into 71 subdomains for parallel solving to improve speed.   2.4.  Initial Conditions and Boundary Conditions  The initial and boundary conditions are as follows.  Simulations were initialized with an instantaneous solution for the static airfoil at zero angle of attack.  At the airfoil surface, a no-slip, adiabatic condition was used in conjunction with a zero normal pressure gradient.  The Neumann boundary conditions at the surface were computed with fourth-order accuracy using equation (2.20) (Gaitonde and Visbal, 1998).  R{R-'*fl{CDkŁ{E:2Ł{NDk{Œ#Š (2.20) where subscripts 2-5 denote the second through fifth grid points away from the surface in the - direction.  The surface velocities of the airfoil (]Ÿ and ^Ÿ) based on the imposed pitching motion were  ]Ÿ'+4cDct^Ÿ'D+4bDbt (2.21) where xp and yp are the pitch axis coordinates and +4 is the pitching velocity.  For all computations, the airfoil pitched about its quarter-chord (xp = 0.25c and yp = 0.0c).  Freestream conditions were specified at 140c from the airfoil surface.  This rapid stretching, in conjunction with the low-pass, high-order filter, allows for damping of spurious pressure wave reflections off the boundary (Visbal and Gaitonde, 2001).  Spatial periodicity to enforce an azimuthally continuous solution was specified along the five-point grid overlap in the tangential direction, ).   36 A zero gradient condition computed with first-order accuracy was imposed on the additional two planes in the spanwise direction for simulating two-dimensional flow with the current, three-dimensional Navier-Stokes solver. 2.5.  Convergence Studies 2.5.1.  Convergence Study Summary  In order to verify the fidelity of the simulations, grid resolution, time-step, and subiteration convergence studies were conducted for both symmetric and asymmetric trajectories.  This section briefly documents the convergence studies, with a more complete discussion provided in Appendix C.  A summary of the final simulation parameters used for all computations in this work is given in Table 2.1.  Table 2.1.  Final simulation parameters. Parameter Symmetric Traj. Asymmetric Traj. Grid Fine (5x106 points) Fine (5x106 points) Time-step 5.0x10-5 2.5x10-5 Subiterations 7 7 Mach Number 0.015 0.005  2.5.2.  Grid, Time-step, Subiteration Studies  For presenting grid, time-step, and subiteration convergence, the asymmetric pitching case was selected.  The flow conditions for the asymmetric trajectory convergence studies were as follows: Re = 12,000, M$ = 0.005, k = 6.68, !o = 2¡, and S = 30%.  To verify grid independence, the computational domain was coarsened using a bi-cubic spline interpolation approach (Press, 1992; Garmann, 2010) to allow for a global coarsening of the grid while maintaining similar grid point distributions.  Two grids coarser grids were created, where one grid had 75% (Medium) and the other had 50% (Coarse) of the number of grid points in the baseline grid.  The leading edge spacing, trailing edge spacing, initial normal spacing away from  37 the airfoil surface, and uniform wake spacing for the two coarser grids are given in Table 2.2, with the Fine Grid spacing information included for reference.  The resulting Fine, Medium, and Coarse overset grids contained 5 million, 3.7 million, and 2.5 million grid points, respectively.  The time-step was varied between 5.0x10-5 to 1.25x10-5 and the number of subiterations, nsubit, was increased from 3 to 9.  To present convergence, time histories of lift and instantaneous vorticity profiles measured one chord length from the trailing edge are at a phase of 0.5T are shown in Figure 2.7. Table 2.2.  Grid spacing information for three grids of varying resolution. Grid !"LE !"TE Initial # Uniform Wake Spacing Fine 5.0x10-4c 2.5x10-4c 5.0x10-5c 2.5x10-3c Medium 5.8x10-4c 2.9x10-4c 5.7x10-5c 2.9x10-3c Coarse 7.1x10-4c 3.5x10-4c 7.0x10-5c 3.5x10-3c    The lift time histories show negligible quantitative differences between the coarsest and finest parameters for the three studies.  The vorticity profiles show greater sensitivity.  All three grids produce the same maximum value of vorticity at this pitching phase and location in the wake.  The shape of the vorticity profile from the medium and fine grids is identical while the coarse grid produces a shift in the profile in the positive transverse direction and a difference in shape in the transverse range 0.147 < y/c < 0.29 (highlighted in the figure).  Similar observations are also seen with respect to time-step size and number of subiterations.  The two smaller time-steps produce the same vorticity profile while the larger time-step shows produce some deviation in the transverse range 0.2 < y/c < 0.29 (highlighted in the figure).  Three subiterations produce a similar profile shape but lower maximum vorticity (highlighted in the figure).  When the number of subiterations, nsubit, increases to five or more, the profiles collapse.  At phases in which the vorticity profiles have low maximum values (#zc/U$ ! 1), some noise is observed in the profiles.   38 Increasing the number of subiterations reduces this low-level noise (see Figure C.2 and C.9 in the Appendix).  Seven subiterations was considered a good balance between accuracy and speed.  Figure 2.7.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at "/T = 0.50 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M$ = 0.005, k = 6.68, !o = 2¡, and S = 30%. 2.5.3.  Effect of Freestream Mach Number A challenging aspect of this work has been using a compressible flow solver to simulate incompressible flow, which requires a low freestream Mach number.  The reason for the challenge is that it is well known that the solution accuracy of compressible solvers deteriorates as the Mach number decreases towards zero (Volpe, 1993).  Although low Mach number pre-conditioning techniques exist that reduce these errors (Turkel, 1987), it requires that the time marching scheme be modified.  This, however, can be non-trivial for well-established flow solvers, such as FDL3DI.  39 In order to use existing compressible codes, Volpe (1993) showed that having high grid resolution can reduce the loss of accuracy in the context of inviscid, external flow.  Although not the intent of using the overset grid in this work, the low Mach number simulations do benefit from the higher grid resolution.  Young (2005) considered freestream Mach numbers in the range of 0.025-0.200 for plunging airfoils using a Navier-Stokes solver and showed that the force data visually converged towards the incompressible unsteady panel method results as Mach number decreased.  He also observed good visual agreement in the wake structure between incompressible flow visualization experiments and compressible simulations at M$ = 0.050.  A brief overview of the Mach number effect will now be discussed, with additional results in Appendix C. In order to determine the Mach number requirement for simulating the incompressible flow limit for the symmetric and asymmetric trajectories, force and wake vorticity profile convergence studies were performed on a variety of cases using freestream Mach numbers ranging from 0.005 to 0.100.  The Mach number effect on the peak-to-peak amplitude of lift, peak-to-peak amplitude of drag, and average thrust coefficients for sinusoidal pitching is shown in Figure 2.8 for three reduced frequencies (k = 4, 8, and 12) at a Reynolds number of 12,000 and pitching amplitude of 2¡.  At k = 4 and 8, the peak-to-peak amplitude of lift and drag decrease monotonically as M$ decreases monotonically.  At k = 12, the lift and drag amplitudes initially increases as M$ decreases from 0.075 to 0.100, then decreases monotonically as M$ decreases to zero.  When M$ " 0.015, the error between the lift and drag peak values is less than 3%.  It is also observed that the phase of maximum lift and drag converge as Mach number decreases as well, which is visually shown the lift and drag histories in Appendix C.  The thrust coefficient increases as M$ decreases from 0.075 to 0.100 for each reduced frequency (though this is  40 difficult to ascertain for k = 4 due to the scale of the figure), then decreases monotonically as M$ decreases to zero.  The maximum difference in the average thrust between successive Mach numbers is on the third decimal point and only occurred at the highest reduced frequencies considered in this work.  Similar findings are observed in the convergence of peak-to-peak amplitude and average value of lift and drag for the asymmetric case.  Figure 2.8.  Mach number effects on the peak-to-peak amplitude of lift and drag coefficients, and average thrust coefficient for Reynolds number of 12,000 and 2¡ pitching amplitude.  The asymmetric trajectory displayed greater significant sensitivity to the freestream Mach number.  Figure 2.9 shows the effect of Mach number on the time-varying lift for k = 6.68 and S = 38%.  At "/T = 0.19 and 0.81, the airfoil changes directions and experiences a sharp acceleration.   Shortly after the high accelerations, a ÒbumpÓ develops lift time history that becomes more pronounced as M$ decreases.  Inspection of the pressure field (not shown) indicates that a pressure pulse radiates from the airfoil after the acceleration jumps and causes the dynamic response observed in the lift.  The response time decreases as the Mach number decreases due to the traveling speed of the pressure waves (i.e., speed of sound) increasing.  In order to capture the incompressible limit considering the dynamics response in the forces present  41 in the asymmetric cases, the freestream Mach number is selected to be 0.005 for the asymmetric pitching investigation.    Figure 2.9.  Lift history for S = 38% illustrating effect of decreasing freestream Mach number for Re = 12,000, k = 6.68, and !o = 2¡. In addition to the influence of Mach number on the forces when the trajectory was asymmetric, there was a strong dependence of the vortical structure in the wake on Mach number.   Figure 2.10 shows instantaneous spanwise vorticity fields for three Mach numbers (0.100, 0.050, and 0.005) and dye visualization from the experiment by Koochesfahani (1989).  The asymmetry is 38%.  As shown in the figure, at each Mach number one positive vortex and two negative vortices shed per pitching cycle.  However, the negative vortices pair downstream after only one chord length for M$ = 0.100.  It is not until at least M$ = 0.050 that the computations produced a flow pattern similar to the experiment, where the negative vortices stay separate instead of pair.  Further reduction in Mach number to 0.005 produces slight differences in the downstream evolution of the wake structure compared to M$ = 0.050, primarily in the strength of the vorticity which formed into the braid and downstream orientation of the negative  42 vortices.  Thus, selection of M$ = 0.005 as for simulating the incompressible limit for the asymmetric pitching study is warranted.    Figure 2.10.  Spanwise vorticity fields, #zc/U$, S = 38% at %/T = 0.0, illustrating influence of M$.  Dye-visualization from the experiment by Koochesfahani (1989) is included for reference at the bottom.  Bottom figure courtesy of AIAA. One final Mach number comparison of interest is the Mach number effect on the density variation.  As the Mach number decreases towards zero, the change in density from its freestream value, 1', also must decrease towards zero.  To show this, the maximum and minimum density change experienced in the flow, 1'max and 1'min, are plotted versus freestream Mach number for a reduced frequency of approximately 12 and 2¡ amplitude in Figure 2.11.  As shown in the figure, the maximum and minimum density variation decreases with decreasing freestream Mach number.  Once the freestream Mach number reached approximately 0.015, the difference between the maximum and minimum density change is less than 0.1%.  It is  43 particularly interesting that the density variation dependence on freestream Mach number can be well represented by the relationship    Ž['Ž[=aGCaG\=C (2.23) In equation (2.23), 1'o and M$,o denote the reference conditions.  In the figure, the reference freestream Mach number is 0.100.   Figure 2.11.  Density variation towards the incompressible limit with decreasing freestream Mach number. 2.6.  Validation Although the flow solver has been successfully validated for a variety of flow conditions (Gaitonde and Visbal, 1999; Visbal and Gaitonde, 2002; Sherer and Visbal, 2007), it was necessary to perform independent qualitative and quantitative validation of the solver at the low Reynolds number and Mach number conditions of interest using the simulation parameters listed in Table 2.1.  Additional validation of for both symmetric and asymmetric trajectories is provided in Appendix D.   44 2.6.1.  Validation of Sinusoidal Pitching  For the sinusoidal pitching validation studies, the Reynolds number is 12,000 and the freestream Mach number is 0.015 while the pitching amplitudes were 2¡ and 4¡. First, the validation of the unsteady and average forces from the current computations is evaluated.  In order to scrutinize the validity of the time-varying forces, the peak-to-peak amplitude of the time-varying lift coefficient, CL,pp, and drag coefficient, CD,pp, are compared with the corresponding values from incompressible simulations at a Reynolds number of 12,600 by Jee and Moser (2012) and direct force measurement experiments at a Reynolds number of 16,600 by Mackowski and Williamson (2015), shown in Figure 2.12.  In the figure, circles denote 2¡ data and squares represent 4¡ data.  The results from the current simulations and those by Jee and Moser (2012) agree extremely well.  There is only marginal disagreement in the drag amplitude at the higher reduced frequencies.  The computational results disagree with direct force measurements at higher reduced frequency range, though this disagreement may be related to three-dimensional effects not present in either simulation.    Figure 2.12.  Peak-to-peak lift coefficient and peak-to-peak drag coefficient for 2¡ and 4¡ pitching amplitudes at a Reynolds number of approximately 12,000.  Data from the literature is included for comparison.  45 The average thrust coefficient validation of the current simulations is shown in Figure 2.13.  Since the average lift is nearly zero for sinusoidal pitching (10-4 for highest k), it is not considered here.  The uncertainty in thrust coefficient in the figure based on the various computational convergence studies performed is within the symbol thickness.  Data from incompressible computations by Ramamurti and Sandberg (2001), compressible computations by Young and Lai (2004), and control volume thrust estimation experiments by Bohl and Koochesfahani (2009; UnPb) are included for comparison as they are considered more reliable data from the literature.  Overall, the thrust coefficient data from the current simulations agree very well with the data from the literature for both amplitudes, specifically Young and Lai (2004) and Bohl and Koochesfahani (UnPb).  The crossover reduced frequency for both amplitudes are also comparable.  The crossover reduced frequency from the current results is 8.15 and 4.08 for 2¡ and 4¡ amplitudes, respectively, and are known to within 3%.  According to Young and Lai (2004) and Bohl and Koochesfahani (UnPb), the crossover reduced frequency for the two amplitudes are closer to 7.7 and 4.0.  The error based between these crossover reduced frequencies with the current results for the two amplitudes is 5.5% and 2%, respectively.   46  Figure 2.13.  Average thrust coefficient for 2¡ and 4¡ pitching amplitudes at a Reynolds number of approximately 12,000.  Data from the literature is included for comparison.  Next, validation of the overset grid approach is discussed for the 2¡ pitching amplitude only.  Visual comparison between the vortical structures produced by the present computations with experiments by Bohl and Koochesfahani (2009) and incompressible CFD results by Jee and Moser (2012) is shown in Figure 2.14 for reduced frequencies of approximately 5.2 and 12 and Reynolds number of 12,600.  In the figure, the airfoil is at zero angle of attack and is pitching down.  There is good agreement between the three studies, particularly between the two computational results, which looked nearly identical.  In all three investigations, a traditional von K⁄rm⁄n vortex street is produced at the k = 5.2 and a reverse von K⁄rm⁄n street at the k ! 12.  There is also agreement in the streamwise spacing of the vortices of same sign for both reduced frequencies.  The transverse spacing between vortices of opposite sign from the two computations agrees much better with each other than they do with the experiment at k = 5.2, though Jee and Moser (2012) hypothesized the difference is due to test-section wall effects not present in the simulations.  At k ! 12, there is better agreement in the transverse spacing among  47 all three studies.  Both computations are also able to resolve the connecting braids between vortices at k = 5.2 more effectively than the experiment, most likely due to better resolution.  However, these braids visually disappear in the simulations of Jee and Moser (2012) by one chord length downstream of the trailing edge whereas these braids are resolved much further downstream in the current computations.  Figure 2.14.  Instantaneous spanwise vorticity field (#z/cU$) comparison at a reduced frequency of approximately 5.2 (left column) and 12 (right column) between current compressible computations with experiments by Bohl and Koochesfahani (2009) and incompressible computations by Jee and Moser (2012).  The Reynolds number is approximately 12,000 and the pitching amplitude is 2¡ amplitude.  The airfoil is at zero angle of attack and is pitching down.  Top and middle figures courtesy of Cambridge University Press and Elsevier.  Next, quantitative validation based on vortex properties is presented.  The vortex properties considered were the following: transverse spacing (b), streamwise spacing (a), the convection speed (UC), the peak spanwise vorticity (#z,peak), vortex circulation ((), and core radius (rc). The location of the vortex core was defined as the centroid of vorticity within the vortex core, defined by equation (2.23).  The transverse spacing was the lateral distance between vortex centroid of opposite sign at the measurement location.  The convection speed was  48 computed by numerically differentiating the centroid position with respect to time.  The streamwise spacing was calculated by relating the spacing to the vortex convection speed using the relationship a = UC/f, where f was the pitching frequency.  Although this method for calculating streamwise spacing differed from Bohl and Koochesfahani (2009), it allowed the streamwise spacing to be tracked over the entire domain.    bs'ıłbıłhKcs'ıłcıł (2.23) For computing vortex circulation and core radius, Bohl and Koochesfahani used a fixed spatial radius of 1 cm (about 0.083c or 2.5 core radii from the peak) to minimize the contribution of neighboring vortices or braid vorticity and a vorticity cutoff of 2s-1 (or 2.3c/U$) to remove low-level noise.  In order to properly compare the current results with those by Bohl and Koochesfahani, the same fixed spatial radius, vorticity threshold, and measurement location were used with the computational data.  The vortex circulation was calculated by the area integral of vorticity, defined by equation (2.24).  Following Bohl and Koochesfahani (2009), the core radius was computed using the radius of gyration of the vorticity, given by equation (2.25)   œ'ıłubuc (2.24)  rs'rCıłubucıłubuc (2.25) A comparison of the vortex properties versus reduced frequency for 2¡ amplitude is shown in Figure 2.15 between the current computations and data from the literature (Bohl and Koochesfahani, 2009; Naguib, et al., 2011).  The measurement location in the work by Bohl and Koochesfahani was x/c = 0.5 whereas it was at x/c = 1.0 in Naguib, et al., (2011).  Vortex  49 property data was unavailable from Jee and Moser (2012) for comparison.  However, the locations of the vortices were extracted from their instantaneous vorticity fields by this author using an in-house code.  The a and b spacing were then computed based on the extracted vortex locations, and the convection speed was calculated using a and f.  The core radius from Naguib, et al., (2011) was taken as the 1/e point for a Gaussian distribution of vorticity.   As seen in the figure, there is excellent agreement overall between the current simulations and data from the literature in both trend and magnitude of the vortex properties.  At reduced frequencies less than approximately 6, there is slight disagreement in the b spacing (for example, refer to visual discrepancy in k = 5.2 in Figure 2.14).  However, the variation of b with k is steep at these reduced frequencies and slight experimental uncertainty in the value of k could cause better data collapse.  At higher reduced frequencies, the b spacing all approach nearly the same asymptotic value.  There is also close agreement in the b spacing crossover reduced frequency (switch from traditional to reverse von K⁄rm⁄n street) between the current computations, Bohl and Koochesfahani (2009), and Jee and Moser (2012), as seen in Table 2.3.  The a spacing and convections speed are nearly identical between the current computations and the literature.  The current results produce higher peak vorticity than in the literature, but this is most likely a result of the greater spatial resolution of the simulations.  The vortex circulation results show good collapse with Bohl and Koochesfahani (2009) at k " 8 while there is some discrepancy at k > 8.  At the highest reduced frequency, this difference is approximately 0.04cU$, or 10%.  The circulation results are close to the circulation from Naguib, et al., (2011) over the entire k range.  The core radius show disagreement in terms of trend with k compared to results of Bohl and Koochesfahani (2009), though this comparison is deceptive due to the scaling in the figure.  For example, the difference between the computational results and Bohl and Koochesfahani (2009)  50 magnitude at the highest k was approximately 14%.  In the current results, the core radius continually decreases from 0.033c to slightly less than 0.03c over the range of reduced frequencies while the core radius remains Òfairly constantÓ for Bohl and Koochesfahani (2009), with a value of rc/c = 0.032-0.034 at the measurement location in the experiment.    Figure 2.15.  Vortex properties versus reduced frequency for Reynolds number of 12,000 and 2¡ amplitude.  Data from available literature is included for comparison.  Properties calculated at x/c = 0.5; except for Naguib, et al., (2011), which are done at x/c = 1.0.  Table 2.3.  Vortex transverse spacing crossover reduced frequency for Re ! 12,000 and !o = 2¡. Author kcr (Visual) x/c kcr (at x/c) Current 5.8 0.5 6.2 Bohl and Koochesfahani (2009) 5.7 0.5 5.9 Naguib, et al., (2011) 7.1 1.0 7.1 Jee and Moser (2012) 5.7-6.7 0.5 6.4     51 2.6.2.  Validation of Asymmetric Pitching In order to validate the present computations for the asymmetric trajectory study, comparison was made between the computed wake structure and the incompressible flow experiment by Koochesfahani (1989).  In the experiment, the vortical pattern was visualized using dye (see Figure 2.16) at a Reynolds number of approximately 12,000 and reduced frequency of 6.68.  The flow visualization showed a single counter-clockwise rotating vortex and two clockwise rotating vortices shed during each pitching cycle for S = 38%, which evolved very asymmetrically.  The opposite behavior was observed for S = 61% in the experiment.  The computational result at a freestream Mach number of 0.005 for S = 38% is illustrated in Figure 2.16 with the dye-visualization by Koochesfahani (1989).  There are noticeable similarities and differences between the experimental and computational results.  The computational result is able to reproduce very similar flow features as observed in the experiment, such as the single counter-clockwise vortex and two clockwise vortices per pitching cycle.  The computation also produces similar trajectories and re-arrangement of the vortices.  However, there are quantitative differences in the orientation of the vortices, particularly the clockwise vortices. The negative vortices are more vertically aligned in the computation than the experiment.  52  Figure 2.16.  Instantaneous spanwise vorticity field #z/cU$ comparison between incompressible dye-visualization experiment (top) and compressible computation (bottom) at a reduced frequency of 6.68 and S = 38%.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  Top image courtesy of AIAA.  In order to compare the results quantitatively, average and r.m.s. velocity profiles were measured one chord length downstream of the trailing edge.  The +4 trajectory was a smoothed version of trajectory from the experiment.  Figure 2.17 shows these profiles alongside the experimentally measured profiles using Laser Doppler Velocimetry (LDV) by Koochesfahani (1989) and Naguib, et al., (2011).  Included at the top of the figure are the airfoil trajectories from the experiment and the computation.  The computational result agrees qualitatively with the experimental profiles.  Both the simulation and experiment show a velocity excess below the wake centerline and a velocity deficit at the wake centerline in the average streamwise velocity profile.  There is also quantitative agreement.  At y/c < 0.027, the average u profile from the computation collapse onto the experimental data.  Moreover, the computation also reproduces ÒsimilarÓ r.m.s. velocity profiles.  However, the transverse region of of the momentum deficit (highlighted in yellow in the figure), magnitude of r.m.s. velocity, and location of the peaks and troughs in the fluctuations do not agree precisely.  These disagreements are the result of  53 differences in the properties of the vortices, such as their circulation, relative spacing, and core size, at the measurement location.  It should be noted that evolution of vortical structures where there are several vortices shed per cycle are sensitive to initial conditions (Koochesfahani, 1989).  Additional attempts at isolating the cause of the discrepancy are provided in Appendix D.  It was found that uncertainty in the freestream velocity of the experiment could be related to the discrepancy between the computation and experiment.  Additional flow features such as slightly non-zero angle of attack, uncertainty in the actual pitching amplitude, three-dimensional effects, and test-section wall effects could also contribute to the discrepancies.  Further attempts at trying to reproduce the experimental result would be an exhaustive exercise and was not pursued further.    Figure 2.17.  Average and r.m.s. velocity profiles measured at x/c = 1.0 using analytical pitching trajectory with S = 38% and k = 6.68. Trajectory comparison between analytical formula and experimental data is included on top for reference.    54 2.7.  Final Computational Considerations 2.7.1.  Time-Periodic Convergence For all pitching airfoil computations, the entire domain rotated rigidly for twenty or more pitching cycles.  To determine if the flow reached a time-periodic state, two criteria were used.  The first criterion was that the cycle-to-cycle change, in the peak and average force values dropped below 1%.  The second criterion was that the cycle-to-cycle variation of instantaneous velocity and vorticity transverse profiles at x/c = 3 in the wake dropped below 1% of U$ and max(#z(y)), respectively.  Once these criteria were satisfied, a single cycle was computed in which evenly-spaced instantaneous phases were extracted from the computation.  The number of phases extracted for the sinusoidal and asymmetric cases were 96 and 192, respectively. 2.7.2.  Grid Remapping  Since the computational domain rotated rigidly, the wake data was remapped onto a non-moving, Cartesian grid with uniform spacing 1x = 1y = 2.5x10-3c.  The remapping was done by bilinear interpolation using the MATLAB in-built function griddata.  Remapping the data allowed for computing statistical quantities in the wake, wake profile measurements, and vortex tracking.  The vorticity data presented in this document were computed using a fourth-order accurate explicit finite difference scheme.  Averaging of the velocity and vorticity fields in the wake was done by taking the mean at every grid point in the remapped data using either 96 or 192 instantaneous phases.  The r.m.s. of the velocity and vorticity was calculated at every grid point according to  {ƒ†Ÿ'{D{C‡x…<„ (2.26)  55 CHAPTER 3.  EFFECT OF REYNOLDS NUMBER ON AERODYNAMIC LOADS   One of the primary motivations of this work is to study the effect of Reynolds number on the unsteady and average aerodynamic loads and vortical structure in the airfoil produced by an airfoil pitching sinusoidally.  Six Reynolds numbers between 2,000 and 22,000 were used in this work.  They were selected to remain low enough such that laminar flow could be assumed.  Pitching amplitudes varied between 2¡ to 10¡ in 2¡ increments while the reduced frequency ranged from 0 (static airfoil) to 16.0.   These reduced frequencies were chosen to provide a thrust curve covering both negative (i.e., drag) and positive thrust values for each Reynolds numbers.  The freestream Mach number is 0.015 for all data presented in this chapter.  This chapter is organized as follows.  Section 3.1 deals with the Reynolds number effect on the average thrust, Section 3.2 the thrust crossover frequency, Section 3.3 presents the efficiency, Section 3.4 the load fluctuation, and Section 3.5 compares the current computations with classical linear theory for an oscillating flat plate by Garrick (1936; 1957).  3.1.  Average Thrust Coefficient In order to establish a baseline for studying the Reynolds number effect on the loads on a pitching airfoil, it is instructive to first quantify the Reynolds number effect on the static airfoil loads.  These data will later be used to scale the thrust coefficient on the pitching airfoil.  Due to symmetry in the kinematics, the average lift coefficient was zero and thus is not presented here.  The drag coefficient, as expected, depends on Reynolds number.  The time-averaged total (/€\=), pressure-induced (/€\=t), and friction-induced (/€\=y) drag coefficients are shown in Figure 3.1.  As seen in from the figure, the average drag coefficient, pressure drag, and friction drag decrease  56 non-linearly as the Reynolds number increases.  Although the friction drag is greater than the pressure drag over the entire Reynolds number range, its contribution to the total drag decreases from 70% to 55% as the Reynolds number increases from 2,000 to 22,000.  The total, pressure, and friction drag coefficient variation follow an A/ReB power-law reasonably well; the power-law was computed using a least-squares procedure.  The R2 values of the fits were greater than 98% for all three curves.  Remarkably, the friction drag for the airfoil had a Re exponent of -0.558, which is within approximately 1% of the exponent for the flat plate value of -0.5 from laminar boundary layer theory (Blasius, 1908).  The pressure and shear stress distributions leading to these net effects on the average drag are discussed in Appendix E.  Figure 3.1.  Total, pressure-induced, and friction-induced average drag versus Reynolds number for static NACA 0012 at ! = 0¡. The broken lines represent a least squares power-law fit to the data.  Considering the oscillating airfoil, the average thrust coefficient, CT, versus reduced frequency at three Reynolds numbers for the 2¡ pitching amplitude is shown in Figure 3.2.   57 Positive CT indicates thrust while negative CT corresponds to drag.  At k = 0, the airfoil experiences the drag behavior shown in Figure 3.1.  As the reduced frequency increases, the drag decreases in magnitude for a given Reynolds number.  The drag eventually reaches zero and switches to thrust as the reduced frequency increases further.  This behavior is quantitatively similar for all Reynolds numbers, although CT values increased with increasing Re.  As a consequence, the reduced frequency when drag switches to thrust depends on Reynolds number.  This dependency is discussed in greater detail later.  Although this conclusion is intuitive since viscous drag must be overcome to achieve thrust (Koochesfahani, 1989) and is consistent with the static airfoil drag data in Figure 3.1 and observation by others (e.g. Liu and Kawachi, 1999; Buchholz and Smits, 2008; Visbal, 2009; Borazjani and Sotriopoulos, 2008), quantitative data showing the extent of this dependence over the range of reduced frequencies and amplitudes presented in this work have not been available in the literature.  Figure 3.2.  Average thrust coefficient versus reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡.  58 One factor that influences the Reynolds number dependency is the baseline static airfoil drag at each Reynolds number.  By removing the static drag from the average thrust coefficient data, the thrust contribution from the unsteady motion alone, CT,u, can be examined.  As shown in Figure 3.3, CT,u increases as both reduced frequency and Reynolds number increase.  The difference in CT,u among the Reynolds numbers at k < 2 is small.  At high reduced frequency, the difference among the Reynolds numbers is much greater.  It is clear from shifting the CT vs. k curves in Figure 3.2 that the Reynolds number dependence on thrust is not simply due to the static airfoil drag coefficient, CD,o.  Figure 3.3.  Average thrust coefficient, with static drag removed, versus reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡. It is instructive to assess the contribution to CT by the pressure and shear stress distributions independently according to /0'/0t:/0y.  The pressure and friction components for the three Reynolds numbers, with their corresponding contributions by the static airfoil drag  59 removed, are presented in Figure 3.4 in order to isolate the influence of the unsteady motion.  Both components increase in magnitude as the reduced frequency increases for a given Reynolds number, which is consistent with observations by Ramamurti and Sandberg (2001) and Borazjani and Sotriopoulos (2008) that both the pressure and shear stress contributions to CT increase in magnitude with k.  The pressure contribution shows weak dependence on Reynolds number.  As the Reynolds number increases, there is a slight increase in the pressure contribution.  As k increases, the difference in /0t between the Re = 2,000 and 22,000 grows.  Theoretically, /0t should represent the inviscid solution for this pitching amplitude and thus be independent of Reynolds number.  However, slight variations in the streamline pattern due to the presence of the boundary layer could alter the resulting pressure distribution at the surface, and hence show a weak dependence of Re on /0t.  The friction contribution shows a strong dependence on Reynolds number.  The data show that the friction force generated by the pitching motion is non-linearly dependent on Reynolds number.  In fact, it scales with the static airfoil friction drag, /€\=y.  These data quantitatively show that, as the Reynolds number increases, an oscillating airfoil has to overcome both less static drag and less unsteady-motion-produced friction drag to produce net thrust.   60  Figure 3.4.  Average thrust coefficient components, with static drag components removed, versus reduced frequency for different Reynolds numbers (color).  The pitching amplitude is 2¡. The effect of pitching amplitude on thrust coefficient variation versus reduced frequency is shown in Figure 3.5 for Re = 2,000 and 22,000.  As the pitching amplitude increases, the thrust coefficient increases for both Reynolds numbers, though this is expected based on similar trends previously observed by Koochesfahani (1989) and Ramamurti and Sandberg (2001) at Re = 12,000.  At an amplitude of 2¡, the increase in thrust with k is relatively shallow.  As the amplitude increases, the thrust variation with k increases in steepness.   For the particular case of 10¡ amplitude, there is a slight increase in drag for both Reynolds numbers at k = 1 before the drag decreased as k increases.  Interestingly, this ÒdipÓ in CT is greater for Re = 2,000 than Re = 22,000.  Since the thrust increases as the amplitude increases, it is observed that there is amplitude dependence on the drag-to-thrust crossover reduced frequency.  This is quantitatively discussed in the next section.  61  Figure 3.5.  Average thrust coefficient versus reduced frequency for different pitching amplitudes (symbol).  Data for Re = 2,000 on left and for Re = 22,000 on right. Multiple authors (Triantafyllou, et al., 1993; Liu and Kawachi, 1999; Ramamurti and Sandberg, 2001) have shown that the Strouhal number (St = Lf/U, where f is the flapping frequency, L is the reference length, and U is the reference velocity) is a very important parameter for scaling the propulsive characteristics of an oscillating airfoil.  To account for the trailing edge amplitude, the length scale is selected to be the trailing edge peak-to-peak amplitude (i.e., L = ATE) for a pitching airfoil.  By plotting the average thrust coefficient data against Strouhal number for different amplitudes, the thrust data become nearly amplitude independent for the two Reynolds numbers (see Figure 3.6).  The 2¡-6¡ data show good collapse over the entire Strouhal number range for each Reynolds number.  The 8¡ amplitude data for Re = 2,000 also seem to collapse on this curve and only deviate at high St.  At 10¡ amplitude, the Re = 2,000 data deviate primarily at low k.  The 8¡ and 10¡ amplitude data for Re = 22,000, however, produce slightly lower values of CT.  Even though the thrust values are slightly lower at the higher amplitudes at a given Strouhal number, the variation of CT with St is fairly uniform across the different amplitudes considered in the present work.   62  Figure 3.6.  Average thrust coefficient versus Strouhal number for different pitching amplitudes (symbols).  Data for Re = 2,000 on left and for Re = 22,000 on right. A clue about the decrease in CT as the amplitude increases for a constant Strouhal number may be elicited from Figure 3.7, which shows the instantaneous vorticity field around the airfoil for the 10¡ amplitude case at two Reynolds numbers and reduced frequencies of 1 and 4.  In the figure, the airfoil is at its maximum angle of attack (%/T = 0.25).  At k = 1, both Reynolds numbers show a separation region that covers the majority of the aft side of the airfoil.  The higher of the two Reynolds number also shows a great deal of instability in the boundary layer. At a reduced frequency of 4, both Reynolds numbers show separation near the leading edge.  The boundary layer separation at Re = 2,000 remains relatively flat while Re = 22,000 produces a very distinctive vortex.  The leading edge separation for both Reynolds number propagates along the airfoil as the airfoil pitches the opposite way.  This vortex is eventually destroyed when the airfoil pitches in the same direction as the phase in which the vortex formed.   63  Figure 3.7.  Instantaneous spanwise vorticity field (#zc/U$) for two reduced frequencies (1.0 and 4.0) and two Reynolds numbers (2,000 and 22,000).  The pitching amplitude is 10¡ and the airfoil is at its maximum angle of attack position and is about to pitch down (%/T = 0.25). 3.2.  Thrust Crossover Frequency As discussed earlier, one consequence of the Reynolds number effect on the thrust force is that the reduced frequency when the drag switches to thrust depends on Reynolds number.  In order to quantify this dependency, the thrust crossover reduced frequency, kcr, was determined for each amplitude and Reynolds number.  The crossover reduced frequency was found by fitting a polynomial ranging from 1st-order to 3rd-order using two to four points, respectively, closest to CT = 0.  Figure 3.8 shows the thrust crossover reduced frequency variation with respect to Reynolds number for each amplitude.  For the 2¡ amplitude case, the crossover reduced frequency decreases from 14.3 to 6.9 as the Reynolds number increases from 2,000 to 22,000.  When the pitching amplitude increases from 2¡ to 4¡, the crossover reduced frequency decreases nearly uniformly by a factor of 2.  As the pitching amplitude increases further, the crossover  64 reduced frequency decreases at a given Reynolds number towards the inviscid theory prediction of kcr = 1 for an oscillating flat plate (Garrick, 1936) and the inviscid prediction of kcr = 3 for NACA 0012 pitching with 2¡ amplitude (Ramamurti and Sandberg, 2001).   Figure 3.8.  Thrust crossover reduced frequency versus Reynolds number for different pitching amplitudes (symbols).   Since the thrust coefficient showed good data collapse when scaled using Strouhal number based on trailing edge amplitude, the thrust crossover reduced frequency was rescaled as a thrust crossover Strouhal number according to the relationship B4sƒ'¡0¢£sƒPA.  The crossover Strouhal number dependency on Reynolds number for different amplitudes is presented in Figure 3.9.  There is a reasonable collapse, particularly at low amplitudes, of the crossover frequency data when rescaled as a crossover Strouhal number.  At !o " 6¡, the crossover Strouhal number data nearly collapse onto a single curve that is well represented by an A/ReB power law fit.  As the Reynolds number increases from 2,000 to 22,000, the crossover  65 Strouhal number decreases from approximately 0.23 to 0.11 in the low amplitude limit.  As the amplitude increases to 8¡ and 10¡, there is very noticeable deviation from the Òlow amplitudeÓ power-law fit, particularly at the higher Reynolds numbers.  At Re = 22,000, the crossover Strouhal number increases from 0.11 to 0.15 as the pitching amplitude increased from 2¡ to 10¡.  This may be related to the behavior seen in Figure 3.7, where there is noticeable boundary layer separation at low reduced frequency (k = 1) and leading edge separation at the higher reduced frequency (k = 4) at these high amplitude conditions.   Figure 3.9.  Thrust crossover Strouhal number versus Reynolds number for different pitching amplitudes (symbols).   Comparison of the crossover Strouhal number is made with reliable data found in the literature shows good agreement (see Figure 3.10).  Although the majority of data from the literature is for Reynolds number of 12,000 for a NACA 0012, thrust estimation experiments using the mean streamwise velocity component measured in the wake were done on a tear-drop  66 shaped by Godoy-Diana, et al., (2008) at a Reynolds number of approximately 1,200.  Although it has been shown that the thrust estimated using the mean u-profile only over-predicts the thrust (Streitlien, et al., 1998; Bohl and Koochesfahani, 2009), the crossover Strouhal number from Godoy-Diana, et al., (2008) agrees with the current data.    Figure 3.10.  Thrust crossover Strouhal number versus Reynolds number data comparison with literature.  Only the 2¡, 6¡, and 10¡ from the current data set are included.  Note that Godoy-Diana, et al., (2008) used a tear-drop (TD) shaped airfoil. These findings appear to be consistent with observations found in nature among flyers and swimmers.  RŠsen, et al., (2004), performed measurements on a thrush nightingale to study its wing beat kinematics.  The nightingale had a mean chord length of 4.8 cm.  The authors found that the wing beat amplitude and flapping frequency had values of 0.08±0.003 m and 14.4±0.4 Hz over flight speeds ranging from 5 m/s to 10 m/s.  The corresponding Reynolds number based on mean chord length are 16,000 and 32,000, respectively.  The reduced frequency based on mean chord length decreased monotonically from approximately 0.4 to 0.2 as the flight speed  67 increased over the same flight speed range.  These operating reduced frequencies are significantly lower than the current results, which is most likely due to the more complex flapping kinematics and larger peak-to-peak amplitudes used by the nightingale compared to the pitching airfoil in this study.    Borazjani and Sotriopoulos (2008) performed incompressible computations on a mackerel performing carangiform undulations (side-to-side oscillations from a view looking down onto the swimming mackerel) to study the effect of Reynolds number on the operating Strouhal number after noting that there exists a complex relationship between operating Strouhal number and Reynolds number based on the findings by Lauder and Tytell (2006) that the operating Strouhal number of Pacific salmon (Re = 104) is not the Strouhal number of optimal efficiency of 0.25-0.35 (Triantafyllou, et al., 1993) and changes according to swimming speed (i.e., Re).  The authors considered Reynolds numbers of 300 and 4,000, along with the inviscid limit.  Borazjani and Sotriopoulos (2008) showed that as the Reynolds number increased, the thrust coefficient increased.  As the Reynolds number increased from 300 to 4,000, the crossover Strouhal number decreased from 1.08 to 0.6, which are substantially higher than in the current work.  These values were also much higher than prediction of 0.26 by the authorsÕ inviscid simulation.  Interestingly, the inviscid Strouhal number falls within the optimal Strouhal number range according to Triantafyllou (1993).  Although the scaling the thrust coefficient with the Strouhal number produced a dramatic collapse for each Reynolds number, the data formed into distinctive bands that depended on Reynolds number (see Figure 3.11).    68  Figure 3.11.  Average thrust coefficient versus Strouhal number for different Reynolds numbers (color) and pitching amplitudes (symbols).  Only the 2¡, 6¡, and 10¡ are included. Using both airfoil static drag from Figure 3.1 and the crossover Strouhal number from Figure 3.9, it is possible to rescale the data in Figure 3.11 to become nearly Reynolds number independent, as shown in Figure 3.12.  Dividing the thrust data for each Reynolds number by its corresponding static drag value starts each curve at -1 on the ordinate.  Dividing the Strouhal number by the crossover Strouhal number shifts each thrust curve on the abscissa so that drag switches to thrust at St/Stcr = 1.  Overall, there is very good data collapse among the three Reynolds numbers for the different amplitudes.  The 2¡-6¡ data are well-represented by a cubic polynomial fit of the form 0.136x3 + 1.140x2 Ð 0.295x Ð 0.986.  As the pitching amplitude increases, the 8¡ and 10¡ data begin to show some deviation at St/Stcr # 1.5.    69  Figure 3.12.  Average thrust coefficient versus the Strouhal number scaled by the crossover Strouhal number as a function of Reynolds number (color) and pitching amplitude (symbol).  The thrust coefficient is scaled by the static airfoil drag. 3.3.  Propulsive Efficiency For completeness, the effect of Reynolds number on propulsive efficiency is now discussed.  The propulsive efficiency is defined by equation (3.1) and represents the ratio between the propulsive energy and the total energy input into the system.  In terms of standard non-dimensional quantities, the propulsive efficiency is the average thrust coefficient divided by the average input power coefficient.   -.'/0/1.'23/0454670623/849:/;4+546706 (3.1)  70 Most studies found in the literature (Triantafyllou, et al., 1993; Liu and Kawachi, 1999; Mackowski and Williamson, 2015) indicate that the propulsive efficiency scales according to amplitude when plotted versus the Strouhal number based trailing edge-amplitude.   The propulsive efficiency for the three Reynolds numbers from the current computations are plotted versus Strouhal number in Figure 3.13.  The propulsive efficiency data fall into three distinctive curves based on Re.  As the Reynolds number increases, the propulsive efficiency increases.  At a Reynolds number of 2,000, the efficiency has a maximum value of approximately 6.5% to while it is closer to 11% for a Reynolds number of 22,000.  These maximum efficiency values are much lower than the linear theory prediction by Garrick (1936) of 50% for a flat plate pitching about its quarter-chord and the efficiency of plunging airfoils of around 20% (Triantafyllou, et al., 1993).  It is noteworthy that the data do not show a well-defined peak, as seen in Figure 1.8 in Chapter 1, for example.  As the efficiency increases, it reaches a plateau.  At Re = 7,000 and 22,000, it is discernable that is there is a very slight drop in efficiency at St > 0.32, which falls within the optimal Strouhal range of 0.25-0.35 (Triantafyllou, et al., 1993).  Similar conclusions can be drawn from data presented by Mackowski and Williamson (2015) for a pitching NACA 0012 at Re = 16,600, specifically that the propulsive efficiency shows a plateau that did not greatly decrease as St increased.    71  Figure 3.13.  Average propulsive efficiency versus Strouhal number as a function of Reynolds number (color) and pitching amplitude (symbol). 3.4.  Force Fluctuation  In this section, the effect of Reynolds number on the fluctuating force on a pitching airfoil is described.  As with the average force, the static airfoil is first considered to establish a baseline.  Examination of the force time histories showed unsteady oscillations in lift for Re # 12,000.  These oscillations were the result of a wake instability at these Reynolds numbers causing the wake to roll-up into discrete vortices.  This feature will be discussed further in the next chapter.  The lift fluctuation was characterized in terms of the peak-to-peak amplitude, CLpp, and the natural frequency, fn.  These data are shown in Figure 3.14, where the natural frequency is non-dimensionalized by U$/$c to produce an equivalent natural shedding reduced frequency, kn = $fnc/U$.   Figure 3.14 shows that the lift amplitude increases nearly linearly with Reynolds number. It is also observed that there is nearly linear increase in natural shedding reduced frequency from  72 8.5 at a Reynolds numbers of 12,000 to approximately 11 at a Reynolds number of 22,000. This trend is consistent with observation by Huang and Lin (1995), who stated Òthe frequency of vortex shedding increases linearly with the increase of freestream velocity.Ó  There is also excellent agreement between the shedding reduced frequency of the current computations with the data presented by Koochesfahani (1989), Huang and Lin (1995), and Young and Lai (2004).    Figure 3.14.  Reynolds number effect of the fluctuating force characteristics for static NACA 0012 airfoil at ! = 0¡. Now we examination of the force fluctuation of the unsteady airfoil. The lift and drag coefficient fluctuation, /8} and /€}, are defined as the time-varying force with the average value removed.  Note that the average lift is zero for these cases.  The load fluctuations at four reduced frequencies ranging from 0.5 to 8 are shown in Figure 3.15 for 2¡ amplitude and three Reynolds numbers: 2,000, 7,000, and 22,000.  All three Reynolds produce very similar lift and drag fluctuation behavior.  The lift fluctuation shows a single set of peaks while there are two sets of peaks in the drag fluctuation.  As the reduced frequency increases, there is a slight increase in the  73 peak values, as well as a phase shift.  This is not surprising based on linear theory for a flat plate (Theodorsen, 1935).   Figure 3.15.  Lift and drag fluctuations for different reduced frequencies.  The Reynolds numbers are Re = 2,000 (blue), 7,000 (red), and 22,000 (yellow).  The pitching amplitude is 2¡ amplitude.  Note that the scale is different for the different plots.  The lift and drag fluctuation for Re = 22,000 shows a high frequency oscillation compared to the pitching reduced frequency of k = 0.5.  The cause of this behavior is the wake instabilities observed for the static case at this Reynolds number.  Due to the low reduced frequency, the wake undulated about the centerline while carrying the vortices formed from the natural instability.  Similar wake structure observations were made by Koochesfahani (1989) for a Reynolds number of 12,000, 0.5 Hz pitching frequency (k = 0.835), and pitching amplitude of 4¡.  Application of a Fast Fourier Transform (FFT) to the lift fluctuation reveals two harmonics (see Figure 3.16).  The first is the forcing reduced frequency (i.e., k ! 0.5).  The second was the natural shedding reduced frequency of approximately 11 seen in Figure 3.14.   74  Figure 3.16.  Magnitude of lift fluctuation FFT for the Re = 22,000 case.  Instantaneous spanwise vorticity field at %/T = 0.25 is included for visual reference. The load fluctuation can be characterized by the peak-to-peak amplitude (pp) and phase angle, /L.  The phase angle was quantified using the lift fluctuation only with respect to the angle of attack.  The lift phase angle was determined by subtracting the phase at which the minimum lift coefficient occurred from the phase at which the minimum angle of attack occurred.  The lift lagged behind the angle of attack if the phase angle was negative (i.e., increasingly negative /L yields increasing magnitude of phase lag or delay).  This is shown schematically in Figure 3.17.    75  Figure 3.17.  Schematic of the phase angle definition. The peak-to-peak amplitudes and phase angle data for three Reynolds numbers and three pitching amplitudes (2¡, 6¡, and 10¡) are shown in Figure 3.18.  As the reduced frequency increases, the lift and drag amplitudes and phase increases in magnitude for each Reynolds number.  As the Reynolds number increases, the lift and drag amplitudes slightly increase in magnitude while the phase delay slightly decreases in magnitude.  As the amplitude increases, the lift and drag amplitudes increase while the phase lag decreases in magnitude.      76  Figure 3.18.  Peak-to-peak lift amplitude (left), peak-to-peak drag amplitude (middle), and lift phase lag (right).   Different Reynolds numbers are shown with the symbol color and different pitching amplitudes are shown with symbol type. It is instructive to scale the lift and drag amplitudes based on classical linear theory.  As shown in Appendix F, the lift and drag scale with += and +=C, respectively.  The resulting scaled data are shown in Figure 3.19.  The rescaled data show good collapse at the over the majority of the k range for the three amplitudes and Reynolds numbers.  The cause for the data collapse can be inferred from Figure 3.20, which compares the total time-varying lift and drag with the pressure and shear stress contributions for the representative example of k = 6.0 and !o = 2¡ at Re = 2000 and Re = 22,000.  As expected, the variation in the lift is due to the pressure distribution.  Less expected is that the the variation in the drag is also due primarily to the pressure distribution.  The shear stress contribution also did vary slightly with time, but is much lower in magnitude than the pressure contribution.  This behavior occurs for each Re, k, and !o.  As the amplitude increases, there is a slight increase in the fluctuation of shear stress drag.     77  Figure 3.19.  Peak-to-peak lift amplitude normalized by += (left) and peak-to-peak drag amplitude normalized by +=C (right).   Different Reynolds numbers are shown with the color and different pitching amplitudes are shown with symbols.  Figure 3.20.  Lift (top) and drag (bottom) histories, showing contribution from pressure and shear stress distribution compared to the total force for Re = 2,000 (left) and Re = 22,000 (right) at k = 6.0 and !o = 2¡.  78 3.5.  Comparison with Inviscid Linear Theory Finally, it is of interest to see how the results from the viscous computations compare with inviscid theory of Theodorsen (1935) and Garrick (1936; 1957), which are presented in Appendix F.  Caution, however, must be exercised when making this comparison since the assumptions of a flat plate geometry, enforcement of the Kutta condition at the trailing edge, flat wake, and low oscillation amplitudes are not satisfied in the present computations, except for low amplitude.  By removing the static drag from the average thrust coefficient data so that only the CT,u is considered, the different cases may be directly compared with linear theory by Garrick (1957) as illustrated in Figure 3.21.  For both amplitudes, the three Reynolds numbers produce lower thrust than the linear theory over the majority of the reduced frequency range.  The thrust contribution from the pressure distribution with the static pressure drag removed showed better agreement with linear theory (see Figure 3.22).    Figure 3.21.  Average thrust coefficient with static drag removed versus reduced frequency for different Reynolds numbers (color) and pitching amplitudes (symbol).  Data shown is for 2¡ and 10¡ amplitudes.  The linear theory prediction is included for comparison.  79  Figure 3.22.  Average thrust coefficient from pressure with static pressure drag removed versus reduced frequency for different Reynolds numbers (color) and pitching amplitudes (symbol).  Data shown is for 2¡ and 10¡ amplitudes.  The linear theory prediction is included for comparison. Comparison is also made with the lift and drag peak-to-peak amplitudes and the phase lag (see Figure 3.23).  The lift and drag were normalized using the amplitudes as done in Figure 3.19.  The lift and drag scaled by the amplitude agree very well with linear theory over the majority the reduced frequencies.  The most deviation from linear theory occurs in the drag amplitude data at reduced frequencies above 10.  The phase angle agrees only at low reduced frequency.  As the reduced frequency increases, the phase delay from linear theory continue to increase towards an asymptotic limit of 180¡ while the data from the present computations appear to asymptote towards approximately 160¡.  The computations also show an amplitude dependence on the phase angle, which does not occur in linear theory.  80  Figure 3.23.  Peak-to-peak lift amplitude (left), peak-to-peak drag amplitude (middle), and lift phase lag (right).   Different Reynolds numbers are shown with the symbol color and different pitching amplitudes are shown with symbol type.   The linear theory prediction of Theodorsen (1935) and Garrick (1936) is included for reference.       81 CHAPTER 4.  EFFECT OF REYNOLDS NUMBER ON WAKE STRUCTURE  Discussed in this chapter is the effect of Reynolds number on the wake structure produced by an airfoil pitching sinusoidally.  Only the 2¡ amplitude is considered for the majority of the results presented in this chapter, with the effect of amplitude discussed at the end of the chapter.  The reduced frequency again varied from approximately 0 (static airfoil) to 16.0.  The freestream Mach number is 0.015 for all data presented in this chapter.  This chapter is organized as follows.  Section 4.1 provides a discussion of the flow behavior near the trailing edge, Section 4.2, the Reynolds number effect on the wake structure configuration, Section 4.3, the crossover reduced frequency for vortex configuration switch, and Section 4.4, the effect of oscillation amplitude on the vortex configuration. 4.1.  Flow Near Trailing Edge  The flow behavior at the airfoilÕs trailing edge may be quantified using the velocity and vorticity profiles.  ÒMeasurementsÓ are taken at the trailing edge over ±1c from the airfoil surface in the transverse direction.  Streamwise and transverse velocity components as well as spanwise vorticity profiles for different Reynolds number are shown in Figure 4.1 for the static airfoil.  The u-velocity component profiles show the velocity beginning at zero due to no-slip at the surface and increasing in magnitude to approximately U$ for each Re over the Y-extent shown.  As the Reynolds number increases, the visual thickness of the boundary layer decreases.  Another feature is the development of a reverse flow region and slight velocity deficit (see the region highlighted by the magenta box), which grows as the Reynolds number increases.  The z-vorticity profiles also show the flow reversal where there is positive vorticity.  The spanwise  82 vorticity profiles show negative values outside of the reverse flow region, which reaches a peak value near the middle of the boundary layer and decreases in magnitude towards zero as the the outer flow is approached.    Figure 4.1.  Upper surface u velocity and z-vorticity profiles measured at trailing edge for static airfoil at ! = 0¡ for different Reynolds numbers. The boundary layer thickness (BTL) is defined as the distance from the surface for the velocity to reach 99% of the outer flow velocity according to laminar flat plate theory.  Since the inviscid outer flow velocity is not known, a zero vorticity condition is used in which the boundary layer thickness is defined as the distance from the surface in which the vorticity reaches a near-zero value.  Here, the near-zero cutoff value is taken to be %zc/U$ = 0.001.  The BTL at the trailing edge for the static airfoil, wo, on the upper surface is shown in Figure 4.2.  As the Reynolds number increases, the BTL decreases according to 1/Re0.4.    83  Figure 4.2.  Boundary layer thickness as a function of Reynolds number for static NACA 0012 at ! = 0¡.  When the airfoil pitches sinusoidally, the trailing edge translates nearly vertically in the opposite direction as the leading edge (i.e., the trailing edge moves down as the airfoil pitches up).   Streamwise velocity profiles measured at the trailing edge for different Reynolds numbers are shown in Figure 4.3 for k = 6.0 and !o = 2¡ at eight, evenly spaced oscillation phases.  The selection of this reduced frequency is provided in the next section.  The transverse direction is scaled by the static airfoil boundary layer thickness in the figure to make differences in the profiles clearer.  The cyan and magenta lines denote u = 0 and u = U$, respectively, for reference.  Instantaneous vorticity fields with a field of view near the trailing edge at the selected phases are shown in Figure 4.4 for Re = 12,000 as a representative case.  As the trailing edge translates in the negative Y direction (%/T = 0.125-0.25), the u velocity in the flow outside of the boundary layer accelerates to approximately U$.  A small region of reverse flow could be seen at %/T = 0.125 for Re = 22,000 and disappears by %/T = 0.25.    84 As the trailing edge translates in the positive Y direction during the pitch-down (%/T = 0.25-0.75), the u velocity continues accelerating for each Reynolds number beyond U$ until the outer flow reaches its maximum value at %/T ! 0.50, then decelerates.  Flow reversal is observed at %/T = 0.75 and the outer flow again reaches a value of approximately U$.  The region of flow reversal increases in magnitude as the Reynolds number increases.  Flow reversal was also observed experimentally by Bohl (2002) for a NACA 0012 airfoil pitching sinusoidally at k = 8.8 and Re = 12,600 based on single component MTV tagging lines near the trailing edge.  Within the phase range 0.75 < %/T < 0.875 as the airfoil pitched up, the flow reversal occurs over a wider region in the transverse direction.  Also observed is a secondary layer at %/T < 0.875, which increases in as Reynolds number increases.  This secondary layer is a small region of high magnitude negative vorticity near the trailing edge that forms below the positive vortex, as indicated in Figure 4.4 for the case of Re = 12,000.  At %/T < 1.0, the flow reversal has decreased in magnitude while the secondary layer decreases in magnitude as the Reynolds number increases.  In general, the streamwise velocity increases in as the Reynolds number increases during the pitch-down (0.50 < %/T < 0.75).  As the airfoil slows down and switches directions, there is a phase lag in u that decreases as the Reynolds number increases.  This is not unexpected, however.  85   Figure 4.3.  Streamwise velocity boundary profiles, where Y is scaled by the static airfoil BLT, at eight phases of the oscillation cycle for different of Reynolds numbers.  The reduced frequency 6.0 and pitching amplitude is 2¡.  Measurements are made at the trailing edge on the upper surface.  The coordinate system starts at the airfoil surface.  The airfoil position at every other phase is included for reference.  The cyan line denotes u = 0 and the magenta line represents u = U!.     86  Figure 4.4.  Instantaneous spanwise vorticity field ("zc/U!) at 8 phases of the oscillation for k = 6.0, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The green line is the measurement location and the white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  87 4.2.  Wake Structure Configuration for 2¡ Pitching The Reynolds number effect on the wake structure for the static airfoil is examined in Figure 4.5.  Note that the coordinate system starts at the trailing edge at x/c = 0.  As the Reynolds number increases from 2,000 to 7,000, the vorticity shed into the wake remains stable.  The width of the wake decreases as the Reynolds number increases.  As the Reynolds number increases from 7,000 to 12,000, the wake transitions from a stable to unstable state.  The instability in the wake causes the wake to roll-up into discrete vortices.  These vortices form into a traditional von K⁄rm⁄n vortex street that is reminiscent of the flow behind a cylinder at low Re.  As the Reynolds number increases from 12,000 to 22,000, the vortices visually become smaller.  The wavelength and vertical spacing of the vortices also visually decrease.  Another feature is that the formation location of the vortices moves upstream towards the trailing edge as the Reynolds number increases. This Reynolds number effect is consistent with flow visualization by Huang and Lin (1995) who saw a stable wake at Re = 3,195 and an unstable wake at Re  8,000.  A stable wake was also observed in flow field images by Liu and Kawachi (1999) for Re = 7,200.   88  Figure 4.5.  Instantaneous spanwise vorticity fields (!zc/U") for static NACA 0012 at ! = 0¡ and different Reynolds numbers (increasing from the top down).  The organization of the vortices is quantified in terms of the transverse spacing (b) and streamwise spacing (a) to determine the Reynolds number effect on the natural flow of a steady airfoil.  Refer to Figure 1.6 in Chapter 1 for a schematic.  The values of an and bn, which are the streamwise and transverse spacing of the naturally formed vortices in the wake, are shown in Figure 4.6.  The location of the core center of a vortex is defined as the spatial position of the centroid of vorticity contained within the vortex.  Measurements were taken three chord lengths from the airfoil trailing edge such that the vortices were fully formed visually (see Figure 4.5).  As seen in Figure 4.6, all three Reynolds numbers produce negative b due to the formation of a traditional von K⁄rm⁄n street in the wake from natural shedding of the shear layer.  As the Reynolds number increases, the transverse spacing decreases in magnitude.  Physically, this means the positive and negative vortices moved toward the wake centerline (y/c = 0) as the  89 Reynolds number increases.  The streamwise spacing decreases as the Reynolds number increases.    Figure 4.6.  Wake vortex transverse (left) and streamwise (right) spacing dependence on Reynolds number for static NACA 0012 at ! = 0¡. Referring to Figure 4.7, the effect of Reynolds number on the pitching airfoil with 2¡ amplitude is considered.  When the airfoil is pitching with a reduced frequency of 2.0, the shear layer undulates about the wake centerline.  At Re = 2,000, vorticity carried by the wake rolls up and causes the undulating wake to form into a von K⁄rm⁄n street.  At the three higher Re values, the undulating wake carries the natural instabilities that form in the shear layer.  As the reduced frequency increases towards 5.2, the wake structure for Re = 2,000 maintains a von K⁄rm⁄n street whereas more complicated patterns form for the higher Reynolds numbers.  By k = 5.2, all four Reynolds numbers produce a vortex array with a single pair of alternating sign vortices shed each pitching cycle.   90   Figure 4.7.  Instantaneous spanwise vorticity field (!zc/U") for #o = 2¡ and different Reynolds numbers (increasing left to right) and reduced frequencies (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up ($/T = 0.0).  The field of view extends over the range 0 ! x/c ! 3.0 and -0.5 ! y/c ! 0.5.    91 The Reynolds number dependence of the vortex array in the wake for the special case of k = 6.0 is shown in Figure 4.8.  The choice of k = 6.0 is that a switch from von K⁄rm⁄n to reverse von K⁄rm⁄n streets is clearly visible as the Reynolds number increases.  The airfoil is at zero angle of attack and is pitching up (!/T = 0.0).  A single pair of alternating sign vortices shed past the trailing edge each cycle and thin strips of vorticity, called braids, connect the vortex cores.  These braids originated in the boundary layer and were a result of vorticity continuously shedding into the wake before and after the shedding of the primary vortex that forms at the trailing edge.  These braids are weak compared to the vortex core strength and dissipate much more rapidly than the vortices themselves.  As the vortices convect downstream, these braids wrap around the vortex and reinforce it.   The lowest Reynolds number of 2,000 produces large vortices with low levels of vorticity in a von K⁄rm⁄n street orientation (drag).  The connecting braids are visually thick compared to the higher Reynolds numbers. At Re = 7,000, the vortices become smaller and more concentrated compared to Re = 2,000.  Although still in a von K⁄rm⁄n street orientation (drag), the transverse spacing between vortices decreases in magnitude compared to the Re = 2,000.  At Re = 12,000, the vortices are slightly in a reverse von K⁄rm⁄n street arrangement.  At Re = 22,000, the vortex array is clearly in a reverse von K⁄rm⁄n street pattern.  The distribution of vorticity within the vortex core at Re = 2,000 is somewhat non-Gaussian (see Figure 4.9), where the transverse width is greater than the streamwise width.  As the Reynolds number increases, while the vorticity distribution becomes more Gaussian.  92  Figure 4.8.  Instantaneous spanwise vorticity field ("zc/U#) for k = 6.0, $o = 2¡, and different Reynolds numbers (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (!/T = 0.0).  Figure 4.9.  Vorticity profiles normalized by peak vorticity at x/c ! 1.0 for different Reynolds numbers at k = 6.0 and $o = 2¡.  93 The evolution of the vortices was quantified by tracking the vortex centroid in space and time (Bohl and Koochesfahani, 2009).  An initial estimate was made of the vortex location by identifying the location of peak vorticity.  In order to remove false peaks near the trailing edge due to strong vorticity from the shear layer, the !" of each peak was computed using equation (4.1) (Graftieaux, et al., 2001).   !"#$%&'()*+),+(-./'()-+),+(0123 (4.1) In the above expression, '() represents the position vector of point of interest, P, and a neighbor point, M, +) is the Cartesian velocity vector of point M, +4 is a spatial average of the velocity at point P, ./ is the unit vector in the spanwise direction, and N is the number of points of interest that form a box around P.  This is shown schematically in Figure 4.10 for a Gaussian vortex.  The vorticity peak locations with !" > 0.9 were retained.    Figure 4.10.  Schematic of %2 variables for a Gaussian vortex.   94 The centroid was then calculated using equation (4.2), also defined in Chapter 2.  The choice of centroid accommodates the changing shape of the vortices as they convect downstream (Bohl, 2002) and would correspond to the location of highest vorticity if the vorticity had a Gaussian distribution.    56$7/57/89:6$7/:7/ (4.2) For computing the centroid, a spatial radius must be specified.  Since the vortex size varies significantly between the high and low Reynolds number (see Figure 4.9), a fixed radius was not used.  Instead, vorticity profiles in both the x- and y-directions were fit independently with a one-dimensional Gaussian distribution to calculate the Gaussian radius, rG, defined as the 1/e point of the vorticity distribution.  The average was then taken from the fit to provide the average Gaussian radius, rGavg, to provide an estimation of vortex core size.  In order to encircle the vortex core and minimize the contribution from the connecting braids, a spatial radius of 2.5rGavg was selected and only vorticity within the radius was included in the calculation of the centroid.  For the example of k = 6, the spatial radius is shown in Figure 4.11.    95  Figure 4.11.  Instantaneous spanwise vorticity field ("zc/U#) for k = 6.0, $o = 2¡, and different Reynolds numbers, showing the spatial radius (increasing from the top down). The centroid trajectories for the four Reynolds numbers is shown in Figure 4.12.  In each case, the positive-sign vortex originates above the wake centerline while the negative-sign vortex originates below it.  At Re = 2,000, the positive and negative vortices cross in the transverse direction at x/c ! 0.07 such that the negative vortex is above the positive vortex producing a traditional von K⁄rm⁄n street.  The vortices move away from the centerline and reach a nearly steady-state value at x/c ! 0.89.  At Re = 7,000, the vortices crossed the wake centerline at x/c ! 0.06 in a similar way as Re = 2,000.  The vortices initially move away from the centerline, reach a peak, then move toward the centerline and reached a nearly state-state value at x/c ! 1.5.    96   Figure 4.12.  Positive and negative vortex trajectories for different Reynolds numbers at k = 6.0 and $o = 2¡.   At Re = 12,000, the vortices cross the centerline twice, once at x/c ! 0.085 and again at x/c ! 0.61, with a peak in the trajectory occurring at x/c ! 0.35.  The vortices then form a reverse von K⁄rm⁄n street and reach a nearly steady-state at x/c ! 1.5.  More noticeable oscillations in the trajectory are visible over nearly the entire field of view.  At Re = 22,000, the vortices initially move towards the wake centerline and reached a minimum distance from the centerline at x/c ! 0.2.  The vortices then move away from the centerline and form a reverse von K⁄rm⁄n street.  The vortex trajectory shows no initial formation peak nor a steady-state condition at this Re even up to four chord lengths from the airfoil trailing edge (not shown).  The vortices move away from the centerline, reach a maximum, and then gradually move toward the centerline.   97 Oscillations are also observed in the trajectories for the two higher Reynolds numbers, and are less noticeable in the lower Reynolds numbers.  In order to quantify the wake structure evolution, the transverse and streamwise spacings of the vortices are computed for each Reynolds number.  The transverse spacing the vertical distance between the positive vortex and negative vortex from the wake centerline at each streamwise location (i.e., b = yp Ð yn).  The streamwise spacing was calculated using the relationship a = UC/f, where UC is the vortex convection speed and f is the pitching frequency.  Obtaining the a spacing from the convection speed allowed the evolution of a to be tracked throughout the entire domain and differs from how it was done by Bohl (2002) and Bohl and Koochesfahani (2009).  The convection speed of the vortex was computed by smoothing the centroid x-position with a moving window filter, and then numerically differentiating it with a second-order central difference scheme.  The spatial evolution of b and a are shown in Figure 4.13.  At Re = 2,000, the transverse spacing is negative, indicating the traditional von K⁄rm⁄n street arrangement and a velocity-deficit.  As the Reynolds number increases at this reduced frequency, the transverse spacing spacing decreases in magnitude towards zero, indicating an aligned vortex array.  The b spacing then switches sign and increases, producing reverse von K⁄rm⁄n street arrangement and jet-like behavior in the wake.  The streamwise spacing increases with increasing Reynolds number.  For each Reynolds number, the a spacing experiences a peak, which moves towards the trailing edge as the Reynolds number increases.  The streamwise spacing then decreases to a nearly asymptotic value for each Reynolds number.  At Re = 12,000 and 22,000, the streamwise spacing also shows an oscillatory signature.  98  Figure 4.13.  Wake vortex transverse spacing (top) and streamwise spacing (bottom) at k = 6.0 and $o = 2¡ for different Reynolds numbers.   The transverse and streamwise spacing of the vortex centroid for the range of reduced frequencies at the different Reynolds numbers are shown in Figure 4.14.  A dashed line is included in the b spacing plot to give a reference for the aligned vortex array condition (i.e., b = 0).  The b and a spacings were measured at x/c = 1.0.  The selection of this measurement location was based on three reasons: 1). this measurement location was downstream of the initial acceleration region before the vortex reached steady-state (refer to Figure 4.12 for example), 2) one chord length downstream of the trailing edge is a common location to measure velocity profiles for estimating thrust using a control volume analysis (Koochesfahani, 1989; Bohl and Koochesfahani, 2009; Naguib, et al., 2011), and 3) after the initial acceleration, the spacing did not significantly vary as the vortices convect downstream.  The Reynolds number dependence on  99 the peak vorticity, vortex circulation, core radius, and convection speed are provided in the Appendix G.   For each Reynolds number, the b spacing begins with negative values at low k and switches to positive values as the reduced frequency increases.  The b spacing for each Reynolds number reaches nearly the same asymptotic value of approximately 0.060c.  For reference, the trailing edge amplitude is 0.052c.  However, the reduced frequency at which the switch in orientation takes place depends on Reynolds number.  This crossover reduced frequency is discussed in more detail in the next section.  The a spacing decreases as the reduced frequency increases for each Reynolds number while it increases as the Reynolds number increases. Physically, this is a result of the vortices having a higher convection speed at the higher Re.  Figure 4.14.  Wake vortex transverse spacing (top) and streamwise spacing (bottom) dependency on k and Re for $o = 2¡.  Measurements were done at x/c = 1.0.     100 4.3.  Transverse Spacing Crossover Frequency  From the data in Figure 4.14, the effect of Reynolds number on the crossover reduced frequency at which the b spacing switches sign can be quantified.  The crossover reduced frequency dependence on Reynolds number is presented in Figure 4.15.  Data from the literature (Bohl and Koochesfahani, 2009; Naguib, et al., 2011; Monnier, et al., 2015), as well as the drag-to-thrust crossover reduced frequency, are also included in the figure for reference.  In addition to measurements obtained at x/c = 1.0, the crossover reduced frequency was also determined from visual inspection of the flow field where the vortex array was nearly aligned (see Figure 4.16).  As the Reynolds number increases from 2,000 to 22,000, the crossover reduced frequency decreases from kcr ! 9 to kcr ! 5.2.  The decrease in kcr according to Re can be represented by a least squares power-law fit with an exponent of -0.23.  The data also compare quite well with the data found in the literature.  Noteworthy is that the conclusion by Bohl (2002) that the switch in wake orientation occurs before the switch from drag-to-thrust holds for a wide range of Re.  Interestingly, the ratio of the thrust crossover reduced frequency and orientation crossover reduced frequency does not stay at a fixed value; it decreases as Re increases such that thrust occurs at lower kÕs relative to kcr for b = 0 for higher Reynolds compared to lower Reynolds numbers.  101  Figure 4.15.  Crossover reduced frequency for b = 0 dependence on Reynolds number for $o = 2¡.  Data from literature are included for reference.  Note that Monnier, et al., (2015) used a modified NACA series airfoil.  Figure 4.16.  Instantaneous spanwise vorticity field ("zc/U#) for visually aligned vortex array and $o = 2¡ for different Reynolds numbers (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (!/T = 0.0).    102 4.4.  Effect of Pitching Amplitude The effect of pitching amplitude on the flow structure for different Reynolds numbers is shown in Figure 4.17 for k = 3.  The reduced frequency was selected to highlight the formation of a traditional von K⁄rm⁄n street at $o " 4¡ for the Re = 2,000 and 7,000 that did not occur at Re = 12,000 and 22,000.  Flow visualization for additional kÕs is presented in Appendix H.  The airfoil is at zero angle of attack and is pitching up (!/T = 0.0).  At Re = 2,000 for $o = 2¡, the vortex configuration is a traditional von K⁄rm⁄n street.  As the amplitude increases, the vortices become aligned, then switch to a reverse von K⁄rm⁄n street.  At Re = 7,000 for $o = 2¡, the vortex configuration is similar to traditional von K⁄rm⁄n street, but with slightly different features.  As the amplitude increases, the vortex pattern becomes a traditional von K⁄rm⁄n street at the $o = 4¡.  Increasing the amplitude to $o = 6¡, a second vortex sheds each half-cycle that pair downstream to form a reverse von K⁄rm⁄n street.  The wake structure continues forming into a reverse von K⁄rm⁄n street as the amplitude increases further.  As the Reynolds number increases to 12,000 and 22,000, more complex structures are produced in the wake at this k at $o = 2¡-4¡.  As the amplitude increases to $o = 6¡ and 8¡, the wake structure becomes reverse von K⁄rm⁄n-like, with small additional vortices shedding and pairing with the main vortex.  At $o = 10¡, Re = 12,000 produces a reverse von K⁄rm⁄n street and Re = 22,000 generates a von K⁄rm⁄n-like pattern.  Note that a traditional von K⁄rm⁄n street was not observed for Re = 12,000 or 22,000 when $o " 4¡ at this reduced frequency or the others considered in this work.   103  Figure 4.17.  Instantaneous spanwise vorticity field (!zc/U") at k = 3.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (#/T = 0.0).  The field of view extends over the range 0 ! x/c ! 3.0 and -0.4 ! y/c ! 0.4.     104  In order to quantify the effect of pitching amplitude on the wake vortex arrangement, the b and a spacing are again computed.  For this analysis, only the cases in which a single pair of opposite sign vortices shed per cycle were considered.  The transverse and streamwise spacing dependence on pitching amplitudes for different reduced frequencies at Re = 2,000 and 22,000 are shown in Figure 4.18.  As the amplitude increases for each Reynolds number, the transverse spacing increases at a given k.  There is a slight increase in streamwise spacing as the amplitude increases at a given k for each Reynolds number.  The streamwise spacing shows more sensitivity to reduced frequency than either Re or !o, as expected since a is a first order approximation of k.  By rescaling the spacing data with the Strouhal number (see Figure 4.19), the b spacing shows slightly better collapse while the a spacing data forms into distinctive curves separated by amplitude.   Figure 4.18.  Wake vortex transverse (top) and streamwise (bottom) spacing dependence on k, for Re = 2,000 (left) and 22,000 (right) and !o = 2¡-10¡.  Data is presented for only cases with well-formed single vortex pairs.  105  Figure 4.19.  Wake vortex transverse (top) and streamwise (bottom) spacing dependence on k, for Re = 2,000 (left) and 22,000 (right) and !o = 2¡-10¡.  Data is presented for only cases with well-formed single vortex pairs. Using the transverse spacing data from Figures 4.18 and 4.19 the vortex orientation crossover reduced frequency/Strouhal number dependence on Reynolds number for different amplitudes may be quantified.  These crossover data are presented in Figure 4.20 based on a measurement location at x/c = 1.0 and again only for cases in which a single pair of alternating vortices shed each oscillation cycle.  At 4¡ amplitude, the vortex pattern shows a switch in b only for Re = 2,000 and 7,000.  The crossover reduced frequency is lower at 4¡ than 2¡, as expected based on the drag-to-thrust crossover data presented in Chapter 4.  At Re = 12,000 and 22,000, there is no discernable transition from a traditional von K⁄rm⁄n to reverse von K⁄rm⁄n street, causing no crossover to occur (refer to Appendix H).  As the amplitude increases further, only the Re = 2,000 case continues to produce a clear switch in b.  When the crossover reduced frequency is rescaled in terms of a crossover Strouhal number, there is better data collapse such that the crossover Strouhal number decreases as the Reynolds number increases and slightly  106 increases as the amplitude increases.  For Re = 2,000, the crossover Strouhal number range is 0.15-0.17 over the amplitude range of 2¡-10¡ while it is 0.11-0.125 for 2¡-4¡ amplitude at Re = 7,000.  There is also good agreement with the literature.    Figure 4.20.  Wake vortex transverse spacing switch crossover reduced frequency (left) and Strouhal number (right) dependence on Reynolds number for different amplitudes.  Measurements from computations done at x/c = 1.0.  Data from literature are included for reference.  Note that Godoy-Diana, et al., (2008) used for a tear-drop shaped airfoil and Monnier, et al., (2015) used a modified NACA series airfoil. Finally, the vortex configuration may be quantified as a single parameter, the vortex aspect ratio (b/a).   Bohl and Koochesfahani (2009) showed that once the transverse spacing reached an asymptotic value, the aspect ratio continued to exhibit a slight increase as the reduced frequency increased further for Re = 12,600 and !o = 2¡.  At the highest reduced frequency, the aspect ratio reached a value of 0.19 in that study, which was much lower than value of 0.28 for a stable point vortex street (von K⁄rm⁄n and Burgers, 1935).  The vortex aspect ratio from the present computations for each Re and different amplitudes as a function of Strouhal number is shown in Figure 4.21.  It is notable that the the aspect ratio nearly collapses with respect to amplitude for each Reynolds number.  As the Strouhal number increases, the aspect ratio  107 decreases in magnitude, switches sign, then increases in magnitude for each Reynolds number.  At Re = 2,000, this increase is a gentle.  As the Reynolds number increases, the increase in aspect ratio is much ÒsharperÓ as Strouhal number increases at the low Strouhal number range (see Figure 4.22).    Figure 4.21.  Wake vortex configuration aspect ratio (b/a) dependence on St for different Reynolds number and pitching amplitudes.  Measurements from computations done at x/c = 1.0.   108  Figure 4.22.  Wake vortex configuration aspect ratio (b/a) dependence on St for different Reynolds number for !o = 2¡ (circles) and 10¡ (squares).  Measurements from computations done at x/c = 1.0.     109 CHAPTER 5.  EFFECT OF TRAJECTORY ASYMMETRY ON AERODYNAMIC LOADS  In addition to studying the effect of Reynolds number on the average and fluctuating loads on a pitching airfoil, the effect of trajectory asymmetry is also of interest based on previous work (Koochesfahani, 1989).  According to Koochesfahani (1989), Òit was expected that the load history would be strongly affectedÓ due to symmetry breaking in the wake.  However, forces were not quantified in that investigation.   The Reynolds number and pitching amplitude were selected to be 12,000 and 2¡, respectively, to be in agreement with the previous experiment by Koochesfahani (1989).  The baseline reduced frequency was 6.68 to closely match the experiment by Koochesfahani (1989), though the reduced frequency is varied from 4 to 8 in this work to illustrate the effect of k.  The lowest reduced frequency considered corresponds to the lowest k considered by Bohl (2002) in which flow visualization was performed.  The highest reduced frequency was selected to be close to the drag-to-thrust crossover reduced frequency from Chapter 3.  Reduced frequencies greater than 8 were not pursued since Mach number effects become more prevalent for the sinusoidal pitching cases at k ! 10 and the addition of asymmetry would make these effects more challenging to overcome.  The freestream Mach number is 0.005 for all data presented in this chapter.  As stated in Chapter 2, a small amount of smoothing corresponding to 0.5%TPU was applied to the angle of attack history.  The value of S is less than 50% such that the pitch-up is faster than the pitch-down.  Example angle of attack, velocity, and acceleration time histories are provided schematically for different values of S in Figure 5.1 for reference.     110  Figure 5.1.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) for different values of S.   This chapter is organized as follows.  Section 5.1 deals with the effect of asymmetry on the average loads, Section 5.2 provides a discussion of decomposition of the lift into terms analogous to circulatory and non-circulatory components of inviscid linear theory by Theodorsen (1935), and Section 5.3 shows the load fluctuation caused by the asymmetrically pitching airfoil.    111 5.1.  Average Forces First, the influence of reduced frequency on the average force coefficients for the asymmetrically pitching airfoil is discussed.  The time-average lift and thrust coefficients are presented in Figure 5.2.  Positive lift corresponds to the positive Y direction and positive thrust denotes a forward propelling force in the negative X direction (refer to the schematic in Figure 1.11 in Chapter 1).  When the pitching motion becomes asymmetric, the lift becomes non-zero. Sarkar and Venkatraman (2006) also observed non-zero average lift, stating that that the average lift increased with increasing asymmetry.  However, quantitative data were not presented by the authors.   In the present computations, the average lift depends on both reduced frequency and asymmetry.  It is observed that over the range of investigated asymmetries for these flow conditions, the sign of the lift is primarily influenced by the reduced frequency while the variation of the lift with respect to S also depends on k.  At S = 50% (sinusoidal), the average lift is zero for all reduced frequencies, as expected.  At k = 4, the lift increases as the asymmetry increases until it reaches a nearly asymptotic value at S = 38%.  As the reduced frequency increases from 4 to 5.2, there is a noticeable decrease in average lift over the entire investigated S range.  There is also an interesting dip in lift at S = 38% for this reduced frequency, but the cause of this feature is currently unclear.  As the reduced frequency increases to 5.8, the net lift is nearly zero over the entire range of S.  As the reduced frequency increases to 6.68, the lift becomes negative.  As asymmetry increases, the negative lift increases in magnitude until it reaches a maximum at S " 34%, then begins decreasing.  Further increase in reduced frequency to 8 produces negative lift with greater magnitude than at k = 6.68.  The value of S in which maximum negative occurs moves towards 50% as k increases from 6.68 to 8.  112  Figure 5.2.  CL (left) and CT (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.   The effect on thrust is also interesting.  For sinusoidal pitching, the airfoil experiences drag at each reduced frequency.  Recall from Chapter 3 that the crossover reduced frequency for this Reynolds number and pitching amplitude was approximately 8.1.  At a given reduced frequency, CT increases away from the value at S = 50%.  At k # 6.68, the value of CT increases towards zero (i.e., towards positive thrust).  At k = 8, the force switches from drag to thrust at S = 40%.  The increase in CT makes intuitive sense since the Òeffective frequencyÓ of the airfoil increases when pitching asymmetrically.  An additional feature of the thrust is that its differential between S = 50% and S $ 50% is not constant, but it depends on reduced frequency, as shown in Figure 5.3.   It is found that the thrust differential increases non-linearly with both increasing reduced frequency and asymmetry.    113  Figure 5.3.  Thrust differential vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.   Decomposition of the average lift coefficient into its corresponding components caused by the average pressure and shear stress show that both components experience similar variation by k and S seen in the total lift (see Figure 5.4), but the component from the shear stress is two orders of magnitude smaller than the component from the pressure.  Performing the same exercise on the average thrust coefficient shows that pressure contribution experiences drag at k = 4 and thrust at k ! 5.2.  The shear stress contribution produces slightly increasing drag as the reduced frequency and asymmetry increase.   114  Figure 5.4.  Pressure (top) and shear stress (bottom) components of average CL (left) and CT (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.   In order to better understand the cause of the non-zero lift, the contributions from the average positive lift, !"#, and negative lift, !"$, are isolated.  Figure 5.5 provides a schematic of the contributions of positive and negative lift and their respective temporal durations with respect to the oscillation cycle, %"# and %"$.  The angle of attack and pitching acceleration time histories are included in the figure for reference.  The net contribution from the positive and negative lift to the average total lift (seen in Figure 5.2) may be calculated using the cycle average, given equations (5.1) and (5.2).  However, averaging over the oscillation cycle, T, does not provide information about the amount of ÒareaÓ integrated under the positive and negative lift curves.  The amount of positive and negative lift during !"# and !"$ may be quantified using a equations (5.3) and (5.4), respectively, in which the averaging for the positive and negative lift components is done over %"# and %"$, respectively, instead of the oscillation period, T.  These two expressions are called the duration averages and are related to equations (5.1) and (5.2) by equations (5.5)  115 and (5.6).  Note that the lift pre-dominantly is in phase with the pitching acceleration in Figure 5.5.    Figure 5.5.  Schematic illustrating lift averaging definitions for the example of k = 6.68 and S = 38%.  The angle of attack and pitching acceleration are included for reference.  !"#&'%!"#(#)(*+* (5.1)  !"$&'%!"$(#)(*+* (5.2)  !"#&'%"#!"#(#),-(*+* (5.3)  116  !"$&'%"$!"$(#),.(*+* (5.4) where  %"#/%"$&% (5.5)  !"&!"#/!"$&%"#%!"#/%"$%!"$ (5.6) The results of this averaging analysis are shown in Figure 5.6.  In the figure, the ratio of the positive lift and negative lift is taken to illustrate relative differences while the absolute value is shown to display magnitude.  For sinusoidal pitching (S = 50%), these two components are equal and the ratio would be 1.  Using the oscillation cycle time-average of the positive and negative components (i.e., equations (5.1) and (5.2)), similar behavior to that presented in Figure 5.2 is observed.  At k < 5.8, the average positive lift is greater than the negative lift for S < 50%.  When k = 5.8, the positive lift nearly balances the negative lift over the entire range of S investigated here.  At k > 5.8, the negative lift exceeds the positive lift over the investigated range of S.  The ratio of positive to negative lift duration shows that as the asymmetry increases, the duration of negative lift increases for all reduced frequencies.  At a given value of S, negative lift occurs over a larger duration than the positive lift as k increases.  The ratio of the duration average of positive to negative lift shows that more positive lift occurs over the shorter %"# than the negative lift over the longer %"$at each value of S < 50% regardless of the reduced frequency.  Both the duration ratio and duration average ratio show different variation with S compared to the cycle average ratio, making it difficult to ascertain what mechanism causes the trend in lift seen in Figure 5.2.  In order to elucidate the underlying mechanism causing the behavior of the  117 average lift with respect to both S and k, efforts were undertaken to decompose the lift into contributions analogous the circulatory and non-circulatory components of lift from inviscid linear theory, as discussed in the next section.    Figure 5.6.  Ratio of cycle time-average of positive to negative lift (left), ratio of the duration of positive to that of negative lift period (middle), and ratio of duration time-average of positive to negative lift (right). 5.2.  Lift Decomposition Classical theory by Theodorsen (1935) showed that the unsteady lift on a flat plate undergoing harmonic pitching and plunging in an incompressible, inviscid flow can be decomposed into two contributions.  The first source of lift is that associated with the development of circulation around the airfoil, Lc.  In order to incorporate this shedding of vorticity and its influence on the plate, Theodorsen introduced his complex lift attenuation function, C(k) (see Theodorsen (1935) for more details).  The second source of lift is the non-circulatory (or added mass) component, Lnc, which originates from the reaction force on the oscillating plate caused by the displacement of the mass of fluid by the plate.  For high reduced  118 frequency oscillations, the non-circulatory component contributes the majority of the total time-dependent lift coefficient.  For example, consider the ratio of the circulatory and non-circulatory lift coefficients to the total lift at the phase of maximum total lift, as shown in Figure 5.7.  At k < 1.5, the circulatory component provides the majority of the total lift when the total lift is at its maximum.  At k > 1.5, the non-circulatory component provides the majority of the total lift when the total lift is at its maximum.  However, the applicability of TheodorsenÕs lift formula to an asymmetrically pitching airfoil is unclear since the theory was derived for angle of attack variation of the form 0*&01234(.    Figure 5.7.  Ratio of lift circulatory and non-circulatory components to total lift at the phase of maximum lift.  The lift coefficient time history for k = 6.68 (indicated by the blue line) using TheodorsenÕs formula is included in the figure for reference, where the phase of maximum lift is indicated by the green dashed line. More recently, Wang, et al., (2013) introduced a simple formula for calculating the force on a body using only velocity data for incompressible flow at low Reynolds numbers.  Consider a solid body, B, that lies in a control volume, Vf, with control surface ".    119  Figure 5.8.  Schematic of fluid domain, Vf, with control surface " around solid body B. The force acting on the body, 5, is given by  5&6678/9:+;&6<=>+*/>?@>AB+C/678/9D+; (5.7) Where p is the pressure, 8 is the normal vector, 9 is the surface shear stress vector, > is the velocity vector.  By substituting >?@>&EF>/@>GHI into equation (5.7), the force is calculated as  5&<>FEAB+C6<=>+*AB+C/67/<>GI8D+; /9D+;6<>GHI:8+; (5.8) where E is the vorticity vector,  >FE is called the Lamb vector, and > is the velocity magnitude.  According to Wang, et al., (2013), the first term in equation (5.8) is the Lamb vector integral and represents the vortex force (Prandtl, 1918; Saffman, 1992; Wang, et al., 2013), the  120 second term is the unsteady inertial effect, the third and fourth terms are the contribution from the total pressure and shear stress at the control surface ", and the fifth term is the boundary term that depends on the surface velocity of the body.  In an inviscid, irrotational flow, the combination of the second and fifth term is the added mass (Wang, et al., 2013).  By selecting a rectangular control surface that is far from the solid body, the third and fourth terms are negligible.  In the special case of two-dimensional steady flow and assuming >Jof the form >&KL/>M (where >M is the velocity induced by vorticity), equation (5.8) reduces to classical Kutta-Joukowski theorem N&<KLO (Prandtl, 1918; von Karman and Burgers, 1935; Wang, et al., 2013).   For an unsteady thin airfoil, the time-dependent two-dimensional lift force, L(t), reduces to two leading order terms, given by equation (5.9).  In this expression, "# is freestream density, u is the local streamwise component of velocity, v is the local transverse component of velocity, $z is the local spanwise vorticity, t is the physical time, and A is the two-dimensional area of the control volume that extends far from the body.    N*&6<L>EP+QRSTUVWXJYZ[VJ\WU]/<L=^=*+QR_``WaWUbVZTcJ\WU] (5.9) As seen in equation (5.9), the unsteady lift comes from two primary sources.  The first term is the vortex lift and the second term is the lift from the fluid acceleration.  These two components of lift are analogous to the circulatory and non-circulatory components in the inviscid theory by Theodorsen (1935).  Wang, et al., (2013) computationally showed that the equation (5.9) is very effective in accurately estimating both the time-average and unsteady lift for a variety of flows around stationary and flapping flat plates at Re ~ O(102).  This  121 decomposition is exploited here to determine the underlying contributions to the average lift.  In this work, only the vortex lift coefficient, CVL, was directly computed from the instantaneous velocity and vorticity fields due to difficulty in the ability to accurately compute =^H=* using the phases extracted from the computations.  Since only the vortex lift is calculated, the acceleration component of lift is considered the remainder between the total lift and the vortex lift. To calculate the vortex lift coefficient, a rectangular domain was created around the airfoil that extended H ahead of the leading edge, H/2 above and below the symmetry line, and ends at the trailing edge.  This is shown schematically in Figure 5.9.  For this analysis, H = 50c.  The vortex lift was then calculated according to N*&6>EPdQR, where A is the control volume area and grid cell area, dQ, approximated by the inverse of the transformation Jacobian at each grid point (i.e., dQ&e$f&ghij6ihgj).  This procedure was performed for each phase to produce a time history of the vortex lift coefficient.    Figure 5.9.  Schematic of integration domain for computing the vortex lift. An example vortex lift coefficient time history is shown in Figure 5.10 for k = 6.68 and S = 38% as a representative case.  The pitching acceleration and total lift histories are included for reference.  Two features of the vortex lift coefficient are evident from the figure.  The first is that the vortex lift is approximately in phase with the angle of attack as opposed to the total lift, which is approximately in phase with the acceleration.  The second is that the amplitude of the  122 lift is much lower magnitude than the total lift.  The peak-to-peak amplitude of the vortex lift is approximately 2% of the peak-to-peak amplitude of the total lift.    Figure 5.10.  Time history of angle of attack (top left), acceleration (top right), vortex lift coefficient (bottom left), and total lift coefficient (bottom right) for k = 6.68 and S = 38%. Performing the time-average of the unsteady vortex lift for each reduced frequency and asymmetry and plotting the actual average lift, CL, versus the average vortex lift produces a linear relationship with unit slope (see Figure 5.11).  The R2 value of the entire data set was 0.99996, indicating a very strong correlation.  In fact, the average vortex lift coefficient values were within 2% of the actual average lift coefficient.  It can be inferred that the majority of the non-zero time-average lift is a result of the vortex lift component.    123  Figure 5.11.  CL vs. CVL for each reduced frequency and asymmetry.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  This behavior was not expected, though not surprising.  Consideration of the investigation by Wang, et al., (2013) sheds light on this observation.  As stated previously, the authors performed computations of a flat plate undergoing sinusoidal pitching and plunging about a 10¡ angle of attack at a Reynolds number of 300.  The pitching amplitude was 30¡ while the plunging amplitude was c/4.  The authors observed a very non-sinusoidal signature in the total lift, as shown in Figure 5.12 as the solid black line.  The contribution from the vortex lift (i.e., integrated Lamb vector; dashed line) is nearly sinusoidal and while the acceleration contribution (dot-dashed line) shows similar fluctuations as the total lift.  In fact, the authors state that Òthe Lamb vector integral mainly contributes to the mean lift with a considerable phase shift compared to DNS, while the acceleration term significantly modifies the waveform and phase.  The sum of the Lamb vector integral and the acceleration term recovers the true waveform.Ó  Both the vortex lift results by Wang, et al., (2013) and from the current computations have nearly sinusoidal waveforms.    124  Figure 5.12.  The time histories of lift, integrated Lamb vector contribution, and acceleration from computations by Wang, et al., (2013) for a flat plate at non-zero angle of attack undergoing combined pitching and plunging at Re = 300.  Instantaneous fields of spanwise vorticity (top), Lamb vector projected in the vertical direction (middle), and acceleration projected in the vertical direction (bottom) at % = 18.8 (indicated by the red dashed line) are included on the right for reference. Figure courtesy of AIP Publishing, LLC. The vortex lift shown in Figure 5.12 visually does have a peak-to-peak amplitude an order of magnitude greater than in the current computations.   However, the flow considered by the Wang, et al., (2013) shows boundary layer asymmetry and separation due to the non-zero mean angle of attack and large oscillation amplitudes whereas the mean angle of attack is zero and pitching amplitude is small in the current work, and may lead to more subtle asymmetries in the boundary layer.  Regardless, differences in waveform and peak-to-peak amplitude between the vortex lift and the total lift from the present computations should not be considered unexpected based on this discussion. The effect of reduced frequency on the unsteady vortex lift is presented in Figure 5.13, with the angle of attack included for reference.  The reduced frequencies in the figure are 4 and 8 to represent the lowest and highest kÕs considered.  As previously seen in Figure 5.10, the vortex  125 lift is approximately in phase with the angle of attack for each asymmetry and reduced frequency.  In addition, both reduced frequencies show that the maximum vortex lift increases as the asymmetry increases.  However, the two reduced frequencies show noticeable differences.  As the asymmetry increases for k = 4, the minimum vortex lift decreases in magnitude while an increasing amount of the vortex lift becomes positive.  For k = 8, the minimum vortex lift decreases slightly in magnitude while an increasing amount of the vortex lift becomes negative as asymmetry increases.   Figure 5.13.  Angle of attack (top) and vortex lift coefficient (bottom) time histories for k = 4 (left) and 8 (right). A similar averaging procedure done for average total lift was performed for the vortex lift (shown schematically in Figure 5.14) using equations (5.11) - (5.16) to quantify the contribution of positive and negative sign vortex lift.  126  Figure 5.14.  Schematic illustrating vortex lift averaging definitions for the example of k = 6.68 and S = 38%.  The angle of attack is included for reference.  !A"#&'%!A"#(#)(*+* (5.11)  !A"$&'%!A"$(#)(*+* (5.12)  !A"#&'%A"#!A"#(#)k,-(*+* (5.13)  !A"$&'%A"$!A"$(#)k,.(*+* (5.14) where  %A"#/%A"$&% (5.15)  127  !A"&!A"#/!A"$&%A"#%!A"#/%A"$%!A"$ (5.16) The ratio of positive and negative vortex lift components using the cycle time-average showed similar behavior observed in Figure 5.2 and thus is not included in this analysis.  The ratios of the positive vortex lift duration versus negative vortex lift duration and duration average are shown in Figure 5.15.  These data distinctly show similar trends as seen in the cycle time-average of the total lift in Figure 5.2.  At k = 4 for each asymmetry, the duration of the positive vortex lift exceeds the duration of negative vortex lift (i.e., %A"#/%A"$ > 1), as seen visually in Figure 5.14.  As the reduced frequency increases at a given value of S, an increasing duration of the vortex lift is negative until %A"#/%A"$ < 1.  The duration average ratio, !A"#H!A"$, shows the amount of positive vortex lift decreases relative to the negative vortex lift as the reduced frequency increases at a given value of S.  It can be inferred that at k = 4, more positive vortex lift occurs over a large duration compared to the negative vortex lift.  As the reduced frequency increases, an increasing amount of negative vortex lift occurs over an increasing duration compared to the positive vortex lift component.     128  Figure 5.15.  Ratio of the duration of positive to that of negative lift period (left), and ratio of duration time-average of positive to negative lift (right). This can be clearly seen in Figure 5.16, which shows the unsteady vortex lift for both reduced frequencies at S = 38%, with the duration of the positive and negative lifts highlighted schematically.  There is clearly more positive vortex lift visually than negative vortex lift, both in duration and amount, for k = 4 while the opposite is true for k = 8.    129  Figure 5.16.  Vortex lift coefficient time histories for k = 4 (left) and 8 (right), with S = 38%. We note that it has been reported in the literature that natural flyers do flap asymmetrically (Park, et al., 2001; Ros”n, et al., 2004).  Both authors showed that the downstroke fraction of the wingbeat cycle period (i.e., S) of barn swallows and thrush nightingales decreased as the flight speed increased.  This behavior implies an increase in Reynolds number and a decrease in reduced frequency.  The value of S dropped from around 50% at low speed (U " 10m/s) to approximately 40%-45% at higher speeds (U > 10m/s).  This is consistent with the current data in the following ways.  There was positive lift at low k when pitching asymmetrically with S < 50% in the current computations.  It can be inferred that as the reduced frequency decreases to the reduced frequencies considered by both authors (O(10-1)), the lift becomes higher.  Although these results are interesting, the value of average lift values produced by these computations undergoing idealized pitching are small (CL ~ O(10-2)).  Realistic flapping is much more complex and presumably leads to higher lift and thrust coefficients.    130 5.3.  Force Fluctuation For completion, consideration is made on the effect of airfoil trajectory asymmetry to the force fluctuation.  The lift and drag fluctuations, !"M and !lM, for three values of S are shown in Figure 5.17 for a representative reduced frequency k = 6.68.  The angle of attack and pitching acceleration histories for the three cases are included for reference.  The lift and drag time histories show a very strong dependency on S.  For S = 50%, both the lift and drag fluctuation vary sinusoidally with time.  As the asymmetry increases, new features can be seen in the forces.  First, the force fluctuations become non-sinusoidal and follow the pitching acceleration, suggesting that the force fluctuations are primarily a result of non-circulatory forces.  Second, the force fluctuations exhibit an inertial spike just after the jump in acceleration.  Third, the maximum and minimum force magnitudes increase as S increases.   131  Figure 5.17.  Lift and drag fluctuation for different asymmetries at k = 6.68.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.  The angle of attack and pitching acceleration history is included for reference.  The effect of both asymmetry and reduced frequency are quantified by the maximum, minimum, and peak-to-peak amplitude of lift and drag coefficient fluctuation, as shown in Figures 5.18.  As the reduced frequency and asymmetry increases, the lift and drag max/min and peak-to-peak values increase as well.  The additional lift and drag due asymmetry varies non-linearly with S.    132  Figure 5.18.  Maximum (top), minimum (middle), and peak-to-peak (bottom) values of CL (left) and CD (right) vs. S for different reduced frequencies.  The Reynolds number is 12,000 and the pitching amplitude is 2¡.    133 CHAPTER 6.  EFFECT OF TRAJECTORY ASYMMETRY ON WAKE STRUCTURE  The final results that will be discussed relate to the effect of the trajectory asymmetry on qualitative and quantitative characteristics of the wake structure.  As shown by Koochesfahani (1989), the extent of asymmetry can have a profound impact on the number of vortices shed and their orientation due to the vortex interactions in the wake.  For the purposes of this chapter, a limited range of reduced frequencies were considered.  The reduced frequencies are 5.2, 5.8, and 6.68 to show the special cases of a traditional von K⁄rm⁄n street, neutral wake, and reverse von K⁄rm⁄n street when S = 50%.  The Reynolds number is 12,000 and freestream Mach number is 0.005 for all data presented in this chapter.   This chapter is organized as follows.  Section 6.1 contains an examination of the effect of trajectory asymmetry on the wake structure qualitatively Section 6.2 presents the basic wake provides a discussion of the vortex formation at the trailing edge, and Section 6.3 gives a presentation of the effect of trajectory asymmetry on the vortex properties. 6.1. Wake Structure  Qualitative features of the wake structure for an asymmetric trajectory are discussed for different reduced frequencies.  The vortical structures in the wake of the airfoil for the three reduced frequencies and each asymmetric are shown in Figure 6.1.  The wake structures for k = 4 and 8 are provided in Appendix I.  Approximately two chord lengths in the wake is shown in the figure.  In the figure, reduced frequency increases from left to right while asymmetry increases from the top down.  The airfoil is at zero angle of attack and is pitching up (%/T = 0.0) for each instantaneous field shown.  In each case presented, a single positive vortex forms during the  134 pitch-up.  One or more negative vortices shed during the slow pitch-down, where the number of negative vortices depends on reduced frequency and asymmetry, as discussed next.    At low asymmetries (S = 46%-42%), there is a very subtle effect caused by asymmetry for each reduced frequency.  The vortex structure has a slight downward deflection by observing both the positive and negative vortices below the wake centerline (y/c = 0).  The deflection angle of the wake visually increases as the asymmetry increases from S = 46% to 42%.  At k = 5.2 and S = 42%, the negative-sign vorticity braid rolls into a vortex but quickly pairs with the larger negative vortex.  A secondary vortex sheds at S = 42% for k = 6.68 but also pairs with the larger negative vortex downstream.  No secondary vortex is observed for k = 5.8.  At S = 38%, each reduced frequency produces two negative vortices.  As the asymmetry increases further, k = 5.2 and 5.8 continue to produce two negative vortices during the pitch-down, with slightly different orientations.  At k = 6.68, three negative vortices shed into the wake. 135  Figure 6.1.  Instantaneous spanwise vorticity field (!zc/U") for Re = 12,000 and #o = 2¡ for different reduced frequencies (increasing to the right) and asymmetries (increasing from top down).  The airfoil is at zero angle of attack and is pitching up ($/T = 0.0).  The range covered in the figure is 0.5 ! X/c ! 3 and -0.3 ! Y/c ! 0.3. 136 Although asymmetry in the trajectory of the vortices can be seen visually in Figure 6.1, it is helpful to quantify the motion of the vortices by tracking the centroid through space.  Two values of S (50% and 38%) are considered for the three reduced frequencies in this exercise to demonstrate differences in how the vortices move when the airfoil trajectory is asymmetric.  The vortices considered are shown in Figure 6.2.    Figure 6.2.  Schematic of vortex labeling for both S = 50% and S = 38%. The trajectories of the vortex centroids are shown for the three reduced frequencies and two values of S in Figure 6.3.  As shown in the figure, N1 and N2 followed the same trajectory for each reduced frequency when S = 50% due to symmetry in the pitching motion.   When S = 38%, there is clear asymmetry in the vortex trajectories such that a steady state is not obtained over the downstream range considered in this work.  For each reduced frequency, both P1 and N2 convect below the wake centerline over the measurement domain whereas N1 convect above the centerline.  As the k increases, each vortex moves closer to the wake centerline, particularly downstream.  The highest reduced frequency also shows the centroid trajectories changing directions and moving towards the wake centerline between 2.7c and 3c.    137  Figure 6.3.  Positive and negative vortex centroid trajectories at three reduced frequencies numbers for S = 50% and 38%.  The Reynolds number is 12,000 and the pitching amplitude is 2¡. 138 6.2.  Vortex Formation at Trailing Edge In order to better understand the effect of trajectory asymmetry on the wake structure, it is of interest to study the vortex formation near the trailing edge that serves as the initial conditions for the subsequent evolution of the wake.  Figures 6.4-6.9 show the vorticity field near the trailing edge at twelve evenly-spaced oscillation phases for three asymmetries (S = 50%, 38%, and 30%) at a reduced frequency of 6.68.  The figures are divided by pitch-up and pitch-down phases, which are marked on the accompanying angle of attack, pitching velocity, and pitching acceleration.  An additional example for k = 5.2 at the same values of S is provided in Appendix I.  The three selected cases generally show the vorticity evolution associated with the formation of one, two, or three negative vortices to form in the wake.   First, the sinusoidal case will be examined as shown in Figure 6.4 for the pitch-up and Figure 6.5 for the pitch-down.  As the angle of attack begins increasing, a counter-clockwise rotating (positive) vortex forms at the trailing edge.  At !/T = 0.792, a small region of high magnitude negative vorticity is observed below the positive vortex.  As the airfoil continues pitching up, the upper surface boundary layer approaches the negative vorticity near the trailing edge.  By !/T = 0.875, the upper surface boundary layer and trailing edge negative vorticity meet.  The positive vortex subsequently detaches from the airfoil and sheds into the while continuing to be reinforced by the lower surface ÒfeedingÓ shear layer.  At !/T = 0.083-0.125, this ÒfeedingÓ shear layer becomes separated from the vortex and continues to shed into the wake, forming a vorticity braid.  A region of high magnitude negative vorticity that forms at the trailing edge is also observed at these phases.  At !/T = 0.208-0.250, this high magnitude negative vorticity spreads to cover more of the upper surface and some of the lower surface.    139 As the airfoil pitches down, a clockwise (negative) vortex forms at the trailing edge.  At !/T = 0.333, a region of high magnitude positive vorticity forms at the trailing edge.  This positive vorticity meets the lower surface boundary layer by !/T = 0.375 and the negative vortex detaches from the airfoil and sheds into the wake while being reinforced by the high vorticity ÒfeedingÓ shear layer.  The ÒfeedingÓ shear layer continuous shedding into the wake and forms a vorticity braid after being separated from the negative vortex.  As the airfoil reaches its minimum angle of attack, a region of high magnitude positive vorticity is observed to form near the trailing edge, and spreads to the upper and lower surface near the trailing edge.     140  Figure 6.4.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 6.68, S = 50%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 141  Figure 6.5.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 6.68, S = 50%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 142 When trajectory was asymmetric, the vortices formed asymmetrically.  For S = 46% (not shown), the effect of asymmetry on the vortex formation was subtle.  When S = 38%, there was a more obvious effect on the vortex formation (see Figure 6.6-6.7).  As the airfoil pitches up, positive vorticity is observed to form on the upper and lower surface near the trailing edge as in the case of S = 50%.  However, the center of the counter-clockwise (positive) vortex is more upstream that it is for S = 50%.  A slightly larger region of high negative vorticity compared to S = 50% is observed below the positive vortex.  This negative vorticity meets the upper surface boundary layer and the positive sign vortex sheds into the wake.   When the airfoil pitches down, a clockwise (negative) vortex forms at the trailing edge, with regions of negative vorticity on the upper and lower surface near the trailing edge.  A small region of high magnitude positive vorticity forms below the negative vortex and meets the lower surface boundary layer at !/T = 0.240.  The negative vortex then sheds into the wake.  As the airfoil continues pitching down at !/T = 0.344, a second negative vortex is observed near the trailing edge and is connected to the first shed negative vortex by a thin layer of negative vorticity.  At !/T = 0.5, the lower surface boundary layer connects with a small region of positive vorticity at the trailing edge while the second negative vortex sheds into the wake.  At !/T = 0.656-0.708, a third ÒconcentrationÓ of high magnitude negative vorticity is observed near the trailing edge in the feeding shear layer while high magnitude positive vorticity forms at the trailing edge.  The third ÒconcentrationÓ of negative vorticity sheds into the wake and becomes a braid downstream while high positive magnitude vorticity spreads along the upper and lower surface near the trailing edge.   143  Figure 6.6.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 6.68, S = 38%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 144  Figure 6.7.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 6.68, S = 38%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 145 For S = 30%, the positive vortex forms rapidly at the trailing edge during the fast pitch-up and sheds into the wake (Figure 6.8).  As the airfoil switches directions, a negative vortex forms at the trailing edge and quickly sheds into the wake (see Figure 6.9).  As the airfoil continues pitching down, high magnitude negative vorticity from the feeding shear layer continues shedding the wake.  At !/T = 0.50-0.557, a ÒconcentrationÓ of vorticity separates from the shear layer and rolls up into a vortex.  As positive vorticity forms at the trailing edge (!/T = 0.729-0.844), a third ÒconcentrationÓ of vorticity forms in the shear layer and eventually rolls up into a third vortex.   146  Figure 6.8.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 6.68, S = 30%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 147  Figure 6.9.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 6.68, S = 30%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075. 148 6.3.  Vortex Properties In order to quantify the effect of trajectory asymmetry on the wake structure for the three reduced frequencies, the vortex properties were computed.  As observed in Figure 6.1, a variety of vortical configurations were identified: two opposite sign vortices, one positive and two negative, and one positive and three negative.   In order to examine the effect of asymmetry on the vortex properties for different k, only the positive vortex and initially-formed negative vortex are considered (see Figure 6.10 for a schematic).  The vortex arrangement was quantified by the spatial location of the centroid in the transverse direction, the strength by the vorticity peak value and circulation, and the size by the core radius.  The properties for both vortices were quantified at x/c = 1.0.  This location was chosen since it was the measurement location from the experiment by Koochesfahani (1989).  Figure 6.10.  Schematic of vortex labeling for both S = 50% and S = 38% for vortex properties. The properties for the the two vortices of interest are shown in Figure 6.11.  The transverse centroid coordinate, yc, at S = 50% depends on the arrangements of the vortex array at the different reduced frequencies, such as traditional or reverse von K⁄rm⁄n street.  As asymmetry increases, yc decreases for a given k.  At S = 38%, each vortex is below the wake  149 centerline (i.e., yc < 0).  As asymmetry increases further, yc continues decreasing.  The peak vorticity magnitude increases as both reduced frequency and asymmetry increase for vortex P.  For vortex N, the peak vorticity magnitude decreases as asymmetry increases until S = 38%, then begins increasing in magnitude as asymmetry increases further.  It is unclear why this behavior occurs.  As the reduced frequency increases, the circulation magnitude increases in general.  More interesting behavior is seen in the variation with respect to S.  As S decreases, the circulation magnitude of vortex P increases, slightly drops, then increases for k = 5.2.  The opposite trend occurs with vortex N at this reduced frequency, where the circulation magnitude decreases, slightly increases, then decreases.  At k = 5.8 and 6.68, the circulation magnitude of vortex P increases and decreases for vortex N in a monotonic fashion.  Interestingly, the circulation magnitude of vortex N becomes nearly independent of k at S = 34-30%.  The core radius in general decreases as asymmetry increases.  At k = 5.2 and S = 38%, there is a jump in core radius, though it is unclear why this occurs.    150  Figure 6.11.  Vortex properties of positive vortex P (circles) and negative vortex N (triangles) for different values of S at different reduced frequencies.  The Reynolds numbers is 12,000 and the pitching amplitude is 2¡.  The streamwise spacing between vortices of the same sign, a, is calculated at x/c = 1.0 for both vortices using the convection speed and oscillation frequency.  The relative spacing between vortex P and N is !xP-N and is calculated when vortex P is at x/c = 1.0.  These definitions are shown schematically in Figure 6.12.  The transverse spacing is not considered here.  The vortex spacing data versus S is shown in Figure 6.13 for different reduced frequencies.  As the reduced frequency increases, both a and !xP-N decrease.  As S decreases, the a spacing slightly increases for both vortices at a given k while !xP-N decreases.   151  Figure 6.12.  Schematic of vortex spacing definitions for asymmetry trajectory.  Figure 6.13.  Streamwise spacing of centroid between opposite sign vortices (top) and same-sign vortices (bottom) for different values of S at different reduced frequencies.  The Reynolds numbers is 12,000 and the pitching amplitude is 2¡.  152 CHAPTER 7.  CONCLUSIONS AND FUTURE RECOMMENDATIONS The goal of this work was to investigate the effect of both Reynolds number and trajectory asymmetry on the forces and wake structure of a pitching NACA 0012 at low Reynolds number.  In order to accomplish this, high-order, two-dimensional computations were performed over a wide parameter space encompassing Reynolds numbers, reduced frequency, pitching amplitude, and asymmetry.  The wake structure, which is not typically of interest in most computational studies on flapping aerodynamic, was resolved using an overset grid system with very high resolution in the wake capable of capturing the evolution of the vortex array for up to four chord lengths from the airfoilÕs trailing edge.  The high resolution allowed for both qualitative observation and quantitative characterization to be undertaken.  Simulation of the incompressible limit was done using low freestream Mach numbers with the current compressible Navier Stokes solver.  Only two-dimensional simulations were performed, limiting to assumptions of laminar flow which are appropriate for this Reynolds number range.  Initial studies were conducted to determine the freestream Mach number requirements for both the symmetric and asymmetric trajectories.  It was found that freestream Mach numbers lower than those typically used in the literature were required to converge the various parameters of interest, such as force maximum and minimum values, average force values, and r.m.s. vorticity profiles in the wake.  The low Mach number requirements became more challenging when the airfoil trajectory became asymmetric.  However, with appropriate time-step size and enough subiterations, converged solutions were obtained.    153 7.1.  Reynolds Number Dependence The Reynolds number was shown to have a dramatic effect on the average forces, propulsive efficiency, and wake structure at the low Reynolds number range considered in this research.  As the Reynolds number increased, the average thrust coefficient decreased in magnitude from its starting static airfoil drag value, switched sign, and increased.  This increase in thrust caused the drag-to-thrust crossover reduced frequency to decrease monotonically with increasing Reynolds number and was a result of two mechanisms: 1) the static airfoil drag offset, which decreased with increasing Reynolds number, and 2) the shear stress contribution to drag, which increased with increasing reduced frequency and decreased as the Reynolds number increased.  This two-prong effect caused the behavior in the crossover reduced frequency.  It was also observed that the propulsive efficiency, lift peak-to-peak, and drag peak-to-peak also increased with Reynolds number, though the effect of Reynolds number on the lift and drag peak-to-peak amplitude was small.   Different pitching amplitudes ranging from 2¡ to 10¡ were also considered to demonstrate the effect of amplitude.  As the amplitude increased, the thrust increased for each Reynolds number and caused the crossover reduced frequency to decrease.  When re-scaled based on the Strouhal number, there was a pronounced collapse of the crossover data.  The collapse showed that the 2¡-6¡ data approximated the low-amplitude limit, while amplitudes greater than 6¡ showed deviation, particularly at the higher Reynolds numbers.  Scaling the Strouhal number by the crossover Strouhal number and the average thrust by the static drag offset caused the thrust coefficient data to become nearly independent of Reynolds number amplitude, with higher amplitudes showing small deviation.  The peak-to-peak amplitudes of lift and drag increased  154 with pitching amplitude, but could be nearly collapsed onto a single curve using amplitude scaling based on linear theory.   The propulsive efficiency was found to increase as Reynolds number increased.  These data nearly collapsed with amplitude using the Strouhal number.  Also interesting was that the propulsive efficiency increased more gradually than found in the literature and did not reach a well defined peak over St = 0.2-0.4 for each Reynolds number.   The wake data showed strong dependence of wake vortex configuration on Reynolds number.  At a reduced frequency k = 6, increasing Reynolds number from 2,000 to 22,000 showed a switch in the vortex arrangement from a traditional von K⁄rm⁄n street (drag) to a reverse von K⁄rm⁄n street (thrust).  The wake structure had a nearly aligned arrangement at a Reynolds number of 12,000 at this reduced frequency.  It was found that the vortex transverse spacing increased as the Reynolds number increased over the low k range.  As the reduced frequency increased, the transverse spacing switched from negative to positive for each Reynolds number and asymptoted to a nearly common value.   Analysis of the vortex transverse spacing data showed that the reduced frequency for vortex spacing crossover decreased as the Reynolds number increased.  The effect of pitching amplitude showed an increase in transverse and streamwise spacing with increasing amplitude.  However, only the lower Reynolds numbers produced traditional von K⁄rm⁄n streets at "o ! 4.  The resulting crossover reduced frequency decreased with increasing Reynolds number and could also be better scaled using a crossover Strouhal number definition.  It was also found that the vortex aspect ratio nearly collapsed when plotting the resulting data versus Strouhal number for a given Reynolds number.    155 7.2.  Trajectory Asymmetry Trajectory asymmetry showed a dramatic impact on the flow characteristics for the single Reynolds number and pitching amplitude considered.  As the asymmetry increased, the average lift became non-zero while there was a reduction in drag.  Behavior in average lift and the drag reduction depended on both reduced frequency and asymmetry.  In particular, as the reduced frequency increased, the non-zero lift switched in sign from positive to negative.  It was found that the behavior in average lift was directly tied to the domain integrated Lamb vector, called the vortex lift.  As the reduced frequency increased, the time duration and amount of negative vortex lift increased and eventually overtook the positive vortex lift.  Increasing asymmetry caused the peak-to-peak amplitudes of lift and drag to increase while the lift and drag fluctuation was primarily controlled by the pitching acceleration.   Trajectory asymmetry produced interesting flow features in the wake, such as multiple negative vortices shed during the slower pitch-down.  As the trajectory asymmetry increased, the distribution of vorticity in the wake became increasingly asymmetric as well.  This was quantified both in the trajectories of the vortices and the vortex properties, particularly the peak vorticity and circulation.  As asymmetry increased, the positive-sign vortex developed higher peak vorticity and circulation values than than of the negative-sign vortex formed at the trailing edge.  The relative streamwise spacing between opposite sign vortices showed more influence by trajectory asymmetry than the streamwise spacing between same-sign vortices.   7.3.  Future Recommendations Although this work was extensive, there do exist additional effects future studies should address.  Work by Dong, et al., (2006) and Bucholz and Smits (2008) suggest that the thrust coefficient, propulsive efficiency, and implicitly the crossover frequency depend on aspect ratio.   156 Since one of the assumptions in this work was two-dimensional flow, it would be of interest to investigate aspect ratio effects, with the presence of a wall boundary condition in the spanwise direction in particular, on the forces and wake structure for a pitching airfoil for different Reynolds numbers.  Since a single pitch-axis was used in this work, it would also be interesting to see how the forces and wake structure depend on pivot axis.  Since only a limited parameter space was considered for the asymmetric trajectory investigation, the effect of Reynolds number and pitching amplitude would be of interest for future study.  Finally, non-zero mean angles of attack would be of interest since natural flyers most likely utilize a non-zero mean angle of attack during flight.      157        APPENDICES   158 APPENDIX A.  SMOOTHED ASYMMETRIC MOTION TRAJECTORY A.1. Smoothed Trajectory Analytical Derivation The original three-part formulation of the asymmetric motion trajectory is given by equation (A.1), where T is the oscillation period, !" is the pitching amplitude, and S is the symmetry parameter (S = TPU/T).   Part 1: 0 " t < TPU/2    !#$%!"&'()$*+  Part 2: TPU /2 " t < T-TPU /2    !,$%-!"&'().-*$+-./ (A.1) Part 3: T Ð TPU /2 " t < T    !0$%!"&'()*$+-.  In order to apply a Gaussian filter, the angle of attack, !$, must be rewritten as a Fourier series, given by   !$%1"/21345&/)6+$273&'(/)6+$839# (A.2) Since the expression is an odd function, only the bn terms will be non-zero.  Equation (A.3) is then used to solve for the Fourier coefficients.  159  73%/+!$&'(/)6+$:$,;< (A.3) where n is the number of terms in the series solution.  In order to recover the asymmetric trajectory, velocity, and acceleration, at least 100+ terms are used in the series (see Figure A.1).    Figure A.1.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) coefficients for k = 6.68 and S = 38%, illustrating number of Fourier series terms. Using the relationships between the phase angle, !, and time, (i.e., ! = 2#t/T and d! = 2#/Tdt), equation (A.3) may be rewritten as  73%.)!=&'(6=:=,;< (A.4) Substituting the equation (A.1) into (A.4) yields   160  73%.)&'(=/*&'(6=:=>?<-&'(=-)/.-*&'(6=:=,;@>?>?2&'(=-/)/*&'(6=:=,;,;@>? (A.5) Solution of these integrals give the angle of attack as  !=%!"73&'(6=839# (A.6) where  73%+#2+,2+02+A2+B2+C (A.7) Terms T1-T6 in equation (A.7) are written as equations (A.8)-(A.13).  It is worth mentioning that these terms are derived only in terms of the symmetry parameter, S.    +#%/*.-/*6&'().-/*6/ (A.8)  +,%/*.2/*6&'().2/*6/ (A.9)   +0%/.-*.-/.-*6&'().-D62/*6/2&'().2/*6/ (A.10)   +A%/.-*.2/.-*6&'().2D6-/*6/-&'()-.2/*6/ (A.11)  161   +B%/*.-/*6-&'(/)6-&'()-.-D62/*6/ (A.12)  +C%/*.2/*6&'(/)6-&'()-.2D6-/*6/ (A.13) The Gaussian filter, E$%#,;FGH@IJGGKG, is then applied to the angle of attack according to   !$%E$-$L$:$8@8%./)M,H@IN@NG,FG!$:$8@8 (A.14) Substitution of equation (A.6) into (A.14) gives   !$%./)M,H@IN@NG,FGH@,;3G,FOG73&'(6/)$+839#:$8@8 (A.15) The smoothed trajectory is thus given by  !$%!"73H@,;3G,FOPQG&'(6=839# (A.16)   162  Figure A.2.  Angle of attack (top), pitching velocity (middle), and pitching acceleration (bottom) coefficients for k = 6.68 and S = 38%, illustrating influence of ". A.2. Smoothed Trajectory Results The effect of trajectory smoothing on the lift time histories is shown in Figure A.3.  As the smoothing increases, there is an increasingly dampened response in the forces around the phase where the angle of attack changes direction the average forces is a slight decrease in thrust and a slight alteration in CL trend with respect to S curves as smoothing increases, as shown in Figure A.4.  The wake structure was also slightly affected.  Average u, r.m.s. u, and r.m.s. v profiles measured one chord from the trailing edge are shown in Figure A.5.  As the smoothing increases, the average u deficit drops slightly in magnitude.  There is also a slight increase in r.m.s. u and v at y/c = -0.2.  163  Figure A.3.  Time history of lift lift (left) and drag (right) coefficients for k = 6.68 and S = 38%, illustrating influence of ".   Figure A.4.  Average lift (left) and thrust (right) coefficients for k = 6.68 and different asymmetries, illustrating influence of ".  164  Figure A.5.  Average u (left), r.m.s. u (middle), and r.m.s. v (right) for k = 6.68 and S = 38%.   165 APPENDIX B.  COMPUTATIONAL METHOD The numerical approach in this work is implicit large eddy simulation (ILES) using the high-order FDL3DI Navier-Stokes solver developed at AFRL, located at Wright-Pattern Air Force Base (Gaitonde and Visbal, 1998; Visbal and Gaitonde, 1999).  The details of the method are now discussed in this chapter, which is organized as follows: Section B.1 provides the governing equations utilized while B.2 presents the handling of the metric identities.  Section B.3 details the time marching scheme.  Sections B.4 and B.5 present the spatial discretization and filter, respectively.  Section B.6 contains the description of the high-order interpolation process, solver parallelization discussed in Section B.7.  Note that section follows similar discussions as presented by Galbraith (2009), Garmann (2010), and Garmann (2013) and provides a summary of these descriptions.  Finally, only two-dimensional simulations are considered in this work, whereas this discussion corresponds to the more generalized three-dimensional case. B.1. The Governing Equations  The governing equations being numerically solved are the full, unsteady, compressible, three-dimensional, unfiltered Navier-Stokes equations.  To generalize the Navier-Stokes equations and make them more efficient for a wider range of geometries, appropriate coordinate transformations from Cartesian (x, y, z, t) to curvilinear coordinates (#, $, %, &) can be made as follows:  R%$S%STUVUWU$X%XTUVUWU$Y%YTUVUWU$ (B.1) Application of these transformations allows the Navier-Stokes equations to be written in the strong conservative form, shown below in equation (B.2).    166  ZZR[\2Z]^ZS2Z_^ZX2Z`^ZY%.aHZ]bZS2Z_bZX2Z`bZY (B.2) The solution vector of conservative variables, [, is defined by  [%cUcdUceUcfUcgO (B.3) In the above expression, ' is the density, u, v, and w are the Cartesian velocity components, and E is the specific total energy defined by   g%+hh-.i8,2./d,2e,2f, (B.4) where, T is the temperature, ( is the ratio of specific heats, and M) is the freestream Mach number (M) = U)/a), where U) and a) are the freestream velocity and speed of sound, respectively).  The transformation Jacobian is J, where J = *(#, $, %, &)/ *(x, y, z, t).   The inviscid fluxes (]^, _^, & `^,) are defined by   ]^%.\cjcdj2Sklcej2Smlcfj2Snlcg2lj-SNl_^%.\cocdo2Xklceo2Xmlcfo2Xnlcg2lo-XNl`^%.\cpcdp2Yklcep2Ymlcfp2Ynlcg2lp-YNl (B.5)  167 The subscripts x, y, and z denote a partial derivative.  The directional contravariant velocities produced by the coordinate transformation, U, V, and W, in equation (B.5) are  j%SN2Skd2Sme2Snfo%XN2Xkd2Xme2Xnfp%YN2Ykd2Yme2Ynf (B.6) For simplicity, indicial notation is used to represent the Cartesian coordinates and velocities, as well as the curvilinear coordinates, as follows  T#UT,UT0%TUVUWd#Ud,Ud0%dUeUfS#US,US0%SUXUY (B.7) Using this indicial notation, repeated indices i, j, k, and l refer to summations while qrs is the Kronecker delta function.  The viscous fluxes, ]b, _b, and `b, in indicial notation are given by  ]b%.\tSkuRr#SkuRr,SkuRr0SkudsRrs-vrw_b%.\tXkuRr#XkuRr,XkuRr0XkudsRrs-vrw`b%.\tYkuRr#YkuRr,YkuRr0YkudsRrs-vr (B.8) The shear stress components &ij are defined by equation (B.9) upon implementing StokeÕs hypothesis for bulk viscosity (i.e., + = -2/3µ).  Rrs%xZSyZTsZdrZSy2ZSyZTrZdsZSy-/zqrsZS{ZTsZdyZS{ (B.9) In addition, the heat flux tensor (,i) is  vr%-.h-.i8,x|}ZSsZTrZ+ZSs (B.10)  168 In the above expression, µ is the dynamic viscosity while Pr is a constant Prandtl number (Pr = 0.72 for air).  In order to close the set of equations, the Ideal Gas Law (equation (B.11)) and SutherlandÕs law for viscosity (equation (B.12), where Sc = 0.38 for air) are used.  All flow variables are normalized by their freestream value, except pressure which is normalized by twice the dynamic pressure (i.e., ')U)2).  The length scale is the chord length, c, and the physical time has been normalized by U)/c.  l%c+hi8, (B.11)  x%+0~,.2*+2* (B.12) Although these equations correspond to the three-dimensional Navier-Stokes equations, only the limiting case of two-dimensional flow is considered due to prohibitiveness of three-dimensional computations given the very large parameter space considered in this work.  However, the assumption of two-dimensional flow is supported by experiments (Koochesfahani, 1989).   B.2. Treatment of Metric Identities By putting the transformed Navier-Stokes equations back into strong conservative form, four metric identities related to the surface and volume of cells are invoked.  The surface conservation identities are   169  •#%Sk\†2Xk\‡2Yk\…•#%Sk\†2Xm\‡2Ym\…•0%Sn\†2Xn\‡2Yn\… (B.13) while the volume conservation identity, known as the geometric conservation law (GCL; Thomas and Lombard, 1979), is  •A%.\—2S—\†2X—\‡2Y—\… (B.14) These relationships constitute a differential statement of surface and volume conservation for a closed cell (Visbal and Gaitonde, 2002).  Only identities I1-I3 are applicable when the computational domain is time-invariant.  These identities must be properly treated in order to satisfy freestream preservation for the numerically discretized, transformed governing equations.  Although methods proposed by Pulliam and Steger (1978) and Vinokur (1989) are effective in two-dimensional flows using low-order explicit schemes, they cannot be extended to high-order explicit or compact schemes in their present form.  An alternative method for enforcing the metric identities was given by Thomas and Lombard (1979).  This procedure involves rewriting the metric relations into their respective conservative form before spatial discretization, given by  170  Sk%\V‡W…-V…W‡Sm%\V‡T…-V…T‡Sn%\T‡V…-T…V‡Xk%\V…W†-V†W…Xm%\W…T†-W†T…Xn%\T…V†-T†V…Yk%\V†W‡-V‡W†Ym%\W†T‡-W‡T†Yn%\T†V‡-T‡V† (B.15) It has been shown by Visbal and Gaitonde (2002) that using the conservative form of the metric identities satisfies both freestream preservation and metric cancellation on general three-dimensional, curvilinear meshes using both low- and high-order discretization schemes.   In order to numerically satisfy the GCL, Visbal and Gaitonde (2002) split the time-derivative term from equation (B.2) using the chain rule.  ZZR[\%.\Z[ZR2[ZZR.\ (B.16) The first term is the inverse transformation Jacobian (J-1) and is evaluated using the instantaneous grid coordinates while the second term requires special treatment.  Rather than computing the time derivative of the inverse Jacobian at various time-levels, the GCL is directly invoked as  .\—%-S—\†2X—\‡2Y—\… (B.17)  171 where the time-dependent metric terms, #&, $&, & %&, are  S—\%-T—Sk\2V—Sm\2W—Sn\X—\%-T—Xk\2V—Xm\2W—Xn\Y—\%-T—Yk\2V—Ym\2W—Yn\ (B.18) In the above expressions, x&, y&, and z& are the grid speeds that can be determined either analytically or numerically.  Visbal and Gaitonde (2002) showed that using analytical grid speeds provides effective metric cancellation and freestream preservation.  The spatial derivatives of the time-dependent metrics are computed using the high-order compact difference scheme while the second term in (B.16) is moved to the right-hand side of the time-marching algorithm, found in the next section.   B.3. Time-Marching B.3.1: Beam-Warming Method Time accurate solutions of equation (B.2) are obtained numerically using the linearized, approximately factored integration scheme developed by Beam and Warming (1978).   Consider the three-dimensional Euler equations in a Cartesian coordinate system, where the flux terms, ], _, and `, are functions of the solution vector [.  Z[Z$%Z]^[ZT2Z_^[ZV2Z`^[ZW (B.19) Beam and Warming (1978) showed that the solution vector could be advanced in time by  172  –[3%–$.2=ZZ$–[32–$.2=ZZ$[32=.2=–[3@#2ƒ./2=–$,2–$0%a (B.20) The solution vector is in ÒdeltaÓ form, where [3%[6⁄$ and ⁄[%[3‹#-[3.  The current time-level is n+1 and the two previous time-levels are n and n-1.  The time-step is -t while the constants !1 and !2 set the order of accuracy. Various explicit and implicit time integration methods can be recovered through proper selection of !1 and !2.  The Euler, first-order implicit scheme is obtained from ! = 0.0 and ! = 0.5 produces the three-point backward, second-order implicit scheme.   B.3.2. Beam-Warming Method: Linearization and Approximate Factorization To improve the solution efficiency of the algorithm, Beam and Warming (1978) linearized and approximately factored equation (B.20); the outcome is presented in equation (B.21).   •2–$.2=ZZT›^•2–$.2=ZZV−^•2–$.2=ZZW‰^–[3%a (B.21) The terms AI, BI, and CI are the inviscid flux Jacobians in Cartesian coordinates (refer to Beam and Warming, 1978, or Anderson, et al., 1984 for further details) and I represents the identity matrix.  Equation (B.21) allows for the inversion of three one-dimensional block tridiagonal matrices instead of a single, three-dimensional matrix, which is more efficient to solve.  The inversion process then involves three one-dimensional sweeps, as shown in equation (B.22).  In this expression, ⁄[#3and ⁄[,3 represent intermediate steps in the inversion process.    173  •2–$.2=ZZT›^–[3,%a•2–$.2=ZZV−^–[3,%–[3#•2–$.2=ZZW‰^–[3#%–[3<–[3%–[3< (B.22) However, since the inversion process involves the inversion of block matrices, it is still rather cumbersome.  Therefore, it is still desirable to improve the efficiency of the algorithm.  B.3.3. Diagonalization by Pulliam and Chaussee (1981) To further improve computational efficiency of the Beam-Warming method for application to the Navier-Stokes equations written in general, curvilinear coordinates, Pulliam and Chaussee (1981) introduced a method of diagonalizing the Beam-Warming method by using the eigensystem of the inviscid flux Jacobians.    ›^%+†„^U†+†@#−^%+‡„^U‡+‡@#‰^%+…„^U…+…@# (B.23) I In the above expression, „^U† „^U‡, and „^U… denote the diagonal vector of eigenvalues in AI, BI, and CI and +†, +‡, +…, +†@#, +‡@#, and +…@# are matrices whose columns are the eigenvectors of AI, BI, and CI.  By assuming that the matrices +†, +‡, +…, +†@#, +‡@#, and +…@# are locally constant, they may be taken out of the difference operators in equation (B.21).  Implementing this diagonalization procedure in equation (B.22) yields  174  +†•2–$.2=ZZS„^U†+†@#+‡•2–$.2=ZZV„^U‡+‡@#+…•2–$.2=ZZW„^U…+…@#–[3%a (B.24) The diagonalization leads to an uncoupling of the block tridiagonal system into a system of five independent, scalar tridiagonal equations that are solved simultaneously.  This system of equations may be solved using three one-dimensional sweeps, as detailed by equation (B.25).    •2–$.2=ZZT„^U†+†@#–[3,%+†@#a•2–$.2=ZZV„^U‡+‡@#–[3#%+‡@#+†–[3,•2–$.2=ZZW„^U…+…@#–[3<%+…@#+‡–[3#–[3%+…–[3< (B.25) By replacing the block tridiagonal matrix inversions of the Beam-Warming method with scalar tridiagonal matrix inversions, the computational efficiency is dramatically increased.    The diagonalization process, however, is valid only for the Euler equations since the viscous flux Jacobian is not simultaneously diagonalizable with the inviscid flux Jacobian (Pulliam, 1985).  In order to simulate the full Navier-Stokes equations, the viscous terms must be included in some fashion.  Pulliam (1985) discussed different options for including the viscous terms in the diagonalized Beam-Warming method.  The recommended option is to include a diagonal term that is an estimation of the viscous flux Jacobian eigenvalues in the implicit operator.  These estimations („bU†, „bU‡, and „bU†) are given by equation (B.26), where D is a diagonal operator.    175  „bU†%“xcSk,2Sm,2Sn,„bU‡%“xcXk,2Xm,2Xn,„bU†%“xcYk,2Ym,2Yn, (B.26) The estimation of the viscous flux Jacobian only has to be an approximation of the explicit right-hand-side as long as the solution is converged with subiterations.   In order to enhance the overall stability of the numerical scheme, high-order dissipation terms (Jameson, et al., 1981; Pulliam, 1986) are appended to the implicit operator as follows. The dissipation terms are added to the diagonalized implicit operator in each of coordinate direction and leads to a pentadiagonal system of scalar equations, which is less computationally expensive than the original block tridiagonal system.     B.3.4. Final Time-Marching Expression  The final time-marching algorithm written in ÒdeltaÓ form is given by equation (B.27).  The time-marching algorithm is supplemented by Newton-like subiterations to maintain second-order accuracy while reducing errors caused by the linearization, approximation factorization, and diagonalization, and explicit implementation of physical boundary conditions (Visbal, et al., 2003).  Typically, three subiterations are used (Visbal and Rizzetta, 2002; Visbal and Gaitonde, 2002; Visbal, et al., 2003).   176  .\”‹#2⁄R.2=q†,Z]^”Z[-.aHZ]b”Z[\”‹#.\”‹#2⁄R.2=q‡,Z_^”Z[-.aHZ_b”Z[\”‹#.\”‹#2⁄R.2=q…,Z`^”Z[-.aHZ`b”Z[⁄[%-⁄R.2=.\”‹#2.2=[”-.2/=[3-=[3@#⁄R2[”.\—”2q†C]^”-.aH]b”2q‡C_^”-.aH_b”2q…C`^”-.aH`b” (B.27) In equation (B.27), ⁄R is the time-step size, p denotes subiteration step, and ⁄[%[”‹#-[”, where [”‹# corresponds to the p + 1 approximation of the n + 1 time level and p denotes the subiteration level within the current time-step.  In the first subiteration, p = 1 corresponds to [3.  As p # $, [”#[3‹#.  The inviscid flux Jacobians are Z]^~Z[, Z_^~Z[, and Z`^~Z[ while Z]b~Z[, Z_b~Z[, and Z`b~Z[ are the viscous flux Jacobians.  First-order Euler-implicit and second-order three-point backward schemes result when ! = 0.0 and ! = 0.5, respectively.   The . operator on the left- and right-hand-side of the equation represents spatial derivatives computed using finite difference schemes in the direction of the subscript (i.e., #, $, and %) of the indicated order of accuracy, denoted by the number in the superscript.  For example, q†, denotes second-order central difference scheme in the # direction while q†C is the sixth-order compact scheme in the same direction.  Since the solution is converged with subiterations, the lower order of accuracy on the implicit side does not impact the final solution (Garmann, 2013).     177 Not shown are second- and fourth-order, nonlinear artificial dissipation terms (Jameson, et al., 1981; Pulliam, 1986) that are appended to the implicit, left-hand-side for enhancing stability.  This results in a scalar, pentadiagonal system of equations that is more computationally expensive than solving a scalar tridiagonal system of equations but dramatically less expensive than solving the original block tridiagonal system.  The use of subiterations eliminates the impact of the dissipation terms on the final solution, and thus can be chosen specifically for stability (Visbal and Gaitonde, 2002).  The implementation of this time-marching method has been successfully applied to a variety of flows (e.g. Visbal and Rizzetta, 2002; Visbal and Gaitonde, 2002; Visbal, et al., 2003; Rizzetta, et al., 2008).  The second-order time marching scheme is used for this research. B.4. Spatial Discretization As previously mentioned, the spatial derivatives on the explicit side of the time-marching method are computed with a high-order, compact finite difference scheme, first described by Lele (1992) to obtain spectral-like resolution.  For simplicity, the description of the compact difference scheme will be restricted to a one-dimensional mesh of N number of points with uniform grid spacing h (illustrated in Figure B.1).   Figure B.1.  Notation for the interior and boundary points of a 1-D mesh. For a discrete function known at any location (i.e., conservative variable, primitive variable, metric, etc.), /i = /(xi) where Tr‘.U’, the finite difference approximation of /i can be determined up to sixth-order accuracy using a five point interior stencil (i.e., i-2 thru i+2 in the above figure).  Special treatment is required for computing the spatial derivatives at the  178 boundary points.  Note that a coordinate transformation must be used to remap physical, curvilinear coordinates to the computational, rectangular coordinates for calculating the spatial derivatives.  The following discussion will describe interior schemes in more detail.  B.4.1. Interior Points  For the interior points (shown in Figure B.2), the finite difference approximation of variable / at point i can be expressed as a linear combination of the functional values of / at point i and the surrounding points.    Figure B.2.  Notation for the interior point stencil of a 1-D mesh. For a five-point stencil, this combination is represented by (B.28)   ‚™r@#fi2™rfi2‚™r‹#fi%1™r‹#-™r@#/fl27™r‹,-™r@,Dfl (B.28) where /Õ is the derivative of the points and the coefficients a, b, and 0 define the order of accuracy of the compact scheme.  In order to obtain the values of the coefficients, Taylor series expansions for all functional values and their derivatives about point i are first determined, shown in equations (B.29)-(B.31) in compact notation form.    179  ™r‹#-™r@#/fl%fl,3/62.Ł™r,3‹#839< (B.29)  ™r‹,-™r@,Dfl%/fl,3/62.Ł™r,3‹#839< (B.30)  ™r‹#-™r@#%fl,3/6Ł™r,3‹#839< (B.31)  Next, these Taylor series expansions are directly plugged into equation (B.28), leading to the following after simplification  ™rfi2/‚fl,3/62.Ł™r,3‹#839<%12/,37/fl,3/62.Ł™r,3‹#839< (B.32) By setting terms of the appropriate order to zero, the following three equations are produced for determining a, b, and 0:   ƒfl,Œ.-12/‚-7%tƒflAŒ-12Š‚-D7%tƒflCŒ-12.t‚-.Š7%t (B.33) Solutions to these three equations lead to a family of both explicit (0 = 0) and compact (0 % 0).  By setting 0 = 0, the left-hand side of (B.32) is decoupled, which results in an explicit derivative at point i.  Explicit schemes are one order of accuracy lower than the stencil size.  When 0 % 0, the derivative at point i is coupled with the derivatives at points i+1 and iÐ1 and requires the solving of an implicit tridiagonal system.  By setting b = 0 and b % 0, fourth-order and sixth-order compact schemes can be obtained, respectively.  Thus, the compact schemes are  180 one order of accuracy higher than the stencil size.  A compilation of these coefficients and their accompanying finite difference schemes can be found in Table B.1, where ÒEÓ and ÒCÓ denote explicit and compact schemes, respectively, and the number defines the order of accuracy.   Table B.1.  Coefficients for the interior finite difference schemes. Scheme ! a b Stencil E2 0 1 0 3 points E4 0 4/3 -1/3 5 points C4 1/4 3/2 0 3 points C6 1/3 14/9 1/9 5 points  In order to quantify the resolution capability of the compact and explicit schemes, a wavenumber analysis is performed.  Assuming /(x) is periodic over the domain [0, L], the function can be written as a Fourier series,   ™T%™yŸŽı/)łœTšž~,y9@ž~, (B.34) In the above equation, T‘tUš€, ł%-., k is the physical wave number, ™y is the Fourier coefficient.  For convenience, a scaled wavenumber (1 = 2#kh/L) and coordinate (s = x/h) are introduced.  With this transformation, (B.34) can be rewritten as (B.35).   ™¡%™yŸŽıł¢¡ž~,y9@ž~, (B.35) Taking the derivative of (B.35) results in   ™fi¡%ł¢™yŸŽıł¢¡ž~,y9@ž~,%ł¢™¡ (B.36)  181 If a finite difference approximation is made for the exact derivative, it is assumed a modified wavenumber (¢fi) will result in the following relationship between the modified wavenumber and the original Fourier coefficients as  ™y£¤%ł¢™y (B.37)   In order to obtain a relationship between the modified wavenumber and the coefficients of the compact scheme, the Fourier expansion using the scaled coordinates, and its derivative, are inserted into equation (B.28).  Through algebraic manipulation, the modified wavenumber can be related to scaled wavenumber and the compact scheme coefficients by    ¢fi%1I&'(¢27~/I&'(/¢.2/‚45&¢ (B.38) A comparison between the modified wavenumber and the scaled wavenumber gives a means of quantifying the resolving capabilities of the finite difference approximations with respect to exact differentiation.  The modified wavenumber for various finite difference schemes plotted against the scaled wavenumber is given in Figure B.3.  As seen in the figure, the compact schemes (dashed lines) exceed the resolution capability of the explicit schemes (dashed-dotted lines).   Wavenumbers in which the modified wavenumber deviates from the exact differentiation are under-resolved and can lead to errors.    182  Figure B.3.  Modified wavenumber versus scaled wavenumber for different finite difference schemes. It is worth noting that in the original formulation by Lele (1992), eighth- and tenth-order accurate compact schemes can be obtained with a seven-point stencil.  However, these two higher-order compact schemes require the solution of a pentadiagonal system.  The overall resolving efficiency increases from 63% for the sixth-order compact scheme to 70% and 73% for eighth- and tenth-order schemes, respectively.  Thus, from a computational standpoint, one must make trade-off between increased accuracy and decreased efficiency.  According to Garmann (2013), the sixth-order compact scheme offers an attractive balance between accuracy and computational cost.   B.4.2. Boundary Points As previously mentioned, special care must be taken when computing the spatial derivatives at the boundary points due to the lack of a centered stencil (Gaitonde and Visbal, 1998).  In order to maintain the tridiagonal system of the interior compact scheme, special one-sided equations must be used at the boundary points.  For point i = 1 and 2, the formulas are  183 given by equations (B.39) and (B.40), respectively.  In equation (B.45), the coefficients are set to be symmetric about point 2 (i.e., 021 = 022 % 0).    ™#fi2‚#™,fi%.fl1#™#27#™,2¥#™02:#™A2H#™B2¦#™C2E#™§ (B.39)  ‚,#™#fi2™,fi2‚,#™0fi%.fl1,™#27,™,2¥,™02:,™A2H,™B2¦,™C2E,™§ (B.40) A Taylor series expansion is done about points 1 and 2 and coefficients of like-order terms are matched in a similar way as done for the interior scheme.  This produces a series of equations whose solution yields the coefficients for various order of accuracy.  A summary of the coefficients for point 1 and 2 are given in Tables B.2 and B.3, respectively.  Note that f2 and g2 are zero. Table B.2.  Coefficients for boundary point 1. Scheme !1 a1 b1 c1 d1 e1 Stencil C2 1 -2 -2 0 0 0 2 points C3 2 -5/2 2 1/2 0 0 3 points C4 3 -17/6 3/2 3/2 -1/5 0 4 points  Table B.3.  Coefficients for boundary point 2. Scheme !21 !22 a2 b2 c2 d2 e2 Stencil AC4 1/4 1/4 -3/4 0 3/4 0 0 2 points AC5 3/14 3/14 -19/28 -5/42 6/7 -1/14 1/84 3 points  B.5. Spatial Filtering Spatial accuracy was discussed in the previous section.  An equally necessary issue to address is solution stability.  Since the high-order compact finite difference scheme is based on a centered stencil and thus non-dissipative, it is susceptible to numerical instabilities due to the growth of unstable frequency modes (Gaitonde and Visbal, 2000) that originate from mesh non- 184 uniformity, boundary condition approximations, and non-linear flow features.  In this work, a high-order low-pass Pad”-type filter, developed by Gaitonde, et al. (1999) provides dissipation at high wavenumbers only where the spatial discretization already exhibits significant dispersion errors.  The filter is applied to the conservative variables during each subiteration or time-step.  Since the filter is highly discriminating, it only dampens the poorly evolved high-frequency content of the solution (Gaitonde and Visbal, 1999; Gaitonde and Visbal, 2000; Visbal and Rizzetta, 2002; Visbal, et al., 2003).   The high-order filter is one-dimensional and is applied sequentially one direction at a time.   Before filtering a subsequent direction, the solution is updated with the filtered values.  To minimize potential bias, this sequence is alternated between various permutations (Visbal and Gaitonde, 1999).  The filtering sequence alternates once per time-step; thus the same sequence is applied after each subiteration of a specific time-step.  Next, the interior filter will be discussed in greater detail, with a brief description of the boundary filter.   B.5.1. Interior Points  The high-order filter is derived in a similar way as the high-order compact scheme and is based on templates presented by Lele (1992) and Alpert (1981).  The filtered values ™ of the variable ™ are obtained by solving the following tridiagonal system of equations    ‚¨™r@#2™r2‚¨™r‹#%1©/™r‹©-™r@©ž©9< (B.41) In the above expression, ‚¨ is a free parameter that adjusts the amount of filtering.  A 2N-order filter is obtained from a stencil consisting of 2N+1 interior points.  The coefficients in (B.41) are  185 found matching Taylor coefficients.  The Taylor series expansions of the terms on the left- and right-hand sides of equation (B.41) become equations (B.42) and (B.43), respectively.  ‚¨™r@#2™r2‚¨™r‹#%™r2‚¨fl,3/6Ł™r3ž39< (B.42)  1©/™r‹©-™r@©ž©9<%ªfl,3/6Ł™r,3ž39< (B.43) Matching Taylor series coefficients then produces the following relation  t,32‚¨%ª,31©ž39< (B.44) The resulting filter is thus non-dispersive, causing no wavenumbers to be amplified.  If 0f is left as a free parameter, then there exist N+1 equations with N+2 unknowns (i.e., ao, a1, a2, Éan).  In order to close the set of equations to solve for the filter coefficients, a Fourier analysis similar to that done for the compact scheme is performed.  The functional value and its corresponding filtered value have Fourier expansions, given in equations (B.45) and (B.46), respectively, where the scaled wavenumber is used.  ™¡%™yŸŽıł¢¡ž~,y9@ž~, (B.45)  ™¡%™yŸŽıł¢¡ž~,y9@ž~, (B.46)  186 Where are ™y and ™y are the Fourier coefficients of the unfiltered and filtered variables.  The spectral function of the filter equation is then determined substituting these expressions into equation (B.41), centering about ¡%t, and performing algebraic manipulation.  This results in  *]¢%™y™y%1©45&ª¢ž©9<.2/‚¨45&¢ (B.47) In order to close the set of equations, the highest frequency mode (SF(#)) is set to zero, thus producing  *])%-.©1©ž©9<%t (B.48) Thus, the filter coefficients can be solved with the following system of equations (B.49).    *])%tŒt%1<-1#21,-1021A-1Bƒfl,Œ.2/‚¨%1<21#21,21021A21BƒflAŒ/‚¨%1<21#2/,1,2z,102D,1A2«,1BƒflCŒ/‚¨%1<21#2/A1,2zA102DA1A2«A1Bƒfl¬Œ/‚¨%1<21#2/C1,2zC102DC1A2«C1Bƒfl#<Œ/‚¨%1<21#2/¬1,2z¬102D¬1A2«¬1B (B.49) Up to tenth-order accurate filter (N = 5) can be obtained with an eleven-point stencil.  A summary of the coefficient values for various orders of accuracy is presented in Table B.4.   187  Table B.4.  Interior point filter coefficients. Scheme ao a1 a2 a3 a4 a5 Stencil F2 !"#$%# !"#$%# 0 0 0 0 3 points F4 &"'$%( !"#$%# )!"#$%( 0 0 0 5 points F6 !!"!*$%!' !&"+,$%+# )!"#$%( !)#$%+# 0 0 7 points F8 -+".*$%!#( ."!($%!' )!"#$%!' !)#$%!' )!"#$%!#( * 9 points F10 !-+"!#'$%#&' !*&"+*#$%#&' )!"#$%', ),&"-*$%&!# )&"!*$%#&' )!"#$%&!# 11 points   188 The above analysis requires that !"#$%&'%"#$ in order for the spectral function to range from 0 (completely filtered) to 1 (unfiltered).  A special case occurs for &'("#", where explicit filtering is produced due to the left-hand-side of equation (2.46) being decoupled from the right-hand-side.  The impact of filter order and decreasing &' can be seen in Figure B.4, where the black line represents no filter being applied to ).  The greater the deviation from the black line, the greater the attenuation.  By increasing the order of accuracy, the filter becomes more discriminating in wavenumbers being attenuated.  By decreasing the free parameter towards -0.5, the filter increasingly attenuates lower frequencies.  Figure B.4.  Effect of filter order (left) and &' (right) on filter characteristics. It is important to remember that the filter is only supposed to eliminate wavenumbers unresolved by the compact scheme.  Thus, it is necessary that the selection of the filter be such that the resolved wavenumbers by the compact scheme remain unfiltered.  An example of this is presented in Figure B.5, where the compact scheme and filter are sixth- and eighth-order, respectively.  When the compact scheme (dashed-blue line) begins deviating from the exact derivative (solid-blue line), the wavenumbers become unresolved.  To ensure these resolved  189 wavenumbers remain unfiltered, the filter is set to be two-orders of accuracy higher than the compact scheme.  This spatial filter (F8/&'("#*") is used in the current simulations.  Figure B.5.  Resolution comparison between the compact scheme C6 and filter F8 (&' = 0.400). B.5.2. Boundary Points  Since the interior filter requires a large stencil (2N+1 points), special treatment is required at the boundaries where the interior stencil would extend beyond the boundaries.  For example, for an F8 filtering scheme, a nine-point stencil would not be able to properly span points 1 to 4 and points N-3 to N.  In order to address this issue, two approaches have been developed.    The first approach, presented by Gaitonde and Visbal (2000), is to use one-sided filter formulas to maintain the tridiagonal nature of the filter scheme.  However, the spectral function of the one-sided filter formula will become complex when 0 < &' < 0.5.  This could undesirably introduce artificial dispersion, where the real component of the spectral function exceeds unity over a range of wavenumbers and, thus, amplifying them (see Figure B.6).  As &' approaches  190 0.5, the imaginary component and the amount of excess over unity decreases.  This amplification increases with filter order of accuracy and &'.  The second approach is using a lower order central filter and higher value for &' to produce spectral characteristics similar to the higher order interior point filters (see Figure B.7).  At points very close to the boundary, one-sided formulas are still required using the second approach.  Regardless of approach, however, no filtering is applied at the boundary points 1 and N since the values at these points are determined from the boundary conditions.  Figure B.6.  Spectral filter response for the one-sided boundary point 2 formula for various filter orders and &'.  191  Figure B.7.  Spectral filter response for boundary point 2 for a one-sided filter (F4) and centered filter (F8) with similar spectral characteristics. For the current work, the order of the filter was reduced from eighth-order to fourth-order towards the boundary to maintain a centered stencil while the filter parameter &' increased towards 0.5.  Point 2 used a one-sided filter.  The coefficients and filter order are summarized in Table B.5.  Similar filtering schemes were adopted by Garmann (2011), Galbraith (2011), and Garmann (2013).   Table B.5.  Filter scheme utilized in current work. Coordinate  Pt. 2, N-1 Pt. 3, N-2 Pt. 4, N-3 Pt. 5, N-4 ! Filter "f F4 0.450 F4 0.400 F6 0.400 F8 0.400 # Filter "f F4 0.450 F4 0.400 F6 0.400 F8 0.400 $ Filter "f F4 0.450 F4 0.400 F6 0.400 F8 0.400    192  B.6:  Grid Connectivity and High-Order Interpolation Although structured meshes are very useful in their suitability with high-order compact schemes and filters due to their larger required stencils and tridiagonal nature, single structured meshes are often insufficient when investigations are done on complicated geometries or where additional resolution is desired in specific regions.  In addition, structured meshes often inadvertently cluster grid points in regions of less interest, resulting in undesirable computational overhead.   To remedy this situation, multiple grids comprise the computational domain in an overset approach.  Connectivity between the grids must be established in a pre-processing and is done using Pegasus 5 (Suhs, et al., 2002).  First, donor points and receiver points are identified.  Here, donor points are points in which data is interpolated from and receiver points are points in which data is interpolated to.  For a schematic of the interpolation, please refer to Figure B.8.    Figure B.8.  Schematic of donors points (from Mesh 1) and receiver points (from Mesh 2). Next, it is necessary to establish the distance between the receiver point r and the base point of the donor stencil d in the computational space of the donor grid, termed the interpolation offsets (%).  The interpolation offsets are then computed using an inverse isoparametric mapping  193 with eight donor points (two in each coordinate direction). Trilinear interpolation is performed between the donor cells and the receiver points during the solution process.  The low-order interpolation coefficients are then made high-order using AFRL developed BELLERO (Sherer, et al., 2006).  It computes the high-order interpolation coefficient and offsets with the following procedure.  First, the donor point stencil obtained from Pegasus is expanded such that the number of points in each coordinate direction equals the desired interpolation order of accuracy (N).  The receiver point is made as centered as possible within the donor point stencil during the expansion process.  Next, the high-order offsets in each coordinate direction are obtained by solving the following set of isoparametric mapping equations (equation (B.50)) based on the expanded stencil.    +,-.,/.0(12314516789:6;<=,-!,/("89:4;<89:2;< (B.509) In the above expression, 17, 15, and 13 are the high-order interpolation coefficients in each direction as a function of the interpolation offset 0 while ,- and ,/are the physical (x, y, z)-coordinate locations for the donor and receiver points, respectively.  The high-order offsets are solved iteratively using NewtonÕs method as   0>?:(0>@A+A09:B+,-.,/.0 (B.51) The initial guess for 0 is the second-order offset calculated from Pegasus 5.  The final step is the calculation of the interpolation coefficients using the expression  194  1>C(!D8?>9:E!F!DGFG0C!H89:6;<6I4 (B.52) Finally, the interpolated flow variable at the receiver point r is calculated from   )/(12314516789:6;<=)JK?6.LK?4.MK?289:4;<89:2;< (B.53) where )/ is the interpolated scalar at the receiver point, )- is the value from the donor point, and the indices Id, Jd, and Kd represent the location of the base donor point in computational space (Sherer and Visbal 2007).  Sixth-order interpolation is used in the present work. B.7:  Solver Parallelization  Although computational resources have advanced very rapidly within the last decade, high-order calculations on grids consisting of hundreds of thousands to several million grid points over tens and hundreds of thousands of time-steps and subiterations continues to be very computationally expensive.  In order to increase computational efficiency, the calculations of the discretized Navier-Stokes equations are parallelized such that solutions are obtained through simultaneous calculations on multiple, independent processors (Morgan, et al., 2006).  This parallelization is accomplished by decomposing the computational domain of one or more grids into a series of smaller subdomains (or blocks), which assigned to different processors.  Since individual grids are decomposed into smaller, more manageable pieces, the exchange of information between blocks becomes critical to the solution process.   In order to pass information between blocks, a finite-point overlap between blocks is necessary.  To illustrate this procedure, Figure B.9 gives a schematic of an example block  195 decomposition for a one-dimensional mesh.  Consider a single one-dimensional mesh being decomposed about point i.  Due to the five-point stencil of the compact scheme, each block requires a minimum five-point overlap for transferring information (Gaitonde and Visbal, 2000).  This overlap consists of the point of interest and two fringe points (i.e., two points on either side of i; circles with dashed lines around them in Figure B.9).         Figure B.9.  Illustration of block decomposition with five-point overlap. At every time-step, the equations are solved independently and simultaneously in each block.  Once the solution is computed in a given block, it shares data with adjacent blocks in the following way.  The values of the flow variables at points 1 and 2 on Block 2 are set to be equal to points N-4 and N-3, respectively, from Block 1.  Similarly, the values at points N-1 and N on Block 1 are set to the values of points 4 and 5, respectively, from Block 2.  Points 3 on Block 2 and N-2 on Block 1, however, do not communicate directly.  Thus, the donor points (i.e., points N-4, N-3, 4, and 5) act as the boundary conditions for the receiver points (i.e., points 1, 2, N-1,  196 and N), thus ensuring a smooth transition of the overall flow solution across the different blocks of the computational domain while maintaining high-order accuracy.   The synchronization process between block boundary updates is accomplished using a Message Passage Interface (MPI) library (MPI Forum, 2012).  For consistency, boundaries that rely on the physical boundary conditions are also updated immediately after the MPI synchronization.  Although the additional overlap points produce computational overhead, the parallelization is significantly more efficient than solving the equations in the original domain.  Due to the usage of structured grids, this parallelization approach is easily extendable to higher domain dimensions.  In the present FDL3DI Navier-Stokes solver, two MPI and physical boundary condition updates are performed for each subiteration: once after the solution is updated from the time integration and once after the solution is filtered (Visbal and Gaitonde, 2000).   In order to optimize the efficiency of the solution, grid point decomposition should be as well balanced as possible across the various processors.  For the present studies, grid decomposition is performed using BELLERO (Sherer, et al., 2006).  The first step of the decomposition procedure is to determine how many blocks to generate and is done through an iterative procedure.  Next, various combinations of cuts are made in the blocks, also through an iterative procedure.  The combination that results in a minimum surface-to-volume ratio in order to minimize computational overhead is selected.  The third step is to then apply the ÒoptimalÓ cuts to the individual grids.  The final step is to account for any blanking, as well as enforce the minimum stencils necessary for the spatial filter.  Once the computational domain is decomposed, the grid connectivity and appropriate boundary conditions are established for all  197 blocks.  For this work, the entire computational domain is decomposed into seventy-one blocks with a minimum of 100 points in the !- and #- directions.    198 APPENDIX C.  CONVERGENCE STUDIES C.1. Symmetric Pitching Convergence Studies  In order to test the fidelity of the solutions, grid resolution, subiteration, and time-step convergence were conducted for the extreme conditions of Re = 12,000, M& = 0.015, k ! 12.0 and !o = 2¡.  The grid was coarsened from 5 million grid points to 3.6 million (Medium Grid) and 2.5 million (Coarse Grid) with similar grid point distributions as the baseline grid, the time-step was varied between 5x10-5 and 2.5x10-5, the number of subiterations was increased from 3 to 9.  To present convergence, time histories of lift and instantaneous vorticity profiles measured one chord length from the trailing edge, at a phase of 0.625T, are shown in Figure C.1.  The choice over instantaneous profiles over average or r.m.s. is that the average and r.m.s. profiles are indistinguishable for investigating grid, time-step, and subiterations.  Instantaneous profiles provide a more stringent examination of solution convergence.    In general, there are no visual differences, nor great quantitative differences, between the lift histories for the three studies.  The peak and average values of lift and drag also showed negligible quantitative differences.  The effect of grid resolution on the vorticity profiles show the medium and fine grids produce nearly identical profiles while the coarse grid does show some agreement with the other two grids, but the fine and medium grids show much better agreement.   Both time-steps produce the same vorticity profile, showing time-step insensitivity for this case.  The number of subiterations show much more sensitivity.  Three subiterations do not produce a converged solution.  Five or more subiterations result in nearly identical profiles.  At phases where there the maximum vorticity value of the profile is low (max("z(y))c/U&  ~ 1), there is noticeable background noise in the profiles (Figure C.2).   Increasing the number of  199 subiterations reduces the background noise in the profiles.  Based on these data, the fine grid, seven subiterations and 5.0x10-5 time-step size were selected.  Figure C.1.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at #/T = 0.625 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M& = 0.015, k = 11.9, and !o = 2¡, and S = 50%.  Figure C.2.  Instantaneous vorticity measured at x/c = 1.0 at #/T = 0.25 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M& = 0.015, k = 11.9, and !o = 2¡, and S = 50%.  To illustrate the effect of Mach number, three cases were considered: k = 6.68 and !o = 2¡ (St = 0.11), k ! 12.0 and !o = 2¡ (St = 0.20), and k = 6.68 and !o = 6¡ (St = 0.33).  The Mach number varied between 0.005 to 0.100.  For reference, common values found in the literature are  200 0.050 (Young and Lai, 2004) and 0.100 (Visbal, 2009).  Note that the time-step was 5.0x10-5 for all Mach numbers, except for M& = 0.005, which had a time-step of 2.5x10-5 to remove significant background noise seen when the time-step was 5.0x10-5.  To present convergence, time histories of lift and r.m.s vorticity profiles measured one chord length from the trailing edge are shown in Figure C.3.  Note that r.m.s. profiles were selected instead of instantaneous profiles since slight variation in the vorticity shed from the trailing edge by Mach number results in slight differences in the wake structure that cannot readily be compared when using instantaneous profiles.  The r.m.s. profiles show the net effect on the wake structure.   First, the lift histories are considered.  As the Mach number decreases, the peak lift in general decreases until it reasonably converged at M& = 0.005-0.025 for all three cases.  The phase of maximum and minimum lift also shifts towards the left until the lift histories are in phase at Mach numbers in the range 0.005-0.025.  Although the St = 0.33 case has a larger Strouhal number than the St = 0.20 case, it shows a similar trend in lift as the St = 0.11 case, illustrating that the Mach number becomes more pronounced as the reduced frequency increases as opposed to the amplitude increasing.  Similar observations were seen for the drag time histories and are not shown.  The r.m.s. vorticity profiles show that as the Mach number decreases, the maximum r.m.s. vorticity values decrease while the profile begin to collapse.     201  Figure C.3.  Lift time histories (top) and r.m.s vorticity profiles (bottom) measured at x/c = 1.0 for different symmetric pitching cases illustrating solution convergence by Mach number.  The error due to decreasing Mach number was quantified using the RMSD of the lift and drag time histories, maximum and minimum values, and the RMSD of r.m.s. vorticity profiles.  The error calculated as the root-mean-square difference in the lift and drag time histories was computed for successive Mach numbers, M&,1 and M&,2; M&,1 is the lower of the two Mach numbers.  As the Mach number decreases to the range of 0.005-0.015, the error of the differences in the force histories decreases to or below 3%, as shown in Figure C.4.  The error in the maximum and minimum values of lift and drag is shown in Figure C.5, which decreases to or below 3% as the Mach number decreases to the range of 0.005-0.015.  As the Mach number decreases, the average drag values increase (see Table C.1).  The error of the average drag, not shown in the figure, dropped below approximately 4% as the Mach number decreases to the range 0.005-0.015.  However, it should be noted that the average drag values are small and may be susceptible to some numerical uncertainty due to stiffness in the equations as a result of the Mach number decreasing towards zero.  Quantitative convergence of the r.m.s. vorticity profiles  202 measured at different locations is shown in Figure C.6, which shows the error decreases to less than 1.0% of max("z,rms(y)/c/U&) as the Mach number decreases to the range 0.005-0.015.  Figure C.4.  Error percentage for lift and drag histories versus Mach number for sinusoidal pitching cases.  Figure C.5.  Error percentage for lift and drag maximum and minimum values versus Mach number for sinusoidal pitching cases.        203 Table C.1.  Mach number effect on average thrust coefficient for different cases. M! St = 0.11 St = 0.20 St = 0.33 0.005 0.01337 -0.04476 -0.2278 0.015 0.01291 -0.04673 -0.2294 0.025 0.01263 -0.04916 -0.2314 0.050 0.0118 -0.05403 -0.2358 0.075 0.01125 -0.04565 -0.2336 0.100 0.01161 -0.01858 -0.2162   Figure C.6.  Error percentage for r.m.s. vorticity profiles versus Mach number for sinusoidal pitching cases.  The behavior in the unsteady lift as the Mach number decreases Mach number can be seen in the pressure coefficient distribution around the airfoil at different phases, which are shown for k ! 12 in Figure C.7.  At each phase, Mach numbers in the range 0.005-0.025 produce very close pressure distributions, with slight quantitative differences.   Mach numbers in the range 0.050-0.100 show different values of pressure along the surface while there is a notable phase difference as the Mach number increases.  Thus, Mach numbers in the range 0.005-0.025 are considered within the Òincompressible limit.Ó  The selected Mach number for the symmetric pitching simulations was 0.015 to be remain as low as possible while not being as computationally expensive as 0.005.  204  Figure C.7.  Instantaneous surface pressure coefficient distributions for k ! 12 at four phases illustrating effect of Mach number. C.2. Asymmetric Pitching Convergence Studies  A similar exercise as done for the symmetric pitching was performed for the asymmetric trajectory.  The flow conditions were as follows: Re = 12,000, M& = 0.005, k = 6.68, !o = 2¡, and S = 30%.  The grid was coarsened from 5 million grid points to 3.6 million (Medium Grid) and 2.5 million (Coarse Grid) with similar grid point distributions, the time-step was varied between 5x10-5 to 1.25x10-5, the number of subiterations was increased from 3 to 9.  To present convergence, time histories of lift and instantaneous vorticity profiles measured one chord length from the trailing edge, at a phase of 0.50T, are shown in Figure C.8.  In general, are were no visual differences, nor great quantitative differences, between the lift histories for the three studies.  The effect of grid resolution on the vorticity profiles shows the medium and fine grids produce nearly identical profiles while the coarse grid does show some general agreement with the other two grids.  However, the shape of the coarse grid profile result differs from the fine and medium grids in general.   Similar observations were also seen with  205 respect to time-step size.  The two smaller time-steps produce the same vorticity profile while the larger time-step shows general similarities but quantitative differences.  The number of subiterations showed much more sensitivity.  Three subiterations do not produce a converged solution, where as five or more subiterations produce the nearly identical profiles, with an increasing number of subiterations reducing the background noise in the profiles (Figure C.9).  Based on these data, the fine grid, seven subiterations and 1.25x10-5 time-step size were selected.  Figure C.8.  Lift time histories (top) and instantaneous vorticity profiles (bottom) measured at x/c = 1.0 at #/T = 0.50 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M& = 0.005, k = 6.68, and !o = 2¡, and S = 30%.     206  Figure C.9.  Instantaneous vorticity profiles measured at x/c = 1.0 at #/T = 0.25 illustrating solution convergence by grid resolution (left), number of subiterations (middle), and time-step size (right).  Flow conditions: Re = 12,000, M& = 0.005, k = 6.68, and !o = 2¡, and S = 30%. To illustrate the effect of Mach number, three cases were considered for k = 6.68: S = 50%, 38%, and 30%.  The Mach number varied between 0.005 to 0.100.  To present convergence, time histories of lift and fluctuation vorticity profiles measured one chord length from the trailing edge are shown in Figure C.10.  First, the lift histories are considered.  As the Mach number decreases, the peak lift in general decreases until it reasonably converged at M& = 0.005-0.025 for all three cases.  The phase of peak lift also shifted towards the left until the lift histories are in phase, also occurring around at Mach numbers in the range 0.005-0.025.  However, the introduction shows a much greater sensitivity to the Mach number.  Just after the acceleration jump, there is a dynamic response that becomes ÒpeakierÓ as the Mach number decreases.  The time-scale in which this inertial peak occurs decreases as the Mach number decreases.    207  Figure C.10.  Lift time histories (top) and fluctuation vorticity profiles (bottom) measured at x/c = 1.0 for different asymmetric pitching cases illustrating solution convergence by Mach number. Again, the error was quantified using the RMSD of the lift and drag time histories, maximum and minimum values, and the RMSD of r.m.s. vorticity profiles.  The error calculated as the root-mean-square difference in the lift and drag time histories was computed for successive Mach numbers, M&,1 and M&,2; M&,1 is the lower of the two Mach numbers.  As the Mach number decreases to the range of 0.005-0.015, the error of the differences in the force histories decreases to or below 3%, as shown in Figure C.11.  The error in the maximum and minimum values of lift and drag is shown in Figure C.12, which decreases to or below 3% as the Mach number decreases to the range of 0.005-0.015.  As the Mach number decreases, the average lift decreases in magnitude while the average drag increases (see Table C.2 and C.3).  The error of the average lift and drag, not shown in the figure, drops below approximately 3% as the Mach number decreased to the range 0.005-0.015, except for the drag at S = 30%.  However, the drag values are on the order of 10-3.  Quantitative convergence of the r.m.s. vorticity profiles measured at different locations is shown in Figure C.13, which shows the error decreases to less than 3.0% of max("z,rms(y)/c/U&) as the Mach number decreases to the range 0.005-0.015.  208  Figure C.11.  Error percentage for lift and drag histories versus Mach number for sinusoidal pitching cases.  Figure C.12.  Error percentage for lift and drag maximum and minimum values by Mach number for asymmetric pitching cases.           209 Table C.2.  Average lift coefficient for different asymmetries. M! S = 38% S = 30% 0.005 -0.01950 -0.01633 0.015 -0.02037 -0.01670 0.025 -0.02155 -0.01972 0.050 -0.02771 -0.02970 0.075 -0.02919 N/A 0.100 -0.02412 N/A  Table C.3.  Average thrust coefficient for different asymmetries. M! S = 50% S = 38% S = 30% 0.005 0.01337 0.01104 0.0064 0.015 0.01291 0.01051 0.0057 0.025 0.01263 0.01012 0.0051 0.050 0.0118 0.00910 0.0038 0.075 0.01125 0.00908 N/A 0.100 0.01161 0.01030 N/A   Figure C.13.  Error percentage for r.m.s. vorticity profiles versus Mach number for sinusoidal pitching cases.     210 APPENDIX D.  ADDITIONAL VALIDATION D.1. Additional Symmetric Pitching Validation Average and r.m.s. velocity profiles measured one chord downstream of the trailing edge from the current computations are compared with velocity profiles from experiments by Bohl and Koochesfahani (2009) in Figure D.1 at a reduced frequency of approximately 12.  The experimental profiles are shifted vertically so that the maximum u velocity is at y/c = 0.  Overall, the velocity profiles show good agreement.  The average and r.m.s. velocity maximum value between the computation and the experiment is extremely close, primarily in the average streamwise velocity.  Both the experiment and computation also produce three peaks in the urms profile, though the peaks produced by the CFD are slightly higher than those from the experiment.  There is slight discrepancy in the vrms profile, where the CFD shows a clear double peak while the profile from the experiment is relatively flat.   Figure D.1.  Average u (left), r.m.s. u (middle), and r.m.s. v (right) profiles vs. y/c for Reynolds number of approximately 12,000 and 2¡ amplitude.  The reduced frequency is approximately 12.  Data from Bohl and Koochesfahani (2009) is included for comparison.  Profiles are measured at x/c = 1.0 for both the computational and experimental results. A final validation test of the wake resolution for sinusoidal pitching was done by comparing the decay of peak vorticity with the theoretical prediction.  According to Cohn (1999)  211 and Bohl and Koochesfahani (2009), the decay of peak vorticity originates from two mechanisms: viscous diffusion and vortex stretching.  It can be shown that the peak vorticity decay of a single Gaussian vortex in incompressible, viscous flow is  NO.PQR2(NO.SD!*TUV,!,SWX (D.1) where 'z,o is the initial condition of the peak vorticity at location xo, V is the vortex circulation, WX is the vortex convection speed, and ( is the kinematic viscosity,.  The downstream decay of peak vorticity of the positive vortex for a reduced frequency of approximately 12 is presented in Figure D.2.   The decay of peak vorticity from the computation agree very well with theory.    Figure D.2.  Comparison of the downstream decay of peak vorticity between computation and theory.  The initial condition location for the theory is highlighted by the yellow circle.  Flow conditions: Re = 12,000, M& = 0.015, k = 11.9, and !o = 2¡, and S = 50%.     212  D.2. Additional Asymmetric Pitching Validation As shown in Chapter 2, there were quantitative differences in the velocity profiles between the experiment by Koochesfahani (1989) and the computational results for the asymmetrically pitching airfoil with S = 38%.  The differences in the profiles were a result of the properties of the vortices.  A comparison of the centroid location, circulation, core radius, and distance relative to the positive vortex ($xP-N, see Figure D.3) is provided in Tables D.1-D.3, where the experimental data was provided in Naguib, et al., 2011.  The data shows that P1 and N2 are comparable between the experiment (Koochesfahani, 1989; Naguib, et al., 2011) and computations.  The greatest discrepancy between the experiments and computations was with respect to vortex N1.   Table D.1.  Vortex properties for P1 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%.   Property Current Naguib, et al., (2011) yp -0.058 -0.066  %/cU& 0.223 0.23 rc/c 0.0277 0.0287 Table D.2.  Vortex properties for N1 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%.   Property Current Naguib, et al., (2011) yp 0.024 0.08  %/cU& -0.09 -0.12 rc/c 0.023 0.0287 $xP1-N1/a 0.517 0.38  Table D.3.  Vortex properties for N2 comparison between current computations and Naguib, et al., (2011) for k = 6.68 and S = 38%.   Property Current Naguib, et al., (2011) yp -0.18 -0.20  %/cU& 0.11 0. 11 rc/c 0.0268 0.0263 $xP1-N2/a 0.8 0.82  213   Figure D.3.  Schematic of vortex labeling for both S = 50% and S = 38% for vortex properties. In order to isolate potential causes for the disagreement, simulations were performed using the actual pitching trajectory of the experiment in order to remove discrepancy between using analytical formula and the actual trajectory of the experiment.  Due to noise in the experimental trajectory, the pitching velocity was first computed with a second-order central difference scheme and smoothed using a moving average filter.  The filter window consisted of 200 points, or 1% of the period for a reduced frequency of 6.68.  The smoothed pitching velocity was then integrated to produce the smoothed angle of attack history (see Figure D.4).  Figure D.4.  Comparison between the experimental asymmetric trajectory and smoothed trajectory for S = 38%.  214  Uncertainty with the flow conditions of the experiment (Koochesfahani, 1989) could affect how well the computation can reproduce the experiment.  Although the chord length of the airfoil was known, the freestream speed in the experiment was approximately 15.0 cm/s Koochesfahani (1989).  Uncertainty in the value of the freestream velocity impacts both the Reynolds number and reduced frequency of the experiment.  To study the effect of uncertainty in the freestream velocity, three Reynolds numbers (11,000, 12,000, and 13,000) and three reduced frequencies (6.20, 6.40, and 6.68) were chosen.   It was found that the velocity profiles were more sensitive to uncertainty in the reduced frequency was more sensitive than uncertainty in the Reynolds number.  For the nominal Reynolds number of 12,000 as a representative case, the average and r.m.s. velocity profiles measured one chord length downstream of the trailing edge for different reduced frequencies are illustrated in Figure D.5 and compared with the experiment (Koochesfahani, 1989).  As the reduced frequency decreases from 6.68 to 6.2, the average velocity deficit widens while the average velocity surplus narrowed.  The magnitude of the deficit also decreases with decreasing the reduced frequency.  The location of the peaks and troughs in the streamwise velocity fluctuation change as a result of changing reduced frequency.  Part of the differences in the profiles is a result in changes in the trajectories of the three vortices, as shown in Figure D.6, while the other part is related to the vortex properties (not shown).     215  Figure D.5.  Average and r.m.s. velocity profile measured at x/c = 1.0 using actual experimental pitching trajectory for Re = 12,000 at S ! 38%., illustrating effect of k.  Figure D.6.  Wake vortex trajectories for k = 6.20 (left), 6.40 (middle), and 6.68 (right) using actual experimental pitching trajectory for Re = 12,000 at S ! 38%.  There are additional factors that could also contribute to the differences.  For example, the airfoil in the experiment may not have been a perfect NACA 0012 due to machining errors.  Potential uncertainty in both the actual mean angle of attack and pitching amplitude could also influence the vortical structure.  Three-dimensionality caused by axial flow in the vortex core present in the experiment was not captured by the current, two-dimensional computations.    216 Finally, possible flow interactions with test section walls or blockage effects were not modeled in these simulations.  Regardless, the high-order computations accurately capture very similar, albeit not identical, flow features observed in the experiment by Koochesfahani (1989).  Additional tuning of the various flow and kinematic parameters (i.e., Re, k, xp/c, )o, )m, etc.) to achieve exact matching with the experimental profiles would be an exhaustive exercise with little added benefit to this work and thus was not pursued further.    217 APPENDIX E.  SURFACE PRESSURE AND SHEAR STRESS DISTRIBUTIONS The dependence of thrust coefficient on the Reynolds number was a direct result of a Reynolds number effect on the average pressure and shear stress distributions. The surface distributions are non-dimensionalized as the pressure and skin friction coefficient, defined by equations (E.1) and (E.2), respectively.  The non-dimensional freestream pressure (YZ[\ZWZ](D[^_Z]) has been removed from the pressure distribution.  The skin friction coefficient was determined by computing the shear stress at the wall.  The wall shear stress was calculated by `aRbb(c-d-5aRbb, where µ is the dynamic viscosity of the fluid and U is the contravariant velocity tangent to the wall (recall W(efg@ehi@eOj).     kP(lY!YZ\ZWZ] (E.1)  k'(l`aRbb\ZWZ] (E.2) First, consideration is made with respect to the static airfoil at zero angle of attack.   The effect of Reynolds number on the average pressure and skin friction coefficient distributions is shown in Figure E.1. As the Reynolds number increases, there is an increasing amount on suction near the leading edge.  There is also a decreasing amount of suction past X/c = 0.4 as the Reynolds number increased.  However, the pressure distribution shows no indication of approaching the inviscid limit near the trailing edge.  The shear stress for each Reynolds number shows a peak near the leading edge, which decreases near zero towards the trailing edge.  The  218 maximum skin friction coefficient non-linearly decreases in magnitude from 0.22 to 0.0735 as the Reynolds number increased from 2,000 to 22,000.   Figure E.1.  Average pressure coefficient (top) and skin friction coefficient (bottom) for static NACA 0012 at ! = 0¡ for different Reynolds numbers. Now, the pitching airfoil is considered.  Three examples of the pressure and shear stress distribution for Re = 2,000 and 22,000 are shown in Figure E.2 for 2¡ amplitude.  Three reduced frequencies are considered, 0, 6, and approximately 12.  As seen in the figure, the pressure distributions show no drastic effect of k near the leading edge.  Over the majority of the airfoil, there is a slight increase in suction as the reduced frequency increases for the two Reynolds numbers.  Near the trailing edge, however, there is a very strong increase in suction as k increased.  The region near the trailing edge over which the suction occurs also increases as k increases.  Similar results are seen in the shear stress.   From the leading edge to the mid-chord, there is negligible effect of k for the two Reynolds number.  Beyond the mid-chord, the shear stress increases, particularly near the trailing edge.   219  Figure E.2.  Average pressure and skin friction coefficient for a Reynolds number of 2,000 (left) and 22,000 (right) oscillating at different k.  The pitching amplitude is 2¡.  The airfoil is shown at the bottom for spatial reference. By projecting the pressure and shear stress distributions in the streamwise direction, the average drag/thrust could be computed.  The projected pressure coefficient, (nxCp)avg, and skin friction coefficient, (nxCf)avg, for the two Reynolds number and different reduced frequencies are shown in Figure E.3.  If the quantity (nxCp)avg is positive, it is called projected pressure and contributes to thrust when integrated.  If the quantity (nxCp)avg is negative, it is called projected suction and contributes drag.  Conversely, positive (nxCf)avg contributes to drag and negative (nxCf)avg contributes to thrust when integrated.  The projected pressure and skin friction distributions for reduced frequencies of 0, 6, and approximately 12 are plotted in Figure E.3.    220  Figure E.3.  Average projected pressure and skin friction coefficient for a Reynolds number of 2,000 (left) and 22,000 (right) oscillating at different k.  The pitching amplitude is 2¡.  The airfoil is shown at the bottom for spatial reference. The projected pressure between the two Reynolds does not differ significantly for k = 0, except near the leading edge where the (nxCp)avg experiences a peak and between the mid-chord and trailing edge.  The projected friction decreases in magnitude as the Reynolds number increases to mn for the static case as well.  As the reduced frequency increases from k = 0, the peak projected pressure near the leading edge does not significantly increase for either Reynolds number while the projected suction between location of maximum thickness and the trailing edge switches to projected pressure, and thus leading to a contribution of positive CT, by a reduced frequency of 6.  Between the mid-chord and trailing edge, this projected pressure region experiences a peak near X/c = 0.75 for both Reynolds numbers.  The projected pressure drops  221 dramatically downstream of X/c = 0.9 and switches to projected suction.  The peak projected friction near the leading edge does not visibly increase at this k for either Reynolds number, while there is an increase in friction past the peak location.  Near the trailing edge, the projected friction increases dramatically for both Reynolds numbers compared to the static case.  At k ! 12, the peak projected pressure increases slightly for both Reynolds numbers while the secondary peak between the mid-chord and trailing edge increases dramatically to values close to the primary peak.   This secondary region of projected pressure extends over a greater portion of the airfoil past the mid-chord than at k = 6, leading to high thrust.  The projected friction peak value near the leading edge again does not visibly change at this k.  However, the increase in friction beyond the peak increases further for both Reynolds numbers, particularly near the trailing edge.      222 APPENDIX F.  LINEAR THEORY Early linearized theory for a flat plate oscillating in an inviscid, incompressible flow (Glauert, 1929; von K⁄rm⁄n and Burgers, 1935; Theodorsen, 1935; Garrick, 1936) has shown that the lift and drag vary with time as a result of a time-varying angle of attack.  For pure pitching, the time-varying angle of attack is defined in complex notation as o(oSn6pq, where H(!D, )o is the pitching amplitude in radians, ' is the angular frequency in units of rad/s, and t is the physical time.  The angular frequency is related to the non-dimensional reduced frequency by k = 'c/2U&.  Although the closed form solutions unsteady lift and drag were originally derived by Glauert (1929) and von K⁄rm⁄n and Burgers (1935), Theodorsen (1935) and Garrick (1936) independently derived the lift and drag, respectively, and are more traditionally referenced.   It is assumed that the geometry is a flat plate while the flow is inviscid and irrotational such that the leaves the trailing edge smoothly (i.e., the Kutta condition).  A schematic of this is provided in Figure F.1.  This would require the flow to remain attached along the entire surface with no flow reversal.  The theory allows for the plate to undergo both plunging and pitching motions of small amplitude, where the non-dimensional pivot point m is related to the physical pitch-axis location as m = 2xp/c Ð 1.  Pitching about the physical quarter-chord yields m = -0.5.  As the airfoil oscillates, vorticity continuously sheds from the trailing edge in accordance with the Kelvin condition and is carried downstream by the fluid at freestream speed.  This shed vorticity remains planar along the wake centerline.   223  Figure F.1.  Schematic of oscillating plate (Theodorsen, 1935; Garrick, 1936). Employing velocity linearization, g(WZ@rCrf (where s is the velocity potential), the lift is given by equation (F.1).  Observation of the terms in the equation shows that the unsteady lift depends only on the freestream conditions, chord length, and plate kinematics.  To enforce the Kutta condition, Theodorsen introduced a complex attenuation function, C (equation (F.3)).  This function is written in terms of Bessel functions of the first and second kind (J0, J1, Y0, and Y1) that depend only on the reduced frequency.  The behavior of the components F and G to TheodorsenÕs function with the inverse of the reduced frequency can be seen in Figure F.2.  Garrick (1936) gave the unsteady drag force as equation (F.3)  tu(T\ZWZvWZo@vlDl!wokxyz{|}~•z†@T\Zv]*WZo!wvlo‡•…9xyz{|}~•z† (F.1) where k(+@H—(–:–:@ƒS@ƒ:ƒ:!–S–:@ƒS]@ƒ:!–S]!Hƒ:ƒS@–:–S–:@ƒS]@ƒ:!–S] (F.2)   224  Figure F.2.  F and G vs. 1/k (Theodorsen, 1935; Garrick, 1936).  ⁄u(!T\ZvlkWZo@vlDl!wo!vlo]‹|y•…=›•z{−!otu‰z•„−{−“=›•z{− (F.3) Examination of the equations (F.1) and (F.3) shows that the unsteady forces consist of two terms each.  The lift force is generated by a contribution related to the changing circulation around the airfoil and a contribution caused by the inertia of the fluid displaced by the plate.  These two terms are called the circulatory component and non-circulatory component, respectively.  The phasing behavior of the unsteady lift primarily falls one of three characteristic categories.  At low reduced frequency, the lift is primarily in phase with the angle of attack.  At kÕs around 1, the lift is is approximately in phase with the pitching velocity.  At higher reduced frequencies, the lift is dominated by the pitching acceleration.  The drag force can be decomposed into a term related to a suction force at the leading edge caused by vorticity being infinite at the leading edge (von K⁄rm⁄n and Burgers, 1935) and a term related to the projection of the lift created by the pitching motion.  The lift and drag forces are linearly and quadratically dependent on the pitching amplitude, respectively.  Therefore, the lift and drag scale with oS and oS], respectively.   225 Garrick (1936; 1957) also derived a formula for the average thrust on a flat plate pitching in an inviscid, irrotational flow and is written in equation (F.4) in non-dimensional coefficient form.    k”(T‘]oS]+]@—]D‘]@Dl!w]@Dl!+Dl!w!+‘]!Dl!w—‘ (F.4) The average thrust is negative if the fluid resists the airfoil (i.e., drag) and is positive if it propels the airfoil in the negative streamwise direction. In this expression, the thrust depends quadratically on pitching amplitude and non-linearly on the reduced frequency and non-dimensional pitch-axis location, m.  The reduced frequency in which the drag switches to thrust (thrust crossover reduced frequency) based on this formula occurs at a reduced frequency of 1 and is independent of amplitude.  Note that an incorrect formula was given by GarrickÕs 1936 publication (see equation (F.5) while the corrected formula given in equation (F.4) was presented in 1957.    k”(T‘]oS]+]@—]D‘]@Dl!w]@DlDl!w!+‘]!Dl!w—‘ (F.5) Finally, Garrick (1936) provided the propulsive efficiency for a flat plate pitching as inviscid, irrotational flow as  ’‚(+]@—]D‘]@Dl!w]@Dl!+Dl!w!+‘]!Dl!w—‘DlDl!w!Dl!w+Dl!w@—‘Dl!+ (F.6) According to equation (F.6), the propulsive efficiency is independent of amplitude and depends only on reduced frequency and pitch-axis location.    226 APPENDIX G.  ADDITIONAL VORTEX PROPERTY RESULTS G.1. Vortex Properties for 2¡ Amplitude Data The vortex properties as a function of reduced frequency for different Reynolds numbers pitching with 2¡ amplitude are shown in Figure G.1.  The transverse spacing increases from a negative value and switches to a positive value as the reduced frequency increases.  At a given k, the spacing increases as the Reynolds number increases until the reach a nearly asymptotic value of 0.060c, which was independent of Re.  For reference, the trailing edge amplitude was 0.052c.  The streamwise spacing decreases with reduced frequency and increases with Re.  In other words, the vortices bunched together more as the Reynolds number decreases.  The cause for the increase in a by Re can be seen in the convection speed, which increases with increasing Re.  Physically, the vortices at the low Re convect more slowly than at higher Re.   The peak vorticity increases with linearly reduced frequency and non-linearly Reynolds number.  As the Reynolds number increases from 2,000 to 22,000, the peak vorticity increases by a factor of 4.5 at k = 5.2 and 4.2 at k ! 12.  The core radius decreases as reduced frequency and Reynolds number increases, though this was expected.  For each Reynolds number, the core radius decreases as the reduced frequency increases.  The vortex circulation, interestingly, remained relatively unchanged (to within the symbol size) among the different Reynolds numbers over the range of reduced frequencies.    227  Figure G.1.  Vortex properties as a function of reduced frequency for different Reynolds numbers.  The pitching amplitude is 2¡. G.2. Effect of Amplitude  The effect of amplitude on the vortex properties is shown for all four Reynolds numbers in Figures G.2-G.5.  As the pitching amplitude increases, both Reynolds numbers show increase that the peak vorticity, spacing, core radius, spacing, circulation, and convection speed.   The data could also be replotted versus Strouhal number, since it has been shown that the thrust coefficient data shows amplitude collapse with respect to Strouhal number.  This is done for each Reynolds numbers in Figures G.6-G.9.  When the data is rescaled using the Strouhal number, only the circulation collapses and does not occur until Re = 7,000-12,000.  Instead, the peak vorticity and convection speed decrease whereas the transverse spacing, streamwise spacing, and core radius increase as the amplitude increases at a given Strouhal number.    228  Figure G.2.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 2,000.  Figure G.3.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 7,000.  229  Figure G.4.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 12,000.  Figure G.5.  Vortex properties as a function of reduced frequency for different amplitudes.  The Reynolds number is 22,000.  230  Figure G.6.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 2,000.  Figure G.7.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 7,000.  231  Figure G.8.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 12,000.  Figure G.9.  Vortex properties as a function of Strouhal number for different amplitudes.  The Reynolds number is 22,000.  232 APPENDIX H.  ADDITIONAL PITCHING AMPLITUDE FLOW VISUALIZATION RESULTS  Figure H.1.  Instantaneous spanwise vorticity field ('zc/U&) at k = 2.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.  Figure H.2.  Instantaneous spanwise vorticity field ('zc/U&) at k = 2.67 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.  233   Figure H.3.  Instantaneous spanwise vorticity field ('zc/U&) at k = 3.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.  Figure H.4.  Instantaneous spanwise vorticity field ('zc/U&) at k = 3.3 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.  234  Figure H.5.  Instantaneous spanwise vorticity field ('zc/U&) at k = 4.0 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.  Figure H.6.  Instantaneous spanwise vorticity field ('zc/U&) at k = 5.2 for different Reynolds numbers (increasing left to right) and pitching amplitudes (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  The field of view extends over the range 0 " x/c " 3.0 and -0.4 " y/c " 0.4.    235 APPENDIX I.  ADDITIONAL ASYMMTRIC MOTION TRAJECTORY FLOW VISUALIZATION RESULTS  Figure I.1.  Instantaneous spanwise vorticity field ('zc/U&) for k = 4.0 at Re = 12,000 and )o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).   236  Figure I.2.  Instantaneous spanwise vorticity field ('zc/U&) for k = 5.20 at Re = 12,000 and )o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  237  Figure I.3.  Instantaneous spanwise vorticity field ('zc/U&) for k = 5.80 at Re = 12,000 and )o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  238  Figure I.4.  Instantaneous spanwise vorticity field ('zc/U&) for k = 6.68 at Re = 12,000 and )o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).     239  Figure I.5.  Instantaneous spanwise vorticity field ('zc/U&) for k = 8.00 at Re = 12,000 and )o = 2¡ for different asymmetries (increasing from the top down).  The airfoil is at zero angle of attack and is pitching up (*/T = 0.0).  240  Figure I.6.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 5.20, S = 50%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  241  Figure I.7.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 5.20, S = 50%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  242  Figure I.8.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 5.20, S = 38%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  243  Figure I.9.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 5.20, S = 38%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  244  Figure I.10.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-up for k = 5.20, S = 30%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  245  Figure I.11.  Instantaneous spanwise vorticity field (!zc/U") at 12 phases of the pitch-down for k = 5.20, S = 30%, Re = 12,000, and #o = 2¡.  The angle of attack, pitching velocity, and pitching acceleration history is shown at the top with circles indicating the phases at which the vorticity is shown.  The white line is the wake centerline.  The range covered in the figure is 0.8 ! X/c ! 1.2 and -0.075 ! Y/c ! 0.075.  246      BIBLIOGRAPHY    247 BIBLIOGRAPHY   Alpert, P., ÒImplicit Filtering in Conjunction with Explicit Filtering,Ó Journal of Computational Physics, Vol. 44, pp. 212-219, 1981. Anderson, J., Streitlien, K., Barret, D., Triantafyllou, M., ÒOscillating Foils at High Propulsive Efficiency,Ó Journal of Fluid Mechanics, Vol. 360, pp. 41-72, 1998. Anderson, Jr. J., Modern Compressible Flow with Historical Perspective, McGraw-Hill Book Company, 1982. Anderson, Jr. J., Fundamentals of Aerodynamics, McGraw-Hill Book Company, 1984. Ashraf, M., Young, J., Lai, C., ÒReynolds Number, Thickness and Camber Effects on Flapping Airfoil Propulsion,Ó Journal of Fluids and Structures, Vol. 27, pp. 145-160, 2011. Beam, R., Warming, R., "An Implicit Factored Scheme for the Compressible Navier- Stokes Equations," AIAA Journal, Vol. 16, No. 4, pp. 393-402, 1978. Benek, J., Steger, J., Dougherty, F., ÒA Flexible Grid Embedding Technique with Application to the Euler Equations,Ó AIAA 6th Computational Fluid Dynamics Conference, AIAA Paper 1983-1944, 1983. Betz, A., ÒEin Beitrag zur Erkl−rung des SegelfugesÓ, Zeitschrift f¬ur Flugtechnik und Motorluftschiffahrt, Vol. 3, pp. 269-272, 1912. Blasius, H., ÒGrenzschichten in Flsigkeiten mit kleiner reibung,Ó Zeitschrift fr Mathematik und Physik, Vol. 56, pp. 1-37, 1908. Bohl, D., ÒExperimental Study of the 2-D and 3-D Structure of a Concentrated Line Vortex Array,Ó PhD Thesis, Michigan State University, East Lansing, MI, 2002. Bohl, D., Koochesfahani, M., ÒMTV Measurements of the Vortical Field in the Wake of an Airfoil Oscillating at High Reduced Frequency,Ó Journal of Fluid Mechanics, Vol. 620, pp. 63-88, 2009.  248  Borazjani, I., Sotiropoulous, F., ÒNumerical Investgation of the Hydrodynamics of Carangiform Swimming in the Transitional and Inertial Flow Regimes,Ó The Journal of Experimental Biology, Vol. 211, pp. 1541-1588, 2008. Bratt, J., ÒFlow Patterns in the Wake of an Oscillating Aerofoil,Ó A.R.C. Technical Report, R. & M. 2773, 1953. Brunton, S., Rowley, C., ÒEmpiral State-Space representations of TheodorsenÕs Lift Model,Ó Journal of Fluids and Structures, Vol. 38, pp. 174-186, 2013. Buchholz, J., Smits, A., ÒThe Wake Structure and Thrust Performance of a Rigid Low-Aspect-Ratio Pitching Plate,Ó Journal of Fluid Mechanics, Vol. 603, pp. 331-365, 2008.  Chan, W., ÒThe Overgrid Interface for Computational Simulations on Overset GridsÓ, AIAA Paper 2002-3188, 32nd AIAA Fluid Dynamics Conference, 24-26 June, 2002, St. Louis, Missouri. Cleaver, D., Wang, Z., Gursul, I., ÒBifurcating Flow of Plunging Airfoils at High Strouhal number,Ó Journal of Fluid Mechanics, Vol. 708, pp. 349-376, 2012. Cohn, R., ÒEffects of Forcing on the Vorticity Field in Confined Wake,Ó PhD. Thesis, Michigan State University, East Lansing, MI, 1999.   DeLaurier, J., Harris, J., ÒExperiment Study of Oscillating-Wing Propulsion,Ó Journal of Aircraft, Vol. 19, No. 50, pp. 368-373, 1982. Devranjan, S., Jalikop, S., Sreenivas, K., ÒStudy of Flapping Flight using Discrete Vortex Method Based Simulations,Ó International Journal of Modern Physics C, Vol. 24, No. 12, 2013. Dong, H., Mittal, R., Najjar, F., ÒWake topology and hydrodynamic performance of low-aspect-ratio flapping foils,Ó Journal of Fluid Mechanics, Vol. 566, pp. 309-343, 2006. Emblemsvag, J., Suzuki, R., Candler, G.  ÒNumerical Simulation of Flapping Micro Air Vehicles,Ó AIAA Paper 2002-3197, 32nd AIAA Fluid Dynamics Conference and Exhibit 24-26 June, 2002, St. Louis, Missouri. Fish, F., ÒComparative Kinematics and Hydrodynamics of Odontocete Cetaceans: Morphological and Ecological Correlates with Swimming Performance,Ó The Journal of Experimental Biology, Vol. 201, pp. 2867-2877, 1998.    249 Freymuth, P., ÒPropulsive Vortical Signature of Plunging and Pitching Airfoils,Ó AIAA Journal, Vol. 26, No. 7, pp. 881-883, 1988. Fuchiwaki, M., Tanaka, K., ÒTwo-Dimensional Structure of the Wake Behind a Pitching Airfoil with Higher Non-Dimensional Pitching Rate,Ó Journal of Visualization, Vol. 4, No. 4, pp. 323-329, 2001. Gaitonde, D., Shang, J., Young, J. ÒPractical Aspects of Higher-Order Numerical Schemes for Wave Propagation Phenomena,Ó International Journal for Numerical Methods in Engineering, Vol. 45, pp. 1859Ð1869, 1999.  Gaitonde, D., Visbal, M., ÒHigh-Order Schemes for Navier-Stokes Equations: Algorithm and Implementation into FDL3DI,Ó Air Force Research Laboratory AFRL-VA- WP-TR-1998-3060, Wright-Patterson Air Force Base, Ohio, 1998. Gaitonde, D., Visbal, M., ÒFurther Development of a Navier-Stokes Solution Procedure Based on Higher-Order Formulas,Ó AIAA Paper 1999-16436, 37th Aerospace Sciences Meeting and Exhibit, 11-14 January 1999, Reno NV. Gaitonde, D., Visbal, M., ÒPad”-type Higher-Order Boundary Filters for Navier-Stokes equationsÓ, AIAA Journal, Vol. 38, No. 11, pp. 2103-2112, 2000. Gaibraith, M., ÒImplicit Large Eddy Simulation of Low-Reynolds-Number Transitional Flow Past the SD7003 Airfoil,Ó M. S. Thesis, University of Cincinnati, Cincinnati, OH, 2009.  Garmann, D., ÒHigh-fidelity simulations of transitional flow over pitching airfoils, M. S. Thesis, University of Cincinnati, Cincinnati, OH, 2010.  Garmann, D., ÒCharacterization of the vortex formation and evolution about a revolving wing using high-fidelity simulation,Ó PhD. Thesis, University of Cincinnati, Cincinnati, OH, 2013.  Garrick, I., ÒPropulsion of a Flapping and Oscillating Airfoil,Ó NACA Technical Report 567, 1936. Garrick, I., ÒNonsteady Wing Characteristics,Ó High Speed Aerodynamics and Jet Propulsion; Aerodynamic Components of Aircraft at High Speeds, Division F., Vol. VII, Eds. Donovan, A., and Lawrence, H., Princeton University Press, Princeton, NJ., pp. 658-793, 1957. Glauert, H., ÒThe Force and Moment on an Oscillating Airfoil,Ó Aerodynamik und Verwandter Gebiete (ed. Giles, A., Hopf, L., von K⁄rm⁄n, T.) Springer-Verlag, pp. 88-94, 1930.    250 Gieseng, J., ÒNonlinear Two-Dimensional Unsteady Potential Flow with Lift,Ó Journal of Aircraft, Vol. 5, No. 2, pp. 135-143, 1968. Godoy-Diana, R., Aider, J., Wesfreid, W. ÒTransitions in the Wake of a Flapping Foil,Ó Physical Review E, Vol. 77, pp. 1-6, 2008. Hover, F., Haugsdal, O., Triantafyllou, M., ÒEffect of Angle of Attack Profiles in Flapping Foil Propulsion,Ó Journal of Fluid Structures, Vol. 19, pp. 37-47, 2004. Huang, R., Lin, C., ÒVortex Shedding and Shear-Layer Instability of Wing at Low-Reynolds Numbers,Ó AIAA Journal, Vol. 33, No. 8, pp. 1398-1403, 1995. Jameson, A., Schmidt, W., Turkel, E., ÒNumerical Solutions of the Euler Equations by Finite Volume Methods Using RungeÐKutta Time Stepping Schemes,Ó AIAA Paper 1981-1259, AIAA 14th Fluid and Plasma Dynamic Conference, June 23Ð25, 1981, Palo Alto, California. Jee, S., Moser, R., ÒConservative Integral Form of the Incompressible Navier-Stokes Equations for a Rapidly Pitching Airfoil.  Journal of Computational Physics, Vol. 231, pp. 6268-6289, 2012. Jones, K., Dohring, C., Platzer, M., ÒExperimental and Computational Investigation of the Knoller-Betz Effect on Plunging Airfoils,Ó AIAA Journal, Vol. 36, No. 7, pp. 1240-1246, 1998. Katz, J., Weihls, D., "Behavior of Vortex Wakes from Oscillating Airfoils," Journal of Aircraft, Vol. 15, No. 12, pp. 861-863, 1978.  Katzmayr, R., ÒEffect of Periodic Changes of Angle of Attack on Behavior of Airfoils,Ó NACA TM-147, 1922. Koochesfahani, M., ÒVortical Patterns in the Wake of an Oscillating Airfoil,Ó AIAA Journal, Vol. 27, No. 9, pp. 1200-1205, 1989.  Knoller, R., ÒDie Gesetze des Luftwiderstandes,Ó Flug-und Motortechnik (Wien), Vol. 3, No. 21, pp. 1-7, 1909. Lai, J., Platzer, M., ÒJet Characteristics of a Plunging Airfoil,Ó AIAA Journal, Vol. 37, No. 12, pp. 1529-1537, 1999. Laitone, E., ÒWind Tunnel Tests of Wings at Reynolds Numbers below 70,000,Ó Experiments in Fluids, Vol. 23, pp. 405-409, 1997.  251 Lauder, G., Tytell, E., ÒHydrodynamics of Undulatory Propulsion,Ó Fish Physiology, Volume 23, pp. 425-468, 2006. Lele, S., ÒCompact Finite Difference Schemes with Spectral-Like Resolution,Ó Journal of Computational Physics, Vol. 103, No. 1, pp. 16Ð42, 1992.  Liang, C., Ou, K., Premasuthan, S., Jameson, A., Wang., Z., ÒHigh-Order Accurate Simulations of Unsteady Flows Past Pitching and Plunging Airfoils,Ó Computers and Fluids, Vol. 40, No. 1, 2011.   Lighthill, J., Mathematical Biofluiddynamics, SIAM, Philadelphia, Pennsylvania, 1975. Lissaman, P., ÒLow Reynolds Number Airfoils,Ó Annual Review of Fluid Mechanics, Vol. 15, pp. 223-239, 1983.  Liu, H., Kawachi, K., ÒA Numerical Study of Undulatory Swimming,Ó Journal of Computational Physics, Vol. 155, pp. 223-247, 1999. Lu, K., Xie, Y., Zhang, D., Lan, B., ÒNumerical Investigations into the Asymmetric Effects on the Aerodynamic Response of a Pitching Airfoil,Ó Journal of Fluid and Structures, Vol. 39, pp. 76-86, 2013 Mackowski, A., Williamson, C., ÒDirect Measurement of Thrust and Efficiency of an Airfoil Undergoing Pure Pitching,Ó Journal of Fluid Mechanics, Vol. 765, 524-543, 2015. Message Passing Interface Forum, ÒMPI: A Message-Passing Interface Standard,Ó Computer Science Department Technical Report CS-94-230, University of Tennessee, Knoxville, TN, Apr. 1994.  Monnier, B., Naguib, A., Koochesfahani, M., ÒInfluence of Structural Flexibility in the Wake Vortex Pattern of Airfoils Undergoing Harmonic Pitch Oscillation,Ó Experiments in Fluids, Vol. 56, No. 4, pp. 1-17, 2015. Morgan, P., Visbal, M., Rizzetta, D., ÒA Parallel Overset Grid High-Order for Large Eddy Simulation,Ó Journal of Scientific Computing, Vol. 29, No. 6, 2005. Naguib, A., Vitek, J., Koochesfahani, M., ÒFinite-Core Vortex Array Model of the Wake of a Periodically Pitching Airfoil,Ó AIAA Journal, Vol. 49, No. 7, pp. 1542-1550, 2011.  252 Olson, D., Naguib, A., Koochesfahani, M., ÒMeasurement of the wall pressure and shear stress distribution using molecular tagging diagnostics,Ó Experiments in Fluids, Vol. 56, No. 171, 2015. Ou, K., Castonguay, P., Jameson, A. Ò3D Flapping Wing Simulation with High Order Spectral Difference Method on Deformable Mesh,Ó AIAA paper 2011-1316, 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition 4 - 7 January 2011, Orlando, Florida. Park, K., RŠsen, M., Hedenstrım, A., ÒFlight Kinematics of the Barn Swallow (Hirundo Rustica) over a Wide Range of Speeds in a Wind Tunnel,Ó The Journal of Experimental Biology, Vol. 204, pp. 2741-2750, 2001. Pletcher, R., Chen, K., ÒOn Solving the Compressible NavierÐStokes Equations for Unsteady Flows at Very Low Mach Numbers,Ó AIAA Paper 93-3368, July 1993. Pointwise, GridGen User Manual Version 15, 1997. Pope, S., Turbulent Flows, Cambridge University Press, 10th ed., 2013. Prandtl, L. ÒTragflltheorie. I.Mitteilung, Nachrichten der Gesellschaft derWissenschaften zu Gıttingen,Ó Mathematisch-Physischen Klasse der, pp. 151Ð177, 1918. Press, W., Vetterling, W., Teukolsky, S., Flannery, B., Numerical Recipes in Fortran: The Art of Scientific Computing, Cambridge University Press, 2nd ed., 1992. Pulliam, T., Steger, J., ÒOn Implicit Finite-Differing Simulations of Three-Dimensional Flows,Ó AIAA Paper 1978-10, 16th Aerospace Sciences Meeting, 16-18 January 1978, Huntsville, AL. Pulliam, T., "Efficient Solution Methods for the Navier-Stokes Equations," Lecture Notes for the von K⁄rm⁄n Institute for Fluid Dynamics Lecture Series: Numerical Techniques for Viscous Flow Computation In Turbomachinery Bladings, von K⁄rm⁄n Institute, Rhode-St-Genese, Belgium, 1985. Pulliam, T., ÒTime Accuracy and the use of Implicit Methods,Ó AIAA Paper 1993- 3360, 11th Computational Fluid Dynamics Conference, July 1993, Orlando, FL. Pulliam, T., Steger J., ÒImplicit Finite-Difference Simulations of Three- Dimensional Compressible Flow,Ó AIAA Journal, Vol. 18, No. 2, pp. 159-167, 1980.  Pulliam, T., Chaussee, D., "A Diagonal Form of an Implicit Approximate- Factorization  253 Algorithm," Journal of Computational Physics, Vol. 39, No. 2, pp. 347-363, 1981 Pulliam, T., ÒArtificial Dissipation Models for the Euler Equations,Ó AIAA Journal, Vol. 24, No. 12, pp. 1931-1940, 1986. Ramamurti, R., Sandberg, W., ÒSimulation of Flow about Flapping Airfoils Using a Finite Element Incompressible Flow Solver,Ó AIAA Journal, Vol. 38, No. 2, pp. 253-260, 2001. Rizzetta, D., Visbal, M., Morgan, P., ÒA High-Order Compact Finite Difference Scheme for Large Eddy Simulation of Active Flow Control (Invited),Ó AIAA Paper 2008-526, 46th AIAA Aerospace Sciences Meeting and Exhibit 7 - 10 January 2008, Reno, Nevada. Ros”n, M., Spedding, G., Hedenstrım, A., ÒThe Relationship Between Wingbeat Kinematics and Vortex Wake of a Thrush Nightingale. Journal of Experimental Biology, Vol. 207, pp. 4255-4268, 2004. Saffman, P., Vortex Dynamics, Cambridge University Press, New York, 1992. Sarkar, S., Venkatraman, K., ÒNumerical Simulation of Incompressible Viscous Flow Past a Heaving Airfoil,Ó International Journal of Numerical Methods in Fluids, Vol. 29, pp. 1-29, 2006. Sarkar, S. Chajjed, S., Krishnan, A., ÒStudy of Asymmetric Hovering in Flapping Flight,Ó European Journal Mechanics B/Fluids, Vol. 37, pp. 72-89, 2013. Schnipper, T., Andersen, A., Bohr, T., ÒVortex Wakes of a Flapping Foil,Ó Journal of Fluid Mechanics, Vol. 633, pp. 411-423, 2009. Sheldahl, R. Klimas, P., ÒAerodynamic Characteristics of Seven Symmetrical Airfoil Sections Through 180o Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines,Ó Tech. Rep. 80Ð2114, Sandia National Labs, Albuquerque, NM, 1981.  Sherer, S., Scott, J., ÒHigh-order compact finite-difference methods on general overset grids,Ó Journal of Computational Physics Vol. 210, pp. 459-496, 2005.  Sherer, S., Visbal M., Galbraith, M., ÒAutomated PreProcessing Tools for Use with a High-Order Overset-Grid Algorithm,Ó AIAA Paper 2006-1147, 44th AIAA Aerospace Sciences Meeting and Exhibit 9 - 12 January 2006, Reno, NV. Sherer, S., Visbal, M., ÒMulti-Resolution Implicit Large Eddy Simulation using a High-Order Overset Grid Approach,Ó International Journal of Numerical Methods in Fluids Vol. 55, pp. 455-482, 2007.  254 Silverstein A., Joyner, U., ÒExperimental Verification of the Theory of Oscillating Airfoils,Ó NACA Technical Report 673, 1939. Stanek, M., Visbal, M., ÒStudy of the Vortical Wake Patterns of an Oscillating Airfoil,Ó AIAA Paper 1989-0554, 27th Aerospace Sciences Meeting, 9-12 January, 1989 Reno, Nevada. Steger, J., ÒImplicit Finite-Difference Simulations of Flow about Arbitrary Two- Dimensional Geometries,Ó AIAA Journal, Vol. 16, No. 7, pp 679-686, 1978. Streitlien, K., Triantafyllou, G., ÒOn Thrust Estimates for Flapping Foils,Ó Journal Fluids Structures, Vol. 12, pp. 47-55, 1998. Suhs, N., Rogers, S., Dietz, W., ÒPegasus 5: An Automated Pre-Processor for Overset-Grid CFD,Ó AIAA Paper 2002-3186, 32nd AIAA Fluid Dynamics Conference 24-26 June 2002, St. Louis, MO.  Tannehill, J., Anderson, D., Pletcher, R., Computational Fluid Dynamics and Heat Transfer, McGraw-Hill Book Company, 1984. Tay, W., Lim, K., ÒAnalysis of Non-Symmetrical Flapping Airfoils,Ó Acta Mech Sinica, Vol. 25, No. 4, pp. 433-450, 2009. Theodorsen, T. ÒGeneral Theory of Aerodynamic Instability and the Mechanism of Flutter,Ó NACA Technical Report 496, 1935.  Thomas, P., Lombard, C., ÒGeometric conservation law and its application to flow computations on moving grids,Ó AIAA Journal, Vol. 17, No. 10, pp. 1030-1037, 1979. Thomson, W. (Lord Kelvin), ÒOn vortex motion,Ó Transactions of the Royal Society of Edinburgh, Vol. 25, pp. 217Ð260, 1869. Triantafyllou, G., Triantafyllou, M., Grosenbaugh, M., ÒOptimal thrust development in oscillating foils with application to fish propulsion,Ó Journal of Fluids and Structures, Vol. 7, pp. 205-224, 1993. Turkel, E., ÒPreconditioned Methods for Solving the Incompressible and Low Speed Compressible EquationsÓ, Journal of Computational Physics, Vol. 72, pp. 277-298, 1987. Vinokur, M., ÒConservation Equations of Gas-Dynamics in Curvilinear Coordinate Systems,Ó Journal of Computational Physics, Vol.14, pp. 105-125, 1974.  255 Vinokur, M., ÒAn Analysis of Finite-Difference and Finite-Volume Formulations of Conservation Laws,Ó Journal of Computational Physics, Vol. 81, pp. 1-52, 1989. Visbal, M., Gaitonde, D., ÒHigh-Order-Accurate Methods for Complex Unsteady Subsonic Flows,Ó AIAA Journal, Vol. 37, No. 10, pp. 1231Ð1239, 1999.  Visbal, M., Gaitonde, D., ÒVery High-Order Spatially Implicit Schemes for Computational Acoustics on Curvilinear Meshes,Ó Journal of Computational Acoustics, Vol. 8, No., 4, pp. 1259-1286, 2001. Visbal, M., Rizzetta, D., ÒLarge-Eddy Simulation on Curvilinear Grids Using Compact Differencing and Filtering Schemes,Ó Journal of Fluids Engineering, Vol. 124, pp. 836Ð847, 2002.  Visbal, M., Gaitonde, D., ÒOn the Use of High-Order Finite Difference Schemes on Curvilinear and Deforming Meshes,Ó Journal of Computational Physics, Vol. 181, pp. 155-185, 2002. Visbal, M., Morgan, P., Rizzetta, D., ÒAn Implicit LES Approach Based on High-Order Compact Differencing and Filtering Schemes,Ó AIAA Paper 2003-4098, 16th AIAA Computational Fluid Dynamics Conference, 23-26 June 2003, Orlando, Florida.  Visbal, M., ÒHigh-Fidelity Simulation of the Transitional Flow Past a Plunging Airfoil,Ó AIAA Journal, Vol. 47, No. 11, pp. 2686-2697, 2009. Viviand, H., ÒConservative forms of Gas Dynamic Equations,Ó Rech. Aerosp., No. 1974-1, pp. 65-68, 1974. Volpe, G., ÒPerformance of Compressible Flow Codes at Low Mach NumbersÓ, AIAA Journal, Vol. 31, No. 1, pp. 49-56, 1993. Von K⁄rm⁄n, T. 1911 Nachr. Ges. Wissenschaft. Gıttingen Math. Phys. Klasse, 509. Von K⁄rm⁄n, T. 1912 Nachr. Ges. Wissenschaft. Gıttingen Math. Phys. Klasse, 547. Von K⁄rm⁄n, T., Burgers, J. ÒGeneral Aerodynamic Theory Ð Perfect Fluids,Ó Division E, Vol. II, Aerodynamic Theory (ed. Durand, W.), Springer-Verlag, 1935. Von K⁄rm⁄n, T., Sears, W. ÒAirfoil theory for non-uniform motion,Ó Journal Aeronautical Sciences, Vol. 5, No. 10, pp. 379-390, 1938. Wang, S., Zhang, X., He, G., Liu, T., ÒA Lift Formula Applied to low-Reynolds number  256 Unsteady Flows,Ó Physics of Fluids, Vol. 25, 2013. Wang, S., Zhang, X., He, G., Liu, T., ÒEvaluation of lift formulas applied low Reynolds number unsteady flows,Ó AIAA Journal, Vol. 53, No. 1, pp. 161-175, 2014. Wu, T. "Hydromechanics of Swimming of Fishes and Cetaceans," Advances in Applied Mechanics, Vol. 11, pp. 1-63, 1971.  Xiao, Q. Liao, W., ÒNumerical Study on Asymmetric Effect on a Pitching Foil,Ó International Journal of Modern Physics C, Vol. 20, No. 10, pp. 1663-1680, 2009. Young, J., Lai, J., ÒOscillation Frequency and Amplitude Effects on the Wake of a Plunging Airfoil,Ó AIAA Journal, Vol. 42, No. 10, pp. 2042-2052, 2004. Young, J., ÒNumerical Simulation of the Unsteady Aerodynamics of Flapping Airfoils,Ó PhD Thesis, University of New South Wales, New South Wales, Australia, 2005. Yu, Y., Tong, B., ÒA Flow Control Mechanism in Wing Flapping with Stroke Asymmetry During Insect Forward Flight,Ó Acta Mechanica Sinica, Vol. 21, pp. 218-227, 2005. Yu, M., Wang, Z., Hu, H. ÒA Numerical Study of Vortex-Dominated Flow Around an Oscillating Airfoil with High-Order Spectral Difference Method,Ó AIAA Paper 2010-726, AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 4 - 7 January, 2010, Orlando, Florida. Yu, M., Wang, Z., Hu, H., ÒHigh-Fidelity Numerical Simulation of Airfoil Thickness and Kinematic Effects on Flapping Airfoil Propulsion,Ó Journal of Fluids and Structures, Vol. 42, pp. 166-186, 2013. Zhang, Y., Zhengyin, Y., Fei, X., ÒComputational Investigation of a Doubly Hinged Flapping-Wing Model,Ó AIAA Journal, Vol. 50, No. 12, pp. 2643-2654, 2012.