

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


DEVELOPMENT AND VALIDATION OF AN UNSTEADY PANEL CODE TO MODEL AIRFOIL AEROMECHANICAL RESPONSE By AARON M. MCCLUNG Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE July 2004 DEVELOPMENT AND VALIDATION OF AN UNSTEADY PANEL CODE TO MODEL AIRFOIL AEROMECHANICAL RESPONSE By AARON M. MCCLUNG Approved as to style and content by: Eric A. Falk, Chair Andrew S. Arena, Member Gary E. Young, Member Dean, Graduate College ii TABLE OF CONTENTS Page LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii NOTATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv CHAPTER 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2. FUNDAMENTALS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2 Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.3 NavierStokes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.4 Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Potential Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Velocity Potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.2 Superposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Theorems And Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Bernoulii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2 Coefficient of Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Angular Velocity, Vorticity, and Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.1 Motion of a Fluid Element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.2 Angular Velocity and Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 iii 2.4.3 Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.4 Kelvin’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3. PANEL CODES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1 Nonlifting Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1 Discretization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Lifting Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1 Kutta Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.2 Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.3 TimeDependent Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.1 Frame of Reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.2 Wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.3 Unsteady Kutta Condtion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.3.1 BasuHancock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.3.2 Ardonceau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.4 Method of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4. CODE DESCRIPTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 Frame of Reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Gust Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2.1 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2.2 AirfoilGust Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.2.1 Determining Gust Element Condition . . . . . . . . . . . . . . . . . . . . . . 34 4.2.2.2 Case One . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.2.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2.2.4 Case Two . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.3 Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2.4 Gust Influence on the Airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.3 Free Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 iv 4.4 Forced Response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5. CODE VERIFICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1 Wagner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Theodorsen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2.3.1 Pure Pitching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2.3.2 Pure Plunging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2.3.3 Combined Pitching and Plunging. . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2.3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Sears Periodic Gust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3.3.1 Modified NoFlow Boundary Condition . . . . . . . . . . . . . . . . . . . . 75 5.3.3.2 Vortex Sheet Gust Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Kussner’s Sharp Edge Gust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.3.1 Transient Panel Code Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.4.3.2 Single and Double Gust Sheets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.5 Free Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.5.1 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.5.2 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 v 6. FORCED RESPONSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7. SUMMARY AND CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2 Extension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.3 Discussion and Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.4 Contributions of Present Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 APPENDICES A. UVPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 A.1 Revision 120 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 B. COMMON FILES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1 Common Variables Declarations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.1 lengths.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.2 airfoil.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.3 calc.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.4 const.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.5 debug.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.6 file.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.7 forces.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 B.1.8 freeresp.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 B.1.9 freevort.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.10 gau.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.11 graph.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.12 iterative.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.13 motion.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.14 param.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.15 phi.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.16 relax.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.17 strengths.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.18 velocities.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.19 wake.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.20 wakepannel.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 vi B.1.21 graph cons.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 C. INPUT FILES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 C.1 Configuration File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 C.2 Airfoil Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 C.3 Motion History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 C.4 Free Stream Vortices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 D. GRAPHICS ROUTINES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.1 Plotting Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.1.1 Compare Data r10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.2 Animation Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 D.2.1 Animate r21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 vii LIST OF FIGURES Figure Page 3.1 Airfoil modeled with a continuous source sheet. . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Airfoil discretized into constant strength source elements. . . . . . . . . . . . . . . . . . . 15 3.3 Constant Strength Panel Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.1 Influence of a Vortex sheet located in the Freestream flow compared to the Freestream influence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Deformation of a vortex sheet approaching the airfoil leading edge. . . . . . . . . . . 32 4.3 Case One: Gust element straddling the airfoil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Case Two: Gust element endpoint convected into the airfoil. . . . . . . . . . . . . . . . . 34 4.5 Airfoil leading edge vs. the airfoil forwardmost node. . . . . . . . . . . . . . . . . . . . . . 36 4.6 Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+1. . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.7 Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+2. . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.8 Gust element split about the upstream stagnation point. . . . . . . . . . . . . . . . . . . . 38 4.9 Interpolation to determine the time of GustAirfoil impact. . . . . . . . . . . . . . . . . . 39 4.10 Gust element convection along the upper airfoil surface. . . . . . . . . . . . . . . . . . . . . 40 4.11 Gust element convection along the upper airfoil surface. . . . . . . . . . . . . . . . . . . . . 42 4.12 Pitching and Plunging Airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.1 Flat plate at time t = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 viii 5.2 Solutions for the approximate Wagner function, Eq. (5.3), and the approximate Kussner function, Eq. (5.19) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.3 Normalized lift for the NACA 0006, 0010, and 0014 airfoils at α0 = 1deg using a normalized time step of 0.005 compared to Eq. (5.3) . . . . . . . . . . . . . 51 5.4 Normalized lift on a NACA 0010 at α0 = 1, 2, and 4 deg using a normalized time step of 0.010 compared to Eq. (5.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.5 Normalized lift on a NACA 0010 at α0 = 2 deg computed using nondimensionalized time steps of 0.005, 0.075, and 0.010 compared to Eq. (5.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.6 Notation used to describe the Theodorsen pitching and plunging flat plate 53 5.7 Cl vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.8 Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.9 Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.10 Cl vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 58 5.11 Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 58 5.12 Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 59 5.13 Cl vs. Time for a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.14 Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 60 ix 5.15 Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 61 5.16 Cl vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 62 5.17 Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 63 5.18 Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 63 5.19 Cl vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . . . . . 64 5.20 Cmle vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . 65 5.21 Cmea vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . 65 5.22 Cl vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.23 Cmle vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . 67 5.24 Cmea vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . 67 5.25 Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.26 Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . 69 5.27 Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . 69 5.28 Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . . . 70 x 5.29 Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . 71 5.30 Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . 71 5.31 Stationary plate of infinitesimal thickness with periodic transverse gust . . . . . . 73 5.32 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.33 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.34 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.35 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.36 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.37 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.38 Cl vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 xi 5.39 Cmle vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.40 Gust sheet circulation per unit length vs. initial x/c location for a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . 81 5.41 Visualization showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . . . . . 82 5.42 Cl vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.43 Cmle vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.44 Gust sheet circulation per unit length vs. initial x/c location for a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . 84 5.45 Visualization at t = 2.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 85 5.46 Stationary plate of infinitesimal thickness with sharp edge transverse gust . . . . . 86 5.47 Transient lift solutions normilized by the corresponding steady state lift for NACA 0008, 0010, and 0012 airfoils oriented at at α0 = 1deg relative to the timeaveraged freestream computes using a normalized time step of 0.005 compared to Eq. (5.19) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.48 Cl for a single gust sheet of strength γ = −0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gust with an amplitude of ¯ w = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 xii 5.49 Cmle for a single gust sheet of strength γ = −0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gust with an amplitude of ¯ w = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.50 Visualization at t = 2.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 92 5.51 Cl for a pair of gust sheets of strength γ = 0.02 and 0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gusts with amplitudes of ¯ w = 0.01 and 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.52 Cmle for a pair of gust sheets of strength γ = 0.02 and 0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gusts with amplitudes of ¯ w = 0.01 and 0.01. . . . . . . . . . . . . . . . . 93 5.53 Visualization at t = 4.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 94 5.54 Pitch history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.55 Plunge history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.56 Cl history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.57 Cmle history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.58 Plunge vs. Pitch for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.59 Determining modal damping for pitch. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.60 Determining modal damping for plunge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.61 Normalized modal damping vs. freestream velocity for the panel code solution compared to the pk method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 xiii 5.62 Normalized modal frequency vs. freestream velocity for the panel code solution compared to the pk method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.1 Pitch history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Plunge history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Cl history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.4 Cmle history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.5 h/c vs. α for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 xiv NOTATION Φ Velocity Potential ρ Density Γ Total circulation strength γ Circulation strength per unit length λ Source strength per unit length Λ Total source strength n Normal Vector V Velocity Vector q Velocity Vector rij Radius from poing j to point i mc.v. Mass of fluid inside the control volume mout Mass flux out of a control volume min Mass flux into a control volume p Momentum a Acceleration F Force ds Differential along a surface f Body force τij Fluid shear stress p Pressure μ Viscocity coefficient δij Kronecker delta function ω Angular velocity of a fluid element ζ Vorticity of a fluid element L Lift Cl Coefficient of Lift Cd Coefficient of Drag Cm Coefficient of Moment c Chord, length of the airfoil section b Semichord = c 2 s = Ut/b, Semichord location fi Body Force τij Fluid Shear Stress Qh = −L, Generalized force along the +zaxis Qα = My, Generalized moment about the elastic axis m Mass xv h Vertical translation of the airfoil, positive for deflection along the zaxis α Angle between the airfoil centerline and the mean freestream flow Iα Mass moment of inertial per unit span about axis x = ba Sα = mbxα, Static mass imbalance per unit span about axis x = ba ωh = Kh/m, uncoupled natural frequency in bending ωα = Kα/Iα, uncoupled natural frequency in torsion Kh Bending spring stiffness Kα Torsional spring stiffness k = ωb/U, Reduced frequency xvi CHAPTER 1 INTRODUCTION Aeroelastic condsiderations affect a wide range of disciplines. With respect to turbomachinery, particularly the area of highcycle fatigue, aerodynamic forcing of internal components due to rotorstator interactions can significantly impact engine lifecycle and maintenance requirements. To better understand the influence of aerodynamic damping, on highcycle fatigue, the influence of aerodynamic damping on forced structural response must first be be examined. As a first step towards this goal, this thesis develops a computational tool through which the influence of aerodynamic damping can be isolated and systematically studied. 1.1 Goals The goal of this thesis is to develop and validate a computation tool which will enable the systematic investigation into wake induced stuctural responce. The computational tool is based loosely on a HessSmith [5] type unsteady panel code written by Ron Hugo [7, 9] which has been modified to include a freestream gust model and an airfoil structural model. By incorporating the capability to model arbitrary freestream gusts into the unsteady panel code, and coupling the panel code with a structural model, the timedomain response of a body due to an arbitrary freestream disturbance can be computed. 1.2 Organization This thesis is presented in five parts. The first part is an overview of the governing fluid dynamic equations and the derivation of velocity potential which governs the inviscid and 1 incompressible flowfield, as well as the derivation and description of related theorems and concepts which are necessary for the formulation of the numeric solution. The second part describes the formulation of two dimensional panel methods in three sections, starting with the formulation to solve the steadystate flowfield about a nonlifting body, adding the Kutta condition to solve the steadystate flowfield about a lifting body, and then accounting for timedependent effects to solve the timedependent flowfield about a liftingbody undergoing arbitrary motion. The third part of the thesis expands on the timedependent panel method formulation by adding a freestream gust model which represents time dependent freestream pertubations using discrete vortex elements, and a two degree of freedom structural model which allows the responce of an arbitrary body due to aerodynamic forcing to be determined. The fourth part compares the developed panel code against classic analytic solutions for unsteady aerodynamics, and the last part demonstrates the application of the developed panel code to a forced responce problem using a solution which couples the freestream gust model with the structural model. 2 CHAPTER 2 FUNDAMENTALS Before panel codes are discussed, this chapter defines several relations and terms used throughout the later discussion. The first section in the present chapter discusses the derivation of basic governing equations for fluid flow. The second section discusses potential flow and applies basic governing equations to the solution of potential flowfields. The last two sections relate terms and definitions used later in this thesis. 2.1 Governing Equations The fundamental equations governing fluid flow are derived here from the relationships between density, momentum, and energy, and their time rates of change inside a controlvolume. 2.1.1 Continuity The continuity equation relates the time rate of change of mass inside a controlvolume to the mass flux through the controlsurface. The integral form of the continuity equation can be derived by beginning with a statement of mass inside a controlvolume, such as mc.v. = c.v. ρ dV (2.1) Based on Eq. (2.1), the time rate of change of mass inside the controlvolume, ∂mc.v./∂t, is given by ∂mc.v. ∂t = ∂ ∂t c.v. ρ dV (2.2) 3 Mass flux through the controlsurface can also be stated as m˙ in − m˙ out = −c.s. ρqini dS (2.3) If mass is conserved, the net mass flux through the controlsurface must equal the time rate of change of the mass within the controlvolume, leading to the integral form of the continuity equation [8]. ∂ ∂t mc.v. = ∂ ∂t c.v. ρ dV = −c.s. ρqini dS = ˙ min − m˙ out (2.4) The divergence theorem states that given a vector qi, the integral of the normal component of qi relative to the controlsurface equals the integral of the gradient of qi inside the corresponding controlvolume. c.s. qini dS = c.v. ∂ ∂xi qi dV (2.5) By applying Eq. (2.5) to the integral form of the conservation equation, Eq. (2.4), the following simplification can be made c.s. ρqini dS = c.v. ∂ ∂xi (ρqi) dV (2.6) Thus, Eq. (2.4) can be reduced to ∂ ∂t c.v. ρ dV + c.v. ∂ ∂xi (ρqi) dV = 0 (2.7) or c.v. ∂ ∂t ρ + ∂ ∂xi (ρqi) dV = 0 (2.8) Since the volume integral in Eq. (2.8) must equal zero for any arbitrary controlvolume, it must also hold that ∂ ∂t ρ + ∂ ∂xi (ρqi) = 0 (2.9) producing the differential form of the continuity equation [8]. 4 2.1.2 Momentum The momentum equation relates the time rate of change of fluid momentum through a controlvolume to the forces acting on the controlvolume. Momentum is a vector quantity, pj , defined by the product of mass and the corresponding velocity vector. pj = mqj (2.10) For a controlvolume, the summation of forces acting on the volume equal the time rate of change of the controlvolume momentum. c.v. Fj = ∂ ∂t (mqj)c.v. (2.11) When Eq. (2.11) is incorporated with the continuity equation, Eq. (2.4) becomes c.v. Fj = ∂ ∂t (mqj)c.v. = ∂ ∂t c.v. ρqj dV + c.s. ρqjqini dS (2.12) Forces acting on the controlvolume may be either body forces, surface forces, or both. c.v. Fbodyj = c.v. ρfj dV (2.13) c.s. Fsurfacej = c.s. τijni dS (2.14) Thus, substituting Eqs. (2.12) –(2.14) into Eq. (2.11) gives the integral form of the momentum equation. ∂ ∂t c.v. ρqj dV + c.s. ρqjqini dS = c.v. ρfj dV + c.s. τijni dS (2.15) Applying the divergence theorem to Eqs. (2.12) and (2.14) allows simplification of Eq. (2.15) c.s. ρqjqini dS = c.v. ∂ ∂xi (ρqjqi) dV (2.16) 5 c.s. τijni dS = c.v. ∂ ∂xi τij dV (2.17) c.v. ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) − ρfj − ∂ ∂xi τij dV = 0 (2.18) Again, since the volume integral must equal zero for any arbitrary controlvolume, it holds that ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) = ρfj + ∂ ∂xi τij (2.19) producing the differential form of the momentum equation. 2.1.3 NavierStokes If the assumption is made that the fluid is Newtonian (i.e. the stress components τij are linearly related to the derivatives ∂qi/∂xj ), then the following substitution has been widely accepted τij = − p + 2 3 μ ∂qk ∂xk δij + μ ∂qi ∂xj + ∂qj ∂xi (2.20) Thus, Eq. (2.20) can be substituted into Eq. (2.19), giving conservative form of the Navier Stokes relation [8]. ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) = ρfi − ∂ ∂xj p + 2 3 μ ∂qk ∂xk + ∂ ∂xi μ ∂qi ∂xj + ∂qj ∂xi (2.21) 2.1.4 Euler Depending on the flow regime, the NavierStokes equations can be simplified. For example, lowspeed flow about a thin airfoil outside of the boundary layer can be assumed to be incompressible, ρ = constant, and inviscid, μ = 0, if the airfoil is at a conservative angle of attack and large Reynolds Numbers. With these two assumptions, Eq. (2.21) simplifies to the Euler equation. ∂ ∂t qj + ∂ ∂xi (qjqi) = fj + 1 ρ ∂p ∂xj (2.22) 6 2.2 Potential Flow The potential flow assumption is of interest here because it describes the flow regime examined in the current investigation. 2.2.1 Velocity Potential If a flowfield can be considered incompressible, then the continuity equation, Eq. (2.9), simplifies to ∂qi/∂xi = 0. If the flowfield is also inviscid, μ = 0, then vorticity in the flowfield must remain constant with respect to time, ∂ζ/∂t = 0. Given these assumptions, a scaler potential function Φ exists that is a solution to the Laplace equation describing the flowfield ∂2 ∂x2j Φ = 0 (2.23) The potential function, Φ, is often denoted the velocity potential because the velocity field is equal to the gradient of Φ. qj = ∂ ∂xj Φ (2.24) Inversely, the potential at any point, P, in the flowfield can be calculated from any arbitrary reference point, P0, by integrating the velocity field along any path between P0 and P Φ(x1, x2, x3) = P P0 x1 dx1 + x2 dx2 + x3 dx3 = P P0 ∂Φ ∂x1 dx1 + ∂Φ ∂x2 dx2 + ∂Φ ∂x3 dx3 (2.25) Note that with the assumptions of irrotationality and incompressibility, the integrand of Eq. (2.25) is an exact differential, and as such the potential is independent of the integration path. [8] 2.2.2 Superposition Because the velocity potential describes the potential flowfield, and is the solution to the Laplace equation, it holds that [1]: 7 1. Any irrotational incompressible flow has a velocity potential and stream function (for twodimensional flow) that both satisfy Laplace’s equation. 2. Conversely, any solution of Laplace’s equation represents the velocity potential or stream function (twodimensional) for an irrotational, incompressible flow. Since the Laplace equation is a secondorder, linear, partial differential equation, it holds that the sum of two or more particular solutions is also a valid solution. Thus, a complex flowfield with a total potential Φ can be modeled as the superposition of multiple potential solutions, Φk, giving Φ = Φk (2.26) 2.2.3 Boundary Conditions Since solving the Laplace equation is a boundary value problem, applying the correct boundary conditions is essential. The two physical phenomon considered here are the noflow boundary condition at the fluidbody interface, and the farfield condition forcing bodyinduced disturbances to decay to zero strength far from the body. There are two types of boundary condition formulations, the “direct” Neumann boundary condition, and the “indirect” Dirichlet boundary condition. The Dirichlet boundary condition is not explained here because it is not employed in this investigation. See References [1] and [8] for a full explanation. The Neumann boundary condition specifies the normal velocity on the fluidbody boundary must equal zero, ∂Φ ∂n = 0 (2.27) and the potential field due to the presence of the body must be negligible in the farfield (r→∞). lim r→∞ (Φbody) = 0 (2.28) 8 2.3 Theorems And Relations 2.3.1 Bernoulii To compute pressure in a potential flow, the relation between potential and velocity, qj = ∂Φ/∂xj , and the assumption of a conservative body force with potential E, fj = −∂E/∂xj , are substituted in to the Euler equation, Eq. (2.22). ∂ ∂xj E + p ρ + qj 2 2 + ∂Φ ∂t = 0 (2.29) Thus, upon spatial integration E + p ρ + qj 2 2 + ∂Φ ∂t = C(t) (2.30) where C(t) is a spatially independent constant over the entire flowfield, but is a function of time. This is the Bernoulli equation [8]. Because the left hand side of Eq. (2.30) is constant over the entire flowfield at a given point in time, pressure and velocity can be compared at different points in the flow if the potential is known. 2.3.2 Coefficient of Pressure The pressure coefficient is a nondimensional parameter relating pressure between two different locations in a flowfield. Cp = p∞ − p 1 2ρqj 2 ∞ (2.31) Using the Bernoulli equation, Eq. (2.30), the pressure difference in Eq. (2.31) becomes p∞ − p = ρ E + qj 2 2 + ∂Φ ∂t ∞ − ρ E + qj 2 2 + ∂Φ ∂t p (2.32) If care is taken with the choice of reference point, denoted by ∞, such that it exists at a location in the farfield not influenced by any bodyinduced disturbances, then the change in potential with time can be neglected at the reference point. 9 ∂Φ∞ ∂t = 0 (2.33) If the reference point is also chosen such that the difference in the body forces is negligible, E∞ = Ep (2.34) then Eq. (2.32) reduces to p∞ − p = ρ qj 2 ∞ 2 − qj 2 p 2 − ∂Φp ∂t (2.35) Dividing the pressure difference by the free stream dynamic pressure, 1 2ρqj 2 ∞, Eq. (2.35) becomes Cp = 1− qj 2 p qj 2 ∞ − 2 qj 2 ∞ ∂Φp ∂t (2.36) 2.4 Angular Velocity, Vorticity, and Circulation 2.4.1 Motion of a Fluid Element Motion of a fluid element is comprised of translation, rotation, and deformation, where each type of motion is usually caused by different phenomena in the flowfield. Translation is caused by a uniform velocity, where all parts of the element move at the same velocity, disallowing deformation and rotation. Rotation and deformation occur when velocity gradients exist across the element, as can be the case when viscous effects are not negligible. 2.4.2 Angular Velocity and Vorticity The angular velocity of a fluid element relates to the element deformation caused by a velocity gradient. Generally these velocity gradients are caused by shear stresses. The angular velocity of a fluid element, ωi, is defined as the curl of the velocity vector, or ωi = −1 2 εijk ∂qj ∂xk (2.37) 10 Another measure of fluid angular velocity, used to simplify several equations, is vorticity, defined as twice the angular velocity. ζi = 2ωi = −εijk ∂qj ∂xk (2.38) 2.4.3 Circulation Circulation, Γ, is a measure of the vorticity in a fluid region, and equals the integral of the vorticity normal to a surface, S. Γ = S ζini dS (2.39) By substituting the definition of vorticity, Eq. (2.38), into Equation (2.39) and using Stokes Theorm [8], S −εijk ∂qj ∂xk ni dS = C qi dxi (2.40) circulation can be defined as Γ = C qi dxi (2.41) 2.4.4 Kelvin’s Theorem Kelvin’s theorem relates the time rate of change of circulation in a potential flow inside a closed region C. Simply stated, it states that the time rate of change of circulation in a closed fluid region must equal zero, DΓ Dt = 0 (2.42) or the total circulation of a closed fluid region is constant with time. In the case of a lifting body, the body carries some bound circulation related to the body lift. If the body is at steady state, then lift and circulation are constant with time, and Eq. (2.42) is satisfied. For a body that is not at steady state, lift and circulation are functions of time. Therfore, to satisfy Eq. (2.42), another source of equal and opposite circulation must 11 exist in the closed region. From physical observations, the additional circulation is known to be confined to a wake behind the body, Γwake, giving DΓ Dt = (Γbody +Γwake) Δt = 0 (2.43) 12 CHAPTER 3 PANEL CODES This chapter describes the solution of twodimensional potential flowfields using the SmithHess panel method [5]. The description starts with the solution of the flowfield about a nonlifting body, incorporates the Kutta condition to account for bound circulation, and then incorporates timedependent effects to solve for timedependent flowfields using the method of Basu and Hancock [3] as modified by Ardonceau [2]. The solution of the inviscid and incompressible flowfield about a nonlifting body represents the fundamental case to which a panel method can be applied. It also provides a starting point to describe the basic implementation of the panel method which will be expanded upon for the later liftingbody and timedependent solutions. 3.1 Nonlifting Body As described earlier, the inviscid and incompressible flowfield about a nonlifting body can be described by a potential field, which is the combination of the body and freestream potentials. Φ = Φbody +Φ∞ (3.1) To model the body potential, a distributed strength source sheet (of strength λ(s)) is placed along the fluidbody interface, s, as illustrated in Figure 3.1. This allows the body potential at an arbitrary point in the flow, P (x1, x2), to be computed in terms of the potential due to the source sheet. 13 Figure 3.1. Airfoil modeled with a continuous source sheet. Φbody(P) = s λ(s) 2π lnr ds (3.2) Correspondingly, of freestream flow is uniform, parallel to the x1axis, and the origin is a reference point where Φ∞(0, 0) = 0, the potential due to the freestream at point P is Φ∞(P) = qj∞xj (3.3) Substituting Eqs. (3.2) and (3.3) into Eq. (3.1), gives the total potential at point P due to both the body and freestream. Φ(P) = s λ(s) 2π lnr ds + qj∞xj (3.4) The only unknown parameter in Eq. (3.4) is the body source distribution, λ(s). However, by applying the noflow Neumann boundary condition from Eq. (2.27), qjnj = nj ∂Φ ∂xj = 0 (3.5) 14 Figure 3.2. Airfoil discretized into constant strength source elements. to the total potential at the fluidbody interface, nj ∂Φ ∂xj = s nj ∂ ∂xj λ(s) 2π lnr ds + qj∞nj = 0 (3.6) the unknown source distribution can be determined. Unfortunately, solving Eq. (3.6) for the source distribution is a nontrivial exercise for all but the simplest geometries. However, by applying geometric simplifications, determining the body source distribution as a function of body geometry and freestream conditions can be reduced to solving a set of linear equations. 3.1.1 Discretization By discretizing the continuous source distribution, shown in Figure 3.1, into a series of straight segments, or panels, as shown in Figure 3.2, Eq. (3.6) may be reduced to a set of dependent linear equations. For this discussion, each panel represents a unique distributed source element having a constant source strength along the length of the element. A further simplification is made in that the noflow boundary condition is not enforced at all locations 15 on the body. Rather, the noflow boundary conditions are applied to a single location, or collocation point, at the midpoint of each panel, as shown in Figure 3.2. By discretizing the body into N panels, numbered clockwise from panel 1 at the lower body trailingedge to panel N at the upper body trailing edge, the potential at the collocation point of any panel, panel α, can be determined as a function of freestream potential, body geometry, and panel strength distribution along the body. In this manner, the potential on panel α due to a source element on panel β and the freestream is Φαβ = λβ 2π β ln rαβ dsβ + qj∞xjα (3.7) The potential on panel α due to the entire body can be calculated using superposition. Thus, the potential on panel α due to the entire body is the sum of the potential due to the N panels on the body. Φα = N β=1 λβ 2π β ln rαβ dsβ + qj∞xjα (3.8) Applying the noflow boundary condition, Eq. (2.27), to Eq. (3.8) gives the normal velocity on panel α due to the body and freestream. qnα = N β=1 λβ 2π β ∂ ∂nβ ln rαβ dsβ + qj∞njα = 0 (3.9) As in Eq. (3.6), the source strengths in Eq. (3.9) are the unknown. However, because the parameters in the integrand of Eq. (3.9) are based strictly on body geometry, the integral can be replaced by a geometric influence coefficient, aαβ, which represents the geometric influence of panel β on panel α. aαβ = 1 2π β ∂ ∂nα ln rαβ dsβ (3.10) Using the influence coefficient method, the noflow normal condition on panel α, given previously in Eq. 3.9 becomes N β=1 (λβ aαβ) = −qj∞njα (3.11) 16 Figure 3.3. Constant Strength Panel Discretization Equation (3.11) is the basis for a set of linear equations relating the unknown panel source strengths λβ to the noflow boundary condition. This system of equations begins with the influence matrix, Aαβ, which is made up of the influence coefficients, aαβ, based only on the body geometry. Aαβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N a21 a22 ... a2N ... ... ... aN1 aN2 ... aNN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.12) The element strengths for each panel are stored in the column vector xβ. xβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.13) 17 Finally, the column vector Bα represents normal velocity components at the collocation point not induced by the body, such as the freestream normal velocity. Bα = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.14) Combined, these matrices and vectors form a system of equations Aαβxβ = Bα, or ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N a21 a22 ... a2N ... ... ... aN1 aN2 ... aNN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.15) the solution of which is trivial, or nonunique, describing the potentialflow about the nonlifting body. In physical terms, the trivial solution does not include the effects of bound circulation about the body, and therefore does not model lift. 3.2 Lifting Body To model the effects of lift and bound circulation about a body, additional constraints must be considered. To model bound circulation on the body, a set of constant strength vortex panels, each of the same strength, are added to the existing source panel discretization. Since each vortex panel has the same strength, only a single variable must be added to the set of linear equations modeling the nonlifting solution. The additional variable necessitates an additional constraint to solve for the vortex panel strength. This additional constraint is provided by the Kutta condition, which is based on observations of physical flow phenomena about a lifting body, or airfoil, with a sharp trailingedge. 18 3.2.1 Kutta Condition The Kutta condition is a means to relate possible potentialflow solutions about a body to observed physical flow characteristics, thereby generating a unique solution for the flowfield. The general definition of the Kutta condition specifies that the flow must detach from the airfoil at the airfoil trailingedge and that the trailingedge has zero loading [8]. The aforementioned potentialflow solution described for the nonlifting body possess a trailingedge singularity, thus the Kutta condition as specified can not be satisfied at the airfoil trailingedge. For a lifting body, a commonly used first approximation is employed in which the zero loading condition is enforced on the panels adjacent to the airfoil trailingedge. For a discretized airfoil, the condition of zero trailingedge loading is approximately satisfied by specifying equal pressure on the airfoil upper and lower trailingedge panels. The unsteady Bernoulli equation, Eq. (2.30), is used to relate the fluid flow on the upper, u, and lower, l, panels, giving p ρ + q2 j 2 + ∂Φ ∂t l = p ρ + q2 j 2 + ∂Φ ∂t u (3.16) For a steadystate flow, the timedependent potential terms can be neglected in Eq. (3.16), and the condition of equal pressure simplifies to the specification of equal flow velocity on the airfoil upper and lower trailingedge panels. qjl = qju (3.17) Thus, Eq. (3.17) provides the additional constraint necessary to solve for the unique panel strengths on the airfoil in a steadystate flow. 3.2.2 Equations Placing vortex panels along the airfoil does not change the noflow boundary condition described in Eq. (2.27), but the potential at panel α due to panel β and the freestream must now include the influence of the discreyized vortex sheet. 19 Φαβ = λβ 2π β ln rαβ dsβ − γ 2π β θαβ dsβ (3.18) Accordingly, the potential at panel α due to the N source panels, N vortex panels, and the freestream influence along the body becomes Φα = N β=1 λβ 2π β ln rαβ dsβ − γ N β=1 1 2π β θαβ dsβ + qj∞xjα (3.19) Hence, the normal velocity on panel α due to the N panels and freestream can be writen in terms of the potential gradient normal to the body at panel α qnα = N β=1 λβ 2π β ∂ ∂nα ln rαβ dsβ − γ N β=1 1 2π β ∂ ∂nα θαβ dsβ + qj∞njα = 0 (3.20) Again, the integrand for the circulatory term in Eq. (3.20) is based solely on body geometry and therefore may be calculated as a geometric influence coefficient, bα, representing the influence of the N discretized vortex panels on panel α. bα = − N β=1 1 2π β ∂ ∂nα θαβ dsβ (3.21) Substituting the influence coefficients, Eq. (3.10) and Eq. (3.21), the noflow boundary condition gives N β=1 (λβ aαβ) + γ bα = −qj∞njα (3.22) Equation (3.22) still provides N equations, but there are now N + 1 variables (N source strengths, λα, and one vortex strength, γ) describing the potential field about the lifting body. The Kutta condition, Eq. (3.17), provides the N + 1’th condition needed to solve the linear system of equations for the source and vortex strengths. Using the noflow boundary condition to simplify the Kutta condition (i.e. all flow on the trailingedge panels must be tangential) the tangential flow velocity on panel α can calculated 20 in terms of the potential gradient along the body. The tangential flow velocity on panel α due to panel β and the freestream is therefore qsαβ = λβ 2π β ∂ ∂sβ ln rαβ dsβ + γ 2π β ∂ ∂sα θαβ dsβ + qj∞sj (3.23) giving a tangential flow velocity on panel α due to all N body panels of qsα = N β=1 λβ 2π β ∂ ∂sα ln rαβ dsβ + γ N β=1 1 2π β ∂ ∂sα θαβ dsβ + qj∞sj (3.24) Examining Eq. (3.24), two new influence coefficients are introduced, cαβ, the tangential flow component along panel α due to source panel β cαβ = 1 2π β ∂ ∂sα ln rαβ dsβ (3.25) and dα, the tangential flow component along panel α due to the N body vortex panels. dα = N β=1 1 2π β ∂ ∂sα θαβ dsβ (3.26) Rewriting the steadystate Kutta condition, Eq. (3.17), in terms of the geometric influence coefficients, qsl = N β=1 (λβ c1β) + γ d1 = N β=1 (λβ cNβ) + γ dN = qsu (3.27) and rearranging to position the terms on the left hand side, n j=1 (λj (cnj − c1j)) + γ (dn − d1) = 0 (3.28) gives the Kutta condition in a suitable form to incorporate into the system of linear equations, Eq. (3.15). 21 Rewriting Eq. (3.15) to include the vortex influence and the Kutta condition, the Aαβ matrix becomes, Aαβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N b1 a21 a22 ... a2N b2 ... ... ... ... aN1 aN2 ... aNN bN (cN1 − c11) (cN2 − c12) ... (cNN − c1N) (dN − d1) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.29) the xβ vector becomes, xj = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN γ ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.30) and the Bα vector becomes, Bi = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN 0 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.31) Combining Eqs. (3.29), through (3.31) gives the linear system of equations which model the flowfield about the lifting body, ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N b1 a21 a22 ... a2N b2 ... ... ... ... aN1 aN2 ... aNN bN (cN1 − c11) (cN2 − c12) ... (cNN − c1N) (dN − d1) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN γ ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ −qi∞ni1 −qi∞ni2 ... −qi∞niN 0 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.32) 22 providing a unique solution to the flowfield which includes the effects of lift and bound circulation in the solution. 3.3 TimeDependent Solutions The nonlifting and lifting body solutions described in Sections 3.1 and 3.2 provide methods to model steadystate flowfields about a body in a uniform freestream flow. If the body is in motion relative to the freestream, or if the freestream includes perturbations about its timeaverage mean, assumptions neglecting timedependent terms are no longer valid and a timedependent solution methodology must be found. The basic formulation of a timedependent solution is similar to that of the lifting body solution; i.e. the body is discretized using source and vortex panel discretization, the noflow boundary condition provides N linear equations, and the Kutta condition provides the one additional relation necessary to formulate a unique liftingbody solution. The difference between the timedependent and steadystate solutions is the application of the Kutta condition and the incorporation of a model to account for the airfoil wake. The following method describes a solution for the timedependent flowfield about an airfoil in motion relative to the influence of an otherwise uniform freestream flow. 3.3.1 Frame of Reference The choice of coordinate system and reference frame determine the complexity of the mathematical model. For this discussion, the noflow boundary condition is calculated in a bodyfixed coordinate system that is allowed to translate and pitch in the global reference frame with velocity qjrel and pitch rate Ω. 3.3.2 Wake In a viscous solution for attached flow over an airfoil, a low energy boundary layer along the airfoil is shed into the freestream flow from the airfoil trailing edge to form the airfoil wake. The wake represents a pertubation of the freestream flow aft of the airfoil due to the 23 flow about the airfoil. The wake is significant because it can have a profound influence on the flow about the airfoil, even as the wake convects with the freestream. Since viscous effects are neglected in a potentialflow solution, an airfoil wake must be modeled in a way representing the effect of shed bound circulation to satisfy the Kelvin theorm. As such, the wake is often called a “time history” because it represents the change in bound circulation on the airfoil with time. Because the inviscid wake represents the change in bound circulation on the airfoil with time, or the shed circulation, it is possible to model the wake using a set of discrete vortex elements. Basu and Hancock [3] use a set of point vortices and a single constant strength vortex panel to model the inviscid timedependent wake. The strength and orientation of the wake vortex panel, or wake panel, play a key roll in satisfying the timedependent Kutta condition (explained in detail later in this chapter). The the strength of the wake panel is dependent on the amount of circulation shed by the airfoil between time steps, and wake panel orientation is determined by the Kutta condition. After the Kutta condition has been satisfied, calculations necessary to determine the timedependent flowfield solution have been performed, and any necessary post solution calculations have been completed, all wake vortices are convected with the local flow in preparation for the next time step. The wake panel is not convected with the local flow, however, rather the wake panel is replaced by a single point vortex of strength equal to the shed circulation from the previous timestep. This new point vortex is then allowed to convect with the local flow. In this manner, the wak panel and point vortices model shed circulation from the airfoil, which in turn can influence the flow about the airfoil. 3.3.3 Unsteady Kutta Condtion The specification of the timedependent Kutta condition is similar to the steadystate specification described in Section 3.2.1. The difference is in the application of the timedependent Kutta condition, which can no longer neglect timedependent terms in the un 24 steady Bernoulli equation relating pressures on the upper and lower airfoil trailingedge panels. The inclusion of timedependent terms means that the timedependent Kutta condition is a quadratic equation which must be solved iteratively. Two applications of the timedependent Kutta condition are described below, the method of Basu and Hancock [3], and the modification of that method by Ardonceau [2] used in the current investigation. 3.3.3.1 BasuHancock Basu and Hancock propose that “...there is no definitive statement of the Kutta condition for a steady airfoil, each mathematical model requiring its own consistent ‘Kutta’ condition to ensure a unique solution...” [3] Based on that statement, the assumption that the flow separates from the airfoil at the airfoil trailingedge, zero loading exists across the shed vorticity at the trailingedge, and zero loading occurs across the trailingedge elements of the airfoil, Basu and Hancock propose the folowing mathematical model for the Kutta condition. This model determines the orientation, θk, length, Δk, and strength, (γw)k, of the wake panel at time tk. Beginning with the unsteady Bernoulli equation applied at the airfoil trailingedge panels p ρ + q2 j 2 + ∂Φ ∂t l = p ρ + q2 j 2 + ∂Φ ∂t u (3.33) and specifying equal pressure at the trailingedge, pl − pu ρ = q2 ju 2 − q2 jl 2 + ∂Φu ∂t − ∂Φl ∂t = 0 (3.34) a quadratic relation develops for the flow velocity on the upper and lower airfoil trailingedge panels. Because the velocity relation is not linear, an iterative solution is necessary to determine the orientation and strenth of the wake panel which satisfies the Kutta condition. 25 Using the Kelvin theorem, Eq (2.42), the rate of change of circulation about the airfoil must be balanced by the rate of change of the shed circulation in its wake, Δk (γw)k Δt = −∂Γ ∂t = −Γk−1 − Γk Δt (3.35) or the change in circulation about the airfoil from tk−1 to tk must be balanced by an equal and opposite circulation about the wake panel. Δk (γw)k = Γk − Γk−1 (3.36) The rate of change of potential across the airfoil trailingedge is related to the rate of change in circulation by ∂ (Φl − Φu) ∂t = ∂Γ ∂t (3.37) Therefore, substituting Eq. (3.37) into Eq. (3.34) relates the upper and lower trailingedge velocities to the rate of change of circulation about the airfoil. q2 ju − q2 jl 2 + Γk − Γk−1 tk − tk−1 = 0 (3.38) Substituting Eq. (3.36) into Eq. (3.38) gives the circulation strength about the wake panel in terms of trailingedge panels velocities and wake panel length. (γw)k = (tk − tk−1) (q2 l − q2u ) 2Δk (3.39) Wake panel orientation is determined by local velocity on the wake panel, neglecting the effect of the wake panel on itself, tan θk = (q1w)k (q2w)k (3.40) and wake panel length is proportional to the magnitude of the local velocity and the time step. Δk = (qjw)k (tk − tk−1) (3.41) 26 3.3.3.2 Ardonceau Ardonceau proposed a modification to Basu and Hancock’s Kutta condition based on experimental studies [2]. The modified solution method is nearly identical to that of Basu and Hancock, but the wake panel geometry is altered. Instead of allowing the wake panel to change both orientation and length, the wake panel orientation is fixed along the bisector between the upper and lower trailingedge panels. θk = θu + θl 2 (3.42) The length of the Ardonceau wake panel then equals the average of the trailingedge panel velocities porportional to the time step. Δk = 1 2 (qju + qjl)k (tk − tk−1) (3.43) The calculation of wake panel strength is the same as Eq. (3.39). 3.3.4 Method of Solution Regardless of the mathematical formulation of the unsteady Kutta condition, the solution methods are the same. As in the steadystate solutions, the N source strengths, one vortex strength, and freestream along the body are related through the noflow boundary condition which gives a system of N linear equations. As outlined above, however, the Kutta condition becomes a quadratic relation in an unsteady flow which must be solved using an iterative techique. The noflow boundary condition for the timedependent solution also includes induced velocity terms due to body motion relative to the freestream and induced velocity terms due to the airfoil wake. Modifying Eq. (3.22) to include the effects of body rotation, qjrotation = Ω× rα (3.44) 27 body translation, qjtranslation = qjrel (3.45) and the influence of the wake panel and point vortices, qjwake = γw bαN+1 + k−1 β=1 Γβ ∂ ∂nα θαβ 2π (3.46) the time dependent noflow relation becomes N β=1 (λβ aαβ) + γ N β=1 bαβ + γw bαN+1 + k−1 β=1 Γβ ∂ ∂nα θαβ 2π + (qj∞ +Ω× rα + qjrel) njα = 0 (3.47) Rearranging to place all nonsource terms on the right hand side gives N β=1 (λβ aαβ) = −γ N β=1 bαβ−γw bαN+1− k−1 β=1 Γβ ∂ ∂nα θαβ 2π −(qj∞ +Ω× rα + qjrel) njα (3.48) Note that Eq. (3.48) is very similar to Eq. (3.11) but with extra terms on the right hand side. Therefore, rewriting Eq. (3.14) to include the new terms of Eq. (3.48) gives Bi = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −γ N β=1 b1β − γw b1N+1 − k−1 β=1 Γβ ∂ ∂n1 θ1β 2π − (qj∞ +Ω× r1 + qjrel) nj1 −γ N β=1 b2β − γw b2N+1 − k−1 β=1 Γβ ∂ ∂n2 θ2β 2π − (qj∞ +Ω× r2 + qjrel) nj2 ... −γ N β=1 bNβ − γw bNN+1 − k−1 β=1 Γβ ∂ ∂nN θNβ 2π − (qj∞ +Ω× rN + qjrel) njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.49) Substituting Eq. (3.49) for Eqs. (3.14) in Eq. (3.15) gives a linear system of equations for the timedependent solution. 28 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1n a21 a22 ... a2n ... ... ... an1 an2 ... ann ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λn ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −γ N β=1 b1β − γw b1N+1 − k−1 β=1 Γβ ∂ ∂n1 θ1β 2π − (qj∞ +Ω× r1 + qjrel) nj1 −γ N β=1 b2β − γw b2N+1 − k−1 β=1 Γβ ∂ ∂n2 θ2β 2π − (qj∞ +Ω× r2 + qjrel) nj2 ... −γ N β=1 bNβ − γw bNN+1 − k−1 β=1 Γβ ∂ ∂nN θNβ 2π − (qj∞ +Ω× rN + qjrel) njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.50) An iterative solution scheme is used to find the unique solution satisfying both the system of N linear equations and one quadratic relation. The iterative Kutta condition assumes initial values for the wake panel orientation, length, and circulation strength at the initialization of the simulation. The inital values are used in the solution of N linear equations to determine the strength of the body source elements. The calculated body source strengths are then used to recalculate the orientation, length, and circulation strength of the wake panel, and the process is repeated until the orientation, length, and circulation strength of the wake panel meet a given convergence criteria. 29 CHAPTER 4 CODE DESCRIPTION To facilitate investigations into the interaction between elastic airfoil response and arbitrary aerodynamic forcing, two components are added to the timedependent panel code described in Section 3.3. The first component is a gust model originally proposed by Basu and Hancock [3] which uses singularity elements to model the influence of a sharp edge gust. The second component is a structural model which, when coupled with the aerodynamic model, determines the airfoil responce to aerodynamic forcing. This chapter describes the gust and structural models, as well as their implementation and integration into the unsteady panel code. 4.1 Frame of Reference The source and vortex elements modeling the airfoil, wake, and freestream perturbations are tracked in a Lagrangian reference frame. The origin of this reference frame is located at the leading edge of the undisturbed airfoil, with the airfoil trailing edge lying on the x1 axis. 4.2 Gust Model The gust model uses singularity elements to induce velocity perturbations about an otherwise uniform freestream flow. Because this investigation is initially interested in the effect of transverse velocity perturbations, or perturbations perpendicular to the timeaveraged freestream flow, the gust is modeled by a set of vortex sheets. For this discussion, a vortex sheet will be defined as a collection of vortex elements, each sharing at least one end point with a neighboring vortex element. Collectively, these vortex elements produce a continuous 30 Figure 4.1. Influence of a Vortex sheet located in the Freestream flow compared to the Freestream influence. vorticity “sheet” that convects with the freestream flow. Each vortex sheet possess a finite amount of bound circulation that remains constant with time, and is initially distributed evenly along the length of the sheet. The influence of a single gust sheet in the freestream flow is shown in Figure 4.1. By placing gust sheets into the flowfield prior to initialization of the simulation, and specifying a timeinvariant total circulation about each gust sheet, Kelvin’s Theorm is implicitly satisfied. Thus, the Kutta condition discussed in Section 3.3 remains valid. 4.2.1 Deformation The key to properly modeling the transverse gust, such that the gust responds to the influence of the airfoil and its wake, lies in modeling gust convection. To allow the gust sheet to deform and react to the influence of the airfoil, wake, and other elements in the flowfield, each gust sheet is discretized into a finite number of panels, or elements. At the end of each computational time step, the gust sheet is convected by propagating the endpoints of each gust element with the local fluid flow. In this manner, the gust sheet is alowed to deform 31 Figure 4.2. Deformation of a vortex sheet approaching the airfoil leading edge. due to local velocity gradients in the flow. Figure 4.2 illustrates the deformation of a gust sheet as it approaches the airfoil leading edge. Because the total bound circulation about each gust element is timeinvariant, the influence of each gust element on the surrounding fluid flow is a function of element length. 4.2.2 AirfoilGust Interaction Each gust sheet initialized upstream of the airfoil eventually encounters the airfoil as it convects with the freestream flow. Since the gust sheet provides aerodynamic forcing for forcedresponse simulations, proper modeling of airfoilgust interaction is a critical aspect of the overall gust model. To properly model gust sheet influence on the airfoil, the gust sheet must propagate around the airfoil, not propagate through the airfoil. Therefore, the continous gust sheet must be “split” when the gust sheet encounters the forwardmost edge of the airfoil, allowing one section of the sheet to convect across the upper airfoil surface and the other to convect 32 Figure 4.3. Case One: Gust element straddling the airfoil. along the lower airfoil surface. Techniques used to determine if and when a gust sheet must be split, and methods used to split each gust sheet are described below. Two distinct cases can arise when a gust sheet reaches the airfoil. Case One involves a gust element straddling the airfoil leading edge. In this case, one endpoint attempts to convect above the airfoil while the other endpoint convects below, as illustrated in Figure 4.3. To maintain the noflow boundary condition, the gust element straddling the airfoil must be split into two separate elements, one ending on the upper airfoil surface and one ending on the lower airfoil surface. The two “new” elements must also posses a combined bound circulation equal to the bound circulation of the original gust element to satisfy Kelvin’s theorem. Case Two involves a gust element, or pair of elements sharing an endpoint, where the endpoint attempts to convect into the airfoil interior, as illustrated in Figure 4.4. To maintain the noflow boundary condition, the erroneous endpoint must be to the either the upper or lower airfoil surface. In this case, no new gust elements are created, and each affected elements retains its original bound circulation. 33 Figure 4.4. Case Two: Gust element endpoint convected into the airfoil. In each case, once the erroneous elements have been relocated to the airfoil surface, the element endpoints lying on the airfoil surface are convected along the airfoil at the surface tangential velocity, until the gust element propagates past the airfoil trailing edge. 4.2.2.1 Determining Gust Element Condition To ensure a gust sheet does not breach the airfoil interior, the following conditions are checked for each gust element after it is convected in preparation for the next time step. 1. Does the gust element currently terminate on the airfoil surface? 2. Is either element endpoint located between the airfoil leading and trailing edges in the x1 direction? 3. Is one element endpoint located above the airfoil while the other is located below the airfoil in the x2direction? 4. Is either gust element endpoint located inside the airfoil surface? 34 Based on the four conditions, the state of the gust element with respect to the airfoil can be determined. If condition 1 is true, the gust element must be convected along the airfoil at the surface tangential velocity instead of the local flow velocity. If condition 1 is false and conditions 2 and 3 are true, the gust element is an example of Case One and must be split. If condition 1 is false and conditions 2 and 4 are true, the gust element is an example of Case Two and the element endpoint must be relocated to the airfoil surface. If conditions 1 through 4 are false, the gust element is located in the freestream and no gust sheet modifications are required. 4.2.2.2 Case One Case One involves splitting a gust element straddling the airfoil, as illustrated in Figure 4.3, and determining the bound circulation about the split gust elements. Because the airfoil may be at some arbitrary orientation relative to the timeaveraged freestream flow, the airfoil leadingedge node may not be the airfoil node the gust sheet first encounters, therefore the term “forward most” edge, or node, will be defined as the node closest to the gust when the gust impacts the airfoil. In addition, depending on airfoil orientation and the influence of the freestream (including the gust), the upstream stagnation point on the airfoil may not correspond to either the forwardmost airfoil node or the airfoil leading edge. This distinction is subtle, as illustrated in Figure 4.5. The variance in x1 location between the leading edge and stagnation point may only be a few hundredths of a chord length, but the influence of this variance on resulting airfoil forcing can be significant. For example, consider the case where a gust element is split about the airfoil leading edge at time tk+1, but the upstream stagnation point on the airfoil does not correspond to the airfoil leading edge. If the upstream stagnation point is located on the lower airfoil surface, the gust element ending on the lower surface between the airfoil leading edge and the upstream stagnation point will convect towards the airfoil leading edge at time tk+2 instead of towards the airfoil trailing edge, as desired. This process is depicted in Figures 4.6 and 35 Figure 4.5. Airfoil leading edge vs. the airfoil forwardmost node. 4.7. In fact, the lower gust element will eventually propagate around the leading edge and convect towards the airfoil trailing edge along the upper surface. This will stretch the gust element through the airfoil, invalidating the noflow boundary condition. A similar circumstance occurs if the gust element is simply split about the upstream stagnation point. For example, if the upstream stagnation point does not correspond to the leading edge, but rather lies on the lower airfoil surface, the gust element propagating above the airfoil will stretch through the airfoil and end on the lower airfoil surface, as illustrated in Figure 4.8. The upper gust element will eventually propagate around the airfoil leading edge before it propagates towards the trailing edge along the upper airfoil surface, as desired, but this gives the gust element an undue influence on the airfoil as it propagates around the airfoil leading edge and will invalidate the noflow boundary condition. Because the upstream stagnation point and the forwardmost node both exhibit large influences on the gust element, gust elements straddling the airfoil leading edge are split about both the forwardmost airfoil node and the upstream stagnation point. In this manner, if the upstream stagnation point is located on the lower airfoil surface, the lower gust element will 36 Figure 4.6. Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+1. Figure 4.7. Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+2. 37 Figure 4.8. Gust element split about the upstream stagnation point. convect along the airfoil from the upstream stagnation point while the upper gust element will convect along the airfoil from the forwardmost airfoil node, and visa versa for an upstream stagnation point located on the upper airfoil surface. In most cases, this distinction is negligible, but the method ensures that split gust elements will not convect towards the airfoil leading edge, or stretch through the airfoil surface. 4.2.2.3 Implementation As mentioned, once it has been determined that a gust element straddles the airfoil, the element must be split into two “new” elements, one convecting above the airfoil and one convecting below the airfoil. Because the unsteady panel code models the flowfield using discrete time steps, it is unlikely that the instant a gust element impacts the forwardmost airfoil node will correspond exactly to a panel code time step. Therefore, an interpolation routine is employed to accurately determine the instant in time a gust element impacts the forwardmost airfoil node, as illustrated in Figure 4.9. It is necessary to know this time because, for example, if a gust element impact occurs at midway between timesteps, the 38 Figure 4.9. Interpolation to determine the time of GustAirfoil impact. element should be convected along the airfoil during the remaining amount of the time step after being split. In addition to accurately determining the instant in time that a gust element impacts the airfoil, the interpolation routine also provides information regarding what percentage of the original gust element should convect above and below the airfoil. Knowing these percentages is necessary so proper fractions of the original bound circulation can be assigned to each “new” gust element, thereby maintaining a constant total circulation in the flow. 4.2.2.4 Case Two Case Two involves a gust element, or pair of elements, possessing an endpoint that convects into the closed airfoil surface, as illustrated in Figure 4.4. This is the less common of the two cases, and for a simulation with a suitably small time step only occurs if the initial gust sheet contains an element possessing an endpoint close to the x1 axis. Therfore, in an effort to simplify the panel code, this case is controled through well considered initial discretization of the gust sheet. 39 Figure 4.10. Gust element convection along the upper airfoil surface. 4.2.3 Convection For a gust element ending on the airfoil surface, the endpoint on the surface is convected at the surface tangential velocity instead of the local flow velocity. Because the airfoil itself is discretized into a set of discrete panels and the noflow boundary condition is enforced only at each panel midpoint, a gust element endpoint convected at the local flow velocity could convect into the airfoil surface, or off the airfoil surface into the freestream flow. Basu and Hancock [3] calculated the surface tangential velocity at the gust element endpoint by interpolating tangential velocities across adjacent airfoil panels. The interpolated surface tangential velocity value was then multiplied by the local time step to find the distance the element endpoint should convect along the airfoil surface. This method provides a good first aproximation for coarse airfoil discretizations, but fails for finely discretized airfoils in locations where a large velocity gradient exists between adjacent panels, such as at the airfoil leading edge. To acount for large tangential velocity gradients, an alternate method of convecting a gust element along the airfoil surface has been developed. This alternate method estemates the 40 amount time nessisary to convect the gust endpoint along a surface panel based on the length of the surface panel and the surface tangential velocity at the panel midpoint. The estimated time to convect the gust element endpoint to the end of the surface panel is compared to the amount of time remaining in the computational timestep. Based on whether the estimated time is greater than the remaning time step, a decision is made to convect the endpoint a fractional distance along the surface panel, based on the surface tangential velocity and the remaining time step, or to convect the endpoint to the end of the current surface panel, and repeat the time estimation on the next surface panel. For example, to convect the gust element endpoint initially located at some location along Panel a, as depicted in Figure 4.10, the distance between the gust endpoint and the downstream node of Panel a is used with the surface tangential velocity at the midpoint of Panel a to estimate the amount of time necessary to convect the gust endpoint to the downstream node of Panel a. If the estimated time to convect the gust endpoint to the end of Panel a is less than the local time step, Δt, or for convenience, the time remaining, tr, then the gust endpoint is relocated to the downstream airfoil node shared by Panels a and b, and the estimated time is subtracted from tr. In this manner, using the lengths of Panels b, c, and d, along with their respective tangential velocities, the time necessary to convect the gust endpoint across Panels b, c, and d is estimated to be greater than tr, but the time necessary to convect the gust endpoint across only Panels b and c is less than tr. Thus, the gust endpoint is relocated to the shared airfoil node between Panels c and d, and the estimated time to convect the gust endpoint across Panels b and c is subtracted from tr. Because the time necessary to convect across Panel d is greater than tr, the gust endpoint is relocated a fractional distance along Panel d, as determined using the surface tangential velocity at the midpoint of Panel d and tr. 41 Figure 4.11. Gust element convection along the upper airfoil surface. 4.2.4 Gust Influence on the Airfoil Since the gust sheet is composed of singularity elements, each with an influence proportional to 1/r, the influence of a gust element ending on the airfoil depends on the proximity of the element endpoint to an airfoil collocation point. If a gust element ends on a collocation point, r approaches zero and the influence of that gust element becomes infinite. This skews the flowfield solution in a nonphysical manner. Basu and Hancock [3] prevented this possibility by replacing the each gust element ending on the airfoil with a pair of “imaginary” elements, illustrated in Figure 4.11. The two imaginary gust elements share the freestream endpoint with the original gust element, but instead of terminating at some location along an airfoil panel, airfoil panel a, with the original gust element, the imaginary element pair terminate at corresponding endpoints of airfoil panel a. The imaginary elements share the bound circulation of the original gust element in a manner dependent on the location of the original element endpoint on panel a. As such, the influence of the gust continues to propagate across the airfoil surface but the possibility of discontinuities arrising due to a gust element coexisting with an airfoil colocation point is eliminated. 42 Figure 4.12. Pitching and Plunging Airfoil 4.3 Free Response The inclusion of an airfoil structural model enables the unsteady panel code to model timedependent airfoil response due to arbitrary and selfinduced aerodynamic forcing. 4.3.1 Model The airfoil structural model is a two degree of freedom (TDOF) springmass system allowing coupled airfoil motion in rotation and translation, or pitch and plunge. Figure 4.12 shows a basic schematic detailing parameters important to the model. The equations governing twodimensional body motion in terms of sectional characteristic and generalized external forces are m¨h + Sα ¨α + mω2h h =Qh (4.1a) Sα¨h + Iα ¨α + Iαω2α α =Qα (4.1b) For a thin airfoil, the generalized external forces correspond to aerodynamic lift and moment about the elastic axis. 43 Qh = − L (4.2a) Qα =My (4.2b) For compatibility with the developed unsteady panel code, which calculates nondimensional forces and moments through integration of instantanious surfacepressure coefficents, Eq. (4.1) is nondimensionalized with respect to chord, freestream velocity, time, and mass. The resulting nondimensional equations of motion are then rewritten in terms of the nondimensional sectional characteristics, such as density ratio, μ, radius of gyration, rα, static imbalance, xα, reduced bending frequency, kh, reduced pitching frequency, kα, normalized plunge, ˆh , normalized pitch, ˆα, and nondimensional time, ˆt. μˆh + xαμ 2 ˆα + μk2h ˆh = 2 π Cl (4.3a) xαμ 2 ˆh + r2α μ 4 ˆα + r2α μk2α 4 ˆα = 2 π Cmy (4.3b) Expressing Eq. (4.3) in matrix notation, M X + K X = F (4.4) where 44 M = ⎡ ⎢⎣ μ xαμ 2 xαμ 2 r2α μ 4 ⎤ ⎥⎦ (4.5a) K = ⎡ ⎢⎣ μk2h 0 0 r2α μk2α 4 ⎤ ⎥⎦ (4.5b) X = ⎧⎪⎨ ⎪⎩ ˆh ˆα ⎫⎪⎬ ⎪⎭ (4.5c) F = 2 π ⎧⎪⎨ ⎪⎩ Cl Cmy ⎫⎪⎬ ⎪⎭ (4.5d) The second derivative in Eq. (4.4) can be isolated on the LHS, X = M −1 F − M −1 K X (4.6) allowing the equations of motion to be writen as a set of coupled first order ordinary differential equations. X = Y (4.7a) Y = M −1 F − M −1 K X (4.7b) Airfoil orientation and position at time tk+1 is determined by solving Eq. (4.7) with a fourthorder RungeKutta method using nondimensional aerodynamic forces computed at time tk. 4.3.2 Solution The aerodynamic solution and TDOF structural model are coupled directly in the developed unsteady panel method to calculate free and forced response of an arbitrary thin airfoil. The nondimensional forces and moments calulated at each time step are used as inputs to the structural model, predicting airfoil orientation and position at the next timestep. The 45 new airfoil position and orientation are used to calculate new nondimensional aerodynamic forces, and the process is repeated. 4.4 Forced Response By coupling the gust model described in Section 4.2 and the structural model described in Section 4.3, airfoil responce to arbitrary gust induced forcing can be modeled. As will be shown in Chapter 5, the influence of multiple gust sheets can be superimposed to model periodic freestream disturbance having arbitrary shapes, frequencies, and amplitudes. Thus, airfoil responce due to external forcing can be systematically studied by varying the characteristics of the freestream gust. 46 CHAPTER 5 CODE VERIFICATION To verify the accuracy and applicability of the developed panel code, a set of test cases were examined. These test cases compare unsteady panel code simulations with fundamental problems in unsteady aerodynamics having known analytical or computational solutions. In this manner, the accuracy and applicability of the panel code is established prior to its extension to problems of interest not having known solutions. 5.1 Wagner The Wagner problem, one of the fundamental problems in unsteady aerodynamics, explores the lift response of a flat plate to a flowfield which is instantaneously accelerated from one equilibrium state to another. The problem demonstrates the effect of body wake development on lift and moment during transition between equilibrium states. 5.1.1 Description Consider a stationary flat plate, or airfoil of infinitesimal thickness, at some angular orientation relative to a freestream flow, α0, illustrated in Figure 5.1. At time t < 0, the magnitude of the freestream relative to the flat plate is zero, q∞ = 0. Since the noflow boundary condition is implicitly satisfied, the body produces zero lift, and perhaps more importantly, carries zero bound circulation. At time t = 0, the freestream instantaneously accelerates to a finite nonzero velocity, q∞ = c. By applying the unsteady Kutta condition and noflow boundary condition discussed in Section 3.3, a lifting solution can be found for the flat plate. 47 Figure 5.1. Flat plate at time t = 0 It should be recalled from Section 3.3 that the body wake for an inviscid solution represents shed bound vorticity from the body which is necessary to satisfy Kelvin’s theorem. As such, the shed vorticity magnitude in the wake at time t = 0 equals the magnitude of the bound circulation change about the body, but in the opposite direction. The shed circulation caused by the flowfield transition between equilibrium states is often called a “starting vortex” because the magnitude of this vortex is significantly greater then the rest of the wake. Shed vorticity in the wake produces an aerodynamic downwash on the body, influencing the noflow boundary condition. Wake influence on lift is normally of a small magnitude relative to the freestream and the relative body motion, but in the case of a starting vortex where the shed circulation magnitude is on the same order as the bound circulation about the body, the wakeinduced downwash suppresses lift generation on the body. As such, the starting vortex significantly influences lift development on the body until the starting vortex propagates into the far field. 48 5.1.2 Solution By modeling the induced body wake as a continuous vortex sheet of varying strength, originating at the body trailing edge and oriented parallel to the freestream flow, Wagner [4] developed a timeaccurate solution for lift on an instantaneously accelerated flat plate. L = 2πbρq2α0φ (s) (5.1) This solution depends on a modified Bessel function known as the Wagner function, φ (s). φ (s) = 1 2πi ∞ −∞ C (k) k eiks dk (5.2) An approximate representation [4] of the Wagner function has been computed as, φ (s) ≈ 1 − 0.165e −0.0455s − 0.335e −0.3s (5.3) the solution of which is shown in Figure 5.2, along with the solution to the approximate Kussner function described in Section 5.4. 5.1.3 Comparison To assist verification of the developed panel code, lift solutions for thin symmetric airfoils computed using the panel code are compared the Wagner lift solution for a flat plate. Panel code solutions were obtained for instantaneously accelerated NACA 0006, 00010, and 0014 airfoils oriented at α0 = 1, 2, and 4 deg relative to a uniform freestream in the x1 direction. Solutions were computed using nondimensionalized time steps of 0.005, 0.075, and 0.010, corresponding to 4000, 3000, and 2000 computational iterations, respectively. Calculated lift coefficients for each simulation were normalized by corresponding steadystate lift values, allowing a comparison to the Wagner function, Eq. (5.3). Note that differing fundamental assumptions between the panel code and the Wagner solution affect direct comparison of the results. For example, the Wagner solution assumes 49 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 1.2 s φ(s) Approximate Wagner Approximate Kussner Figure 5.2. Solutions for the approximate Wagner function, Eq. (5.3), and the approximate Kussner function, Eq. (5.19) the body wake is a continous vortex sheet convecting at the mean freestream velocity, while the panel code discretizes the wake into a set of discrete vortices convecting at the local velocity. Also, the Wagner solution models a flat plate with negligible thickness, while the panel code models a thin symmetric airfoil. Figure 5.3 compares panel code solutions for airfoils of different thicknesses to the Wagner function. The panel code solutions in Figure 5.3 are computed for NACA 0006, 0010, and 0014 airfoils at α0 = 1 deg using a normalized time step of 0.005. As airfoil thickness decreases, the panel solutions approach the Wagner function. Figure 5.4 compares panel code solutions for a single airfoil at several orientation angles to the Wagner function. Panel code solutions in Figure 5.4 are computed for a NACA 0010 airfoil at α0 = 1, 2, and 4 deg using a normalized time step of 0.010. As Figure 5.4 shows, airfoil orientation does not have a discernable effect on normalized lift. Figure 5.5 compares panel code simulations for a single airfoil thickness and orientation but at varying normalized time steps. Panel code solutions in Figure 5.5 are computed for a 50 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t NACA 0006 NACA 0010 NACA 0014 Wagner Figure 5.3. Normalized lift for the NACA 0006, 0010, and 0014 airfoils at α0 = 1 deg using a normalized time step of 0.005 compared to Eq. (5.3) 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t α = 1 α = 2 α = 4 Wagner Figure 5.4. Normalized lift on a NACA 0010 at α0 = 1, 2, and 4 deg using a normalized time step of 0.010 compared to Eq. (5.3) 51 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t dt = 0.005 dt = 0.075 dt = 0.010 Wagner Figure 5.5. Normalized lift on a NACA 0010 at α0 = 2 deg computed using nondimensionalized time steps of 0.005, 0.075, and 0.010 compared to Eq. (5.3) NACA 0010 airfoil at α0 = 2 deg using nondimensionalized time steps of 0.005, 0.075, and 0.010. The panel code solutions show no significant dependence on the selected normalized time steps. Since neither time step nor orientation significantly affects the panel code solutions, differences between the panel code and Wagner solutions can be attributed primarily to airfoil thickness effects. Despite their differences, however, good overall agreement exists between the two lift solutions. 52 Figure 5.6. Notation used to describe the Theodorsen pitching and plunging flat plate 5.2 Theodorsen The Theodorsen problem, or the problem of a periodically pitching and plunging airfoil in an otherwise steady uniform freestream flow, demonstrates the effect of body motion and timedependent wake on unsteady lift and moments. 5.2.1 Description For an airfoil translating and rotating relative to an otherwise uniform freestream flow, induced flow perturbations near the airfoil surface due to its relative motion can be significant. For an airfoil undergoing periodic translational and rotational relative motion, the influence of the induced surface flow perturbation is a function of motion frequency and amplitude. In addition to motion induced flow perturbations, wake circulation will induce velocity perturbations which also influence the unsteady airfoil lift and moment as a function of motion frequency and amplitude. Figure 5.6 illustrates common notation used to describe the Theodorsen problem. 53 5.2.2 Solution Using conformal mapping techniques and a wake model similar to that employed in the Wagner solution, Theodorsen developed an analytic solution for lift and moment on a flat plate undergoing periodic pitching and plunging. The solution relates induced lift and moment on a flat plate to reduced frequency of the relative motion. L = LNC + LC C (k) (5.4) M = MNC +MC C (k) (5.5) Note that the Theodorsen lift and moment solutions separate noncirculatory terms, the irrotational component of lift and moment due to body motion, LNC = πρb2 ! ¨h + Uα˙ − baα¨ " (5.6) MNC = πρb2 ba¨h − Ub 1 2 − a α˙ − b2 1 8 − a2 ¨α (5.7) from circulatory terms, the rotational component of lift and moment necessitated by the Kutta condition. LC = 2πρUb ˙h + Uα + b 1 2 − a α˙ (5.8) MC = 2πρUb2 a + 1 2 ˙h + Uα + b 1 2 − a α˙ (5.9) The distinction between noncirculatory and circulatory terms is of importance because cirulatory terms depend on motion reduced frequency, as related through the Theodorsen function. C (k) = H2 1 (k) H2 1 (k) + iH2 0 (k) (5.10) 54 Combining Eqs. (5.6) through (5.8) with Eqs. (5.4) and (5.5) produces timedependent Theodorsen lift and moment equations for a flat plate undergoing periodic pitching and plunging relative to the freestream flow. L = πρb2 ! ¨h + Uα˙ − baα¨ " + 2πρUb ˙h + Uα + b 1 2 − a α˙ C (k) (5.11) My =πρb2 ba¨h − Ub 1 2 − a α˙ − b2 1 8 − a2 ¨α + 2πρUb2 a + 1 2 ˙h + Uα + b 1 2 − a α˙ C (k) (5.12) 5.2.3 Comparison To further assist verification of the unsteady panel code, namely the effects of relative body motion, computed solutions for thin symetric airfoils undergoing periodic pitching, periodic plunging, and periodic pitching and plunging are compared to the corresponding Theodorsen solution for a flat plate. As with comparisons to the Wagner function, the effects of thickness and wake model limit direct comparison between the panel code and analytic solutions. 5.2.3.1 Pure Pitching Panel code solutions for NACA 0006, 0010, and 0014 airfoils pitching relative to the freestream flow were computed for reduced frequencies of k = 0.25 and 0.75 and amplitudes of ¯α = 1, 2, and 4 deg about the airfoil quarterchord location. Timedependent lift and moment for a flat plate undergoing similar motion were also computed using Eqs. (5.11) and (5.12). Figures 5.7 through 5.9 demonstrate the effect of airfoil thickness on panel code lift and moment solutions, as compared to the Theodorsen solution. Panel code solutions in Figures 5.7 through 5.9 were computed for NACA 0006, 0010, and 0014 airfoils pitching at a reduced frequency of k = 0.25 and an amplitude of ¯α = 2deg. 55 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.7. Cl vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg For pure pitching at small reduced frequencies, airfoil thickness exhibits a small influence on the phase between the panel code and Theodorsen lift solutions as well as the amplitude ratio of the two solutions. However, both amplitude and phase of the panel code solution approach the Theodorsen solution as airfoil thickness decreases. Figures 5.10 through 5.12 shows the relative agreement between the panel code lift and moment to the Theodorsen solution, for a range of pitching amplitudes. Panel code solutions for Figures 5.10 through 5.12 were computed for a NACA 0010 airfoil pitching at a reduced frequency of k = 0.25, and pitching amplitudes of ¯α = 1, 2, and 4 deg. For pure pitching at a constant reduced frequency, pitching amplitude does not appear to exhibit a significant influence on either the phase between the panel code and Theodorsen lift solutions or the lift ratio between the two solutions. The lift ratio, computed as the maximum panel code lift coefficent divided by the maximum Theodorsen lift coefficient, remains constant around 1.08 for the pitching amplitudes computed. 56 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmle t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.8. Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmea t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.9. Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg 57 0 2 4 6 8 10 12 14 16 18 20 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Cl t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.10. Cl vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 Cmle t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.11. Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg 58 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.12. Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg Figures 5.13 through 5.15 show the relative agreement of the panel code lift and moment solutions to the Theodorsen solution, for a range of reduced frequencies. Panel code solutions in Figures 5.13 through 5.15 are a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with a pitching amplitude of ¯α = 2 deg. For pure pitching at a constant amplitude, reduced frequency does not appear to exhibit an influence on the phase between the panel code and Theodorsen lift solutions, but does appear to influence the amplitude ratio between the two solutions. It appears that the pase between the panel code and Theodorsen lift solutions remains constant as reduced frequency varies, however, the amplitude ratio between the panel code and Theodorsen lift solution decreases as reduced frequency increases. 59 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.13. Cl vs. Time for a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cmle t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.14. Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 60 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.15. Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 61 0 2 4 6 8 10 12 14 16 18 20 −0.025 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 Cl t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.16. Cl vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 5.2.3.2 Pure Plunging The unsteady panel code was used to generate solutions for NACA 0006, 0010, and 0014 airfoils plunging at reduced frequencies of k = 0.25 and 0.75 with plunging amplitudes of ¯h = 0.010, 0.025, and 0.050 relative to the mean freestream flow. Timedependent lift and moment for a flat plate undergoing similar motion were also computed using Eqs. (5.11) and (5.12). The halfchord was used to calculate moments about the elastic axis. Figures 5.16 through 5.18 demonstrate the effect of airfoil thickness on panel code lift and moment solutions, as compared to the Theodorsen solution. Panel code solutions in Figures 5.16 through 5.18 were computed for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and amplitude of ¯h = 0.025. For pure plunging at small reduced frequencies, as in the case of pure pitching, airfoil thickness exhibits a small influence on the phase between the panel code and Theodorsen lift solutions as well as the amplitude ratio of the two solutions. However, both amplitude and phase of the panel code solution approach the Theodorsen solution as airfoil thickness decreases. 62 0 2 4 6 8 10 12 14 16 18 20 −0.01 −0.005 0 0.005 0.01 Cmle t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.17. Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 0 2 4 6 8 10 12 14 16 18 20 −6 −4 −2 0 2 4 6 x 10−3 Cmea t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.18. Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 63 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cl t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.19. Cl vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 Figures 5.19 through 5.21 show the relative agreement between the panel code lift and moment solution to the Theodorsen solution for different plunging amplitudes. Panel code solutions in Figures 5.19 through 5.21 were computed for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050. As is the case for pure pitching, for pure plunging at a constant reduced frequency, pitching amplitude does not appear to exhibit a significant influence on either the phase between the panel code and Theodorsen lift solutions or the lift ratio between the two solutions. The amplitude ratio remains constant around 0.99 for the plunging amplitudes computed. Figures 5.22 though 5.24 show relative agreement between the panel code lift and moment solution to the Theodorsens solution at different reduced frequencies. Panel code solutions in Figures 5.22 though 5.24 were computed for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and amplitude of ¯h = 0.025. As is the case for pure pitching, for pure plunging at a constant amplitude, reduced frequency does not appear to exhibit an influence on the phase between the panel code and Theodorsen lift solutions, but does appear 64 0 2 4 6 8 10 12 14 16 18 20 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 Cmle t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.20. Cmle vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 0 2 4 6 8 10 12 14 16 18 20 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04 0.05 Cmea t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.21. Cmea vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 65 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.22. Cl vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 to influence the amplitude ratio between the two solutions. It appears that the pase between the panel code and Theodorsen lift solutions remains constant as reduced frequency varies, however, the amplitude ratio between the panel code and Theodorsen lift solution decreases as reduced frequency increases. 66 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cmle t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.23. Cmle vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.24. Cmea vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 67 0 2 4 6 8 10 12 14 16 18 20 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Cl t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.25. Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 5.2.3.3 Combined Pitching and Plunging Because solutions for an airfoil undergoing combined pitching and plunging motion represent a superposition of solutions for pure pitching and pure plunging, which have already been examined, this section will use a subset of the previously examined test cases to demonstrate that the panel code properly models combined pitching and plunging. It is expected that the same observations on the effects of amplitude and reduced frequency for the pure pitching and pure plunging cases will hold for the combined pitching and plunging. Figures 5.25 through 5.27 demonstrate the relative agreement between the panel code lift and moment solutions and the Theodorsen solution. Panel code solutions in Figures 5.25 through 5.27 are for a NACA 0010 airfoil pitching and plunging about the quarterchord at a reduced frequency of k = 0.25 with amplitudes of ¯α = 1, 2, and 4 deg, and ¯h = 0.025. Figures 5.28 through 5.30 demonstrate the relative agreement between the panel code lift and moment solutions and the Theodrsen solution. Panel code solutions in Figures 5.28 through 5.30 are for a NACA 0010 airfoil pitching and plunging about the quarterchord at 68 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 Cmle t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.26. Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.27. Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 69 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.28. Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. a reduced frequency of k = 0.25 with a amplitudes of ¯α = 1 deg and ¯h = 0.010, 0.025, and 0.050. 5.2.3.4 Discussion As observed in the pure pitching and pure plunging examples, the panel code solution showed small variations in both phase and amplitude as compared to the Theodorsen solution. It was also shown that these small variations were dependent only on airfoil thickness and reduced frequency. Since the variations do not do not show a dependence on motion amplitude, the variation between the two solutions may be attributed to differences inherent in the wake models. As described earlier, the Theodorsen solution models the shed bound vorticity as a continuous vortex sheet of variable strength, released from the undisturbed airfoil trailing edge. The vortex sheet then convects with the timeaveraged freestream flow, essentially confining the vortex sheet to the x1 axis. The panel code models the shed circulation as a set of discrete vortex elements, released from the airfoil trailing edge location at each time step. The discrete 70 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmle t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.29. Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmea t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.30. Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. 71 vortex elements are allowed to convect with the instantaneous local flowfield, and thus the wake is allowed to deform in time. By convecting the wake at the local velocity and allowing the wake to deform, the panel code wake model induces a different influence on the airfoil than the Theodorsen wake model, which would be dependent on reduced frequency. 72 Figure 5.31. Stationary plate of infinitesimal thickness with periodic transverse gust 5.3 Sears Periodic Gust The Sears periodic gust problem examines timedependent lift and moment generated on a stationary airfoil under the influence of a timeaveraged uniform freestream flow with sinusoidal transverse velocity perturbations, or transverse gusts. 5.3.1 Description Consider a stationary airfoil immersed in a freestream flow where the timeaveraged flow is aligned with the airfoil chord. If the airfoil is symmetric, it will not generate lift. However, if a transverse velocity perturbation is introduced into the freestream, the velocity perturbation will induce a local timedependent angle of attack change on the airfoil that will induce unsteady lift. The Sears gust problem investigates the influence of periodic transverse velocity perturbations on unsteady lift and moment generated on a flat plate. 73 5.3.2 Solution By assuming the transverse gust is not influenced by the presence of the airfoil, i.e. the “frozen gust” approximation, the downwash induced by the gust on the airfoil as a function of time can be writen as w = ¯weiω(t− x U ) = ¯weik(s−x∗) (5.13) Using a wake model similar to that of Wagner and Theodorsen, the relation for wake induced downwash on the airfoil as a function of reduced gust frequency is given by the Sears function. S (k) = 2 πk [H2 0 (k) − iH2 1 (k)] (5.14) Modifying the noflow boundary condition to include both gust and wake induced downwash, Sears lift and moment solutions for a flat plate under the influence of a sinusoidal transverse gust become L = 2πρU ¯ wbeiωt S (k) (5.15) My = L 1 2 + a b (5.16) 5.3.3 Comparison The panel code can implement two separate methods to model a periodic transverse gust. The first method, refered to later as the modified panel code, accounts for the influence of the periodic gust directly by modifying the implementation of the noflow boundary condition on the airfoil surface. The modified boundary condition includes the influence of the velocity pertubation by replacing the constant freestream velocity term with a time and position dependent function. This time and position dependent function must then be included in every other calculation which depends on the freestream velocity, such as the computation of unsteady surface pressures and the convection routines used to convect wake elements with the local flowfield. The application of this method to other problems, such as forced responce, is limited because gust influence on the airfoil is directly modeled as a function of 74 time and location in the flowfield, and as such does not allow for gust deformation due to body or wake influences. The second method employs the gust model descibed in Section 4.2. To use this discrete gust model, the continuous periodic gust is discretized into a set of gust sheets which propagate across the airfoil at the local flow velocity. This method does not require modifications to the original noflow boundary condition since the gust sheet is composed of constant strength vortex elements whose influence was included in the original noflow boundary condition. Since the gust sheets convect with the local flowfield, this method allows for gust deformation due to body and wake influences. The comparison of this method to the Sears solution is only limited by the discretization of the continuous periodic gust into a corresponding gust sheet representation. As such, care must be taken in choosing the method of discretization, since different representations of the same continuous gust will result in different lift and moment solutions. 5.3.3.1 Modified NoFlow Boundary Condition Figures 5.32 through 5.37 demonstrate the effect of airfoil thickness on panel code lift and moment solutions as compared to the Sears lift and moment solutions for reduced frequencies of k = 0.25, 1.0, and 4.0. Panel code solutions in Figures 5.32 through 5.37 were computed for NACA 0006, 0010, 0012, and 00014 airfoils under the influence of a continuous sinusoidal gust having an amplitude of ¯ w = 0.01 using the modified noflow boundary condition. Airfoil thickness exhibits a small influence on the phase between the modified panel code and Sears lift solutions as well as the amplitude ratio of the two solutions. The difference in amplitude of the lift solution computed by the panel code is attributed to the effects of airfoil thickness, because the solutions approach, but do not reach the Sears solution as airfoil thickness decreases. It is interesting to note that the amplitude ratio, computed as the maximum panel code lift coefficent divided by the maximum Sears lift coefficent, remains constant at roughly 0.7 for a NACA 0010 airfoil regardless of reduced frequency. The 75 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cl t Sears k = 0.25 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.32. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 influence of reduced frequency on the phase of the solutions is more dramatic. The phase of the modified panel code solution lags the Sears solution at small reduced frequencies, and shifts such that it leads at higher reduced frequencies. Because the amplitude ratio does not appear to be influenced by reduced frequency, it is assumed that differences in wake models, as discussed in Section 5.2.3.4, are responsible for the phase shift with reduced frequency. 5.3.3.2 Vortex Sheet Gust Model Figures 5.38 and 5.39 compares lift and moment coefficents calculated by the panel code utilizing the freestream gust model to a corresponding Sears solution. Panel code solutios in Figures 5.38 and 5.39 were computed for a NACA 0010 airfoil oriented at α0 = 0.0 deg to the timeaveraged freestream. The gust, having a reduced frequency of k = 1.0 and amplitude of ¯ w = 0.01, was modeled using a set of six gust sheets per gust period for five and a half periods upstream of the airfoil. The strength of each gust sheet is based on the velocity pertubation at the gust sheet’s initial x1 location in the flowfield. Figure 5.40 shows 76 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmle t Sears k = 0.25 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.33. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cl t Sears k = 1.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.34. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 77 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 Cmle t Sears k = 1.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.35. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 0 1 2 3 4 5 6 7 8 9 10 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 Cl t Sears k = 4.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.36. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 78 0 1 2 3 4 5 6 7 8 9 10 −0.01 −0.005 0 0.005 0.01 0.015 Cmle t Sears k = 4.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.37. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 circulation strength per unit length about each gust sheet relative to the initial x1 location of the sheet. One drawback to the use of the discrete gust model is that the gust sheets do not produce a sinusoidal velocity pertubation. The induced pertubation due to the gust sheets closely resembles a set of sharp edge gusts, as shown in Figure 5.41. Figure 5.41 combines a visualization of the location of the airfoil, wake, and gust sheets in the top panel, the instantaneous pressure coefficients along the airfoil in the lower
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Development and Validation of an Unsteady Panel Code to Model Airfoil Aeromechanical Response 
Date  20040701 
Author  McClung, Aaron 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  An unsteady panel code has been developed as a computational tool to investigate the influence of aerodynamic damping on airfoil aeromechanical response. The inclusion of a freestream gust model enables the panel code to simulate freestream gust perturbations of arbitrary shape and magnitude. The inclusion of a two degree of freedom structural model also enables the panel code to model structural response due to both selfinduced and gustinduced aerodynamic forcing. Together, the gust and structural models allow aeromechanical response to be modeled. Panel code solutions compared favorably to classic problems in unsteady aerodynamics, however, panel code solutions utilizing the freestream gust model show a significant dependence on freestream gust model parameters and gust discretization. Aeromechanical response was demonstrated using the panel code, however, further refinement of the gust model is necessary prior to application of the panel code to investigations of interest. 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  DEVELOPMENT AND VALIDATION OF AN UNSTEADY PANEL CODE TO MODEL AIRFOIL AEROMECHANICAL RESPONSE By AARON M. MCCLUNG Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE July 2004 DEVELOPMENT AND VALIDATION OF AN UNSTEADY PANEL CODE TO MODEL AIRFOIL AEROMECHANICAL RESPONSE By AARON M. MCCLUNG Approved as to style and content by: Eric A. Falk, Chair Andrew S. Arena, Member Gary E. Young, Member Dean, Graduate College ii TABLE OF CONTENTS Page LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii NOTATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv CHAPTER 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2. FUNDAMENTALS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2 Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.3 NavierStokes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.4 Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Potential Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Velocity Potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.2 Superposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Theorems And Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Bernoulii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2 Coefficient of Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Angular Velocity, Vorticity, and Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.1 Motion of a Fluid Element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4.2 Angular Velocity and Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 iii 2.4.3 Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.4 Kelvin’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3. PANEL CODES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1 Nonlifting Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1 Discretization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Lifting Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1 Kutta Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.2 Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.3 TimeDependent Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.1 Frame of Reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.2 Wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.3.3 Unsteady Kutta Condtion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.3.1 BasuHancock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.3.2 Ardonceau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.4 Method of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4. CODE DESCRIPTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 Frame of Reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Gust Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2.1 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2.2 AirfoilGust Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2.2.1 Determining Gust Element Condition . . . . . . . . . . . . . . . . . . . . . . 34 4.2.2.2 Case One . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.2.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2.2.4 Case Two . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.3 Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2.4 Gust Influence on the Airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.3 Free Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 iv 4.4 Forced Response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5. CODE VERIFICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1 Wagner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Theodorsen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.2.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2.3.1 Pure Pitching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2.3.2 Pure Plunging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2.3.3 Combined Pitching and Plunging. . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2.3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Sears Periodic Gust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3.3.1 Modified NoFlow Boundary Condition . . . . . . . . . . . . . . . . . . . . 75 5.3.3.2 Vortex Sheet Gust Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Kussner’s Sharp Edge Gust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4.2 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4.3.1 Transient Panel Code Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.4.3.2 Single and Double Gust Sheets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.5 Free Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.5.1 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.5.2 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 v 6. FORCED RESPONSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7. SUMMARY AND CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2 Extension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 7.3 Discussion and Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.4 Contributions of Present Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 APPENDICES A. UVPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 A.1 Revision 120 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 B. COMMON FILES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1 Common Variables Declarations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.1 lengths.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.2 airfoil.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.1.3 calc.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.4 const.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.5 debug.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.6 file.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 B.1.7 forces.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 B.1.8 freeresp.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 B.1.9 freevort.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.10 gau.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.11 graph.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 B.1.12 iterative.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.13 motion.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.14 param.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.15 phi.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 B.1.16 relax.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.17 strengths.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.18 velocities.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.19 wake.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 B.1.20 wakepannel.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 vi B.1.21 graph cons.inc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 C. INPUT FILES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 C.1 Configuration File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 C.2 Airfoil Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 C.3 Motion History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 C.4 Free Stream Vortices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 D. GRAPHICS ROUTINES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.1 Plotting Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.1.1 Compare Data r10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 D.2 Animation Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 D.2.1 Animate r21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 vii LIST OF FIGURES Figure Page 3.1 Airfoil modeled with a continuous source sheet. . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Airfoil discretized into constant strength source elements. . . . . . . . . . . . . . . . . . . 15 3.3 Constant Strength Panel Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.1 Influence of a Vortex sheet located in the Freestream flow compared to the Freestream influence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Deformation of a vortex sheet approaching the airfoil leading edge. . . . . . . . . . . 32 4.3 Case One: Gust element straddling the airfoil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Case Two: Gust element endpoint convected into the airfoil. . . . . . . . . . . . . . . . . 34 4.5 Airfoil leading edge vs. the airfoil forwardmost node. . . . . . . . . . . . . . . . . . . . . . 36 4.6 Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+1. . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.7 Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+2. . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.8 Gust element split about the upstream stagnation point. . . . . . . . . . . . . . . . . . . . 38 4.9 Interpolation to determine the time of GustAirfoil impact. . . . . . . . . . . . . . . . . . 39 4.10 Gust element convection along the upper airfoil surface. . . . . . . . . . . . . . . . . . . . . 40 4.11 Gust element convection along the upper airfoil surface. . . . . . . . . . . . . . . . . . . . . 42 4.12 Pitching and Plunging Airfoil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.1 Flat plate at time t = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 viii 5.2 Solutions for the approximate Wagner function, Eq. (5.3), and the approximate Kussner function, Eq. (5.19) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.3 Normalized lift for the NACA 0006, 0010, and 0014 airfoils at α0 = 1deg using a normalized time step of 0.005 compared to Eq. (5.3) . . . . . . . . . . . . . 51 5.4 Normalized lift on a NACA 0010 at α0 = 1, 2, and 4 deg using a normalized time step of 0.010 compared to Eq. (5.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.5 Normalized lift on a NACA 0010 at α0 = 2 deg computed using nondimensionalized time steps of 0.005, 0.075, and 0.010 compared to Eq. (5.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.6 Notation used to describe the Theodorsen pitching and plunging flat plate 53 5.7 Cl vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.8 Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.9 Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.10 Cl vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 58 5.11 Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 58 5.12 Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg . . . . . . 59 5.13 Cl vs. Time for a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.14 Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 60 ix 5.15 Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 61 5.16 Cl vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 62 5.17 Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 63 5.18 Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . 63 5.19 Cl vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . . . . . 64 5.20 Cmle vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . 65 5.21 Cmea vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 . . . . . . . . . . . . . . . . . 65 5.22 Cl vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.23 Cmle vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . 67 5.24 Cmea vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 . . . . . . . . . . . . . . . . . . . . . . . 67 5.25 Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.26 Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . 69 5.27 Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. . . . . . . . . . . . . . . . . . . . . . . . 69 5.28 Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . . . 70 x 5.29 Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . 71 5.30 Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. . . . . . . . . . . . . . . . . 71 5.31 Stationary plate of infinitesimal thickness with periodic transverse gust . . . . . . 73 5.32 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.33 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.34 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.35 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.36 Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.37 Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.38 Cl vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 xi 5.39 Cmle vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.40 Gust sheet circulation per unit length vs. initial x/c location for a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 6 times the reduced frequency. . . . . . . . . . . . . . . . . . 81 5.41 Visualization showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . . . . . 82 5.42 Cl vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.43 Cmle vs. Time for the Sears solution compared to the panel code solution for NACA 0010 airfoil under the influence of a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.44 Gust sheet circulation per unit length vs. initial x/c location for a periodic freestream gust with a reduced frequency of k = 1.0 and gust amplitude of ¯ w = 0.01, sampled at 4 times the reduced frequency. . . . . . . . . . . . . . . . . . 84 5.45 Visualization at t = 2.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 85 5.46 Stationary plate of infinitesimal thickness with sharp edge transverse gust . . . . . 86 5.47 Transient lift solutions normilized by the corresponding steady state lift for NACA 0008, 0010, and 0012 airfoils oriented at at α0 = 1deg relative to the timeaveraged freestream computes using a normalized time step of 0.005 compared to Eq. (5.19) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.48 Cl for a single gust sheet of strength γ = −0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gust with an amplitude of ¯ w = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 xii 5.49 Cmle for a single gust sheet of strength γ = −0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gust with an amplitude of ¯ w = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.50 Visualization at t = 2.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 92 5.51 Cl for a pair of gust sheets of strength γ = 0.02 and 0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gusts with amplitudes of ¯ w = 0.01 and 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.52 Cmle for a pair of gust sheets of strength γ = 0.02 and 0.02 propagating across a NACA 0010 airfoil oriented at α0 = 0.0 to the timeaveraged freestream computed using a time step of 0.010 compared to the Kussner sharp edge gusts with amplitudes of ¯ w = 0.01 and 0.01. . . . . . . . . . . . . . . . . 93 5.53 Visualization at t = 4.0 showing the location of the airfoil, wake, gust sheets, and selected x2 velocities in the top panel, instantaneous Cp vs. x/c in the lower left panel, and Cl vs. t in the lower right panel. . . . . . . . . . . . . . . . 94 5.54 Pitch history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.55 Plunge history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.56 Cl history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.57 Cmle history for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.58 Plunge vs. Pitch for the panel code free response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.59 Determining modal damping for pitch. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.60 Determining modal damping for plunge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.61 Normalized modal damping vs. freestream velocity for the panel code solution compared to the pk method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 xiii 5.62 Normalized modal frequency vs. freestream velocity for the panel code solution compared to the pk method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.1 Pitch history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Plunge history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Cl history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.4 Cmle history for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.5 h/c vs. α for the panel code forced response simulation above, at, and below the predicted flutter boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 xiv NOTATION Φ Velocity Potential ρ Density Γ Total circulation strength γ Circulation strength per unit length λ Source strength per unit length Λ Total source strength n Normal Vector V Velocity Vector q Velocity Vector rij Radius from poing j to point i mc.v. Mass of fluid inside the control volume mout Mass flux out of a control volume min Mass flux into a control volume p Momentum a Acceleration F Force ds Differential along a surface f Body force τij Fluid shear stress p Pressure μ Viscocity coefficient δij Kronecker delta function ω Angular velocity of a fluid element ζ Vorticity of a fluid element L Lift Cl Coefficient of Lift Cd Coefficient of Drag Cm Coefficient of Moment c Chord, length of the airfoil section b Semichord = c 2 s = Ut/b, Semichord location fi Body Force τij Fluid Shear Stress Qh = −L, Generalized force along the +zaxis Qα = My, Generalized moment about the elastic axis m Mass xv h Vertical translation of the airfoil, positive for deflection along the zaxis α Angle between the airfoil centerline and the mean freestream flow Iα Mass moment of inertial per unit span about axis x = ba Sα = mbxα, Static mass imbalance per unit span about axis x = ba ωh = Kh/m, uncoupled natural frequency in bending ωα = Kα/Iα, uncoupled natural frequency in torsion Kh Bending spring stiffness Kα Torsional spring stiffness k = ωb/U, Reduced frequency xvi CHAPTER 1 INTRODUCTION Aeroelastic condsiderations affect a wide range of disciplines. With respect to turbomachinery, particularly the area of highcycle fatigue, aerodynamic forcing of internal components due to rotorstator interactions can significantly impact engine lifecycle and maintenance requirements. To better understand the influence of aerodynamic damping, on highcycle fatigue, the influence of aerodynamic damping on forced structural response must first be be examined. As a first step towards this goal, this thesis develops a computational tool through which the influence of aerodynamic damping can be isolated and systematically studied. 1.1 Goals The goal of this thesis is to develop and validate a computation tool which will enable the systematic investigation into wake induced stuctural responce. The computational tool is based loosely on a HessSmith [5] type unsteady panel code written by Ron Hugo [7, 9] which has been modified to include a freestream gust model and an airfoil structural model. By incorporating the capability to model arbitrary freestream gusts into the unsteady panel code, and coupling the panel code with a structural model, the timedomain response of a body due to an arbitrary freestream disturbance can be computed. 1.2 Organization This thesis is presented in five parts. The first part is an overview of the governing fluid dynamic equations and the derivation of velocity potential which governs the inviscid and 1 incompressible flowfield, as well as the derivation and description of related theorems and concepts which are necessary for the formulation of the numeric solution. The second part describes the formulation of two dimensional panel methods in three sections, starting with the formulation to solve the steadystate flowfield about a nonlifting body, adding the Kutta condition to solve the steadystate flowfield about a lifting body, and then accounting for timedependent effects to solve the timedependent flowfield about a liftingbody undergoing arbitrary motion. The third part of the thesis expands on the timedependent panel method formulation by adding a freestream gust model which represents time dependent freestream pertubations using discrete vortex elements, and a two degree of freedom structural model which allows the responce of an arbitrary body due to aerodynamic forcing to be determined. The fourth part compares the developed panel code against classic analytic solutions for unsteady aerodynamics, and the last part demonstrates the application of the developed panel code to a forced responce problem using a solution which couples the freestream gust model with the structural model. 2 CHAPTER 2 FUNDAMENTALS Before panel codes are discussed, this chapter defines several relations and terms used throughout the later discussion. The first section in the present chapter discusses the derivation of basic governing equations for fluid flow. The second section discusses potential flow and applies basic governing equations to the solution of potential flowfields. The last two sections relate terms and definitions used later in this thesis. 2.1 Governing Equations The fundamental equations governing fluid flow are derived here from the relationships between density, momentum, and energy, and their time rates of change inside a controlvolume. 2.1.1 Continuity The continuity equation relates the time rate of change of mass inside a controlvolume to the mass flux through the controlsurface. The integral form of the continuity equation can be derived by beginning with a statement of mass inside a controlvolume, such as mc.v. = c.v. ρ dV (2.1) Based on Eq. (2.1), the time rate of change of mass inside the controlvolume, ∂mc.v./∂t, is given by ∂mc.v. ∂t = ∂ ∂t c.v. ρ dV (2.2) 3 Mass flux through the controlsurface can also be stated as m˙ in − m˙ out = −c.s. ρqini dS (2.3) If mass is conserved, the net mass flux through the controlsurface must equal the time rate of change of the mass within the controlvolume, leading to the integral form of the continuity equation [8]. ∂ ∂t mc.v. = ∂ ∂t c.v. ρ dV = −c.s. ρqini dS = ˙ min − m˙ out (2.4) The divergence theorem states that given a vector qi, the integral of the normal component of qi relative to the controlsurface equals the integral of the gradient of qi inside the corresponding controlvolume. c.s. qini dS = c.v. ∂ ∂xi qi dV (2.5) By applying Eq. (2.5) to the integral form of the conservation equation, Eq. (2.4), the following simplification can be made c.s. ρqini dS = c.v. ∂ ∂xi (ρqi) dV (2.6) Thus, Eq. (2.4) can be reduced to ∂ ∂t c.v. ρ dV + c.v. ∂ ∂xi (ρqi) dV = 0 (2.7) or c.v. ∂ ∂t ρ + ∂ ∂xi (ρqi) dV = 0 (2.8) Since the volume integral in Eq. (2.8) must equal zero for any arbitrary controlvolume, it must also hold that ∂ ∂t ρ + ∂ ∂xi (ρqi) = 0 (2.9) producing the differential form of the continuity equation [8]. 4 2.1.2 Momentum The momentum equation relates the time rate of change of fluid momentum through a controlvolume to the forces acting on the controlvolume. Momentum is a vector quantity, pj , defined by the product of mass and the corresponding velocity vector. pj = mqj (2.10) For a controlvolume, the summation of forces acting on the volume equal the time rate of change of the controlvolume momentum. c.v. Fj = ∂ ∂t (mqj)c.v. (2.11) When Eq. (2.11) is incorporated with the continuity equation, Eq. (2.4) becomes c.v. Fj = ∂ ∂t (mqj)c.v. = ∂ ∂t c.v. ρqj dV + c.s. ρqjqini dS (2.12) Forces acting on the controlvolume may be either body forces, surface forces, or both. c.v. Fbodyj = c.v. ρfj dV (2.13) c.s. Fsurfacej = c.s. τijni dS (2.14) Thus, substituting Eqs. (2.12) –(2.14) into Eq. (2.11) gives the integral form of the momentum equation. ∂ ∂t c.v. ρqj dV + c.s. ρqjqini dS = c.v. ρfj dV + c.s. τijni dS (2.15) Applying the divergence theorem to Eqs. (2.12) and (2.14) allows simplification of Eq. (2.15) c.s. ρqjqini dS = c.v. ∂ ∂xi (ρqjqi) dV (2.16) 5 c.s. τijni dS = c.v. ∂ ∂xi τij dV (2.17) c.v. ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) − ρfj − ∂ ∂xi τij dV = 0 (2.18) Again, since the volume integral must equal zero for any arbitrary controlvolume, it holds that ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) = ρfj + ∂ ∂xi τij (2.19) producing the differential form of the momentum equation. 2.1.3 NavierStokes If the assumption is made that the fluid is Newtonian (i.e. the stress components τij are linearly related to the derivatives ∂qi/∂xj ), then the following substitution has been widely accepted τij = − p + 2 3 μ ∂qk ∂xk δij + μ ∂qi ∂xj + ∂qj ∂xi (2.20) Thus, Eq. (2.20) can be substituted into Eq. (2.19), giving conservative form of the Navier Stokes relation [8]. ∂ ∂t (ρqj) + ∂ ∂xi (ρqjqi) = ρfi − ∂ ∂xj p + 2 3 μ ∂qk ∂xk + ∂ ∂xi μ ∂qi ∂xj + ∂qj ∂xi (2.21) 2.1.4 Euler Depending on the flow regime, the NavierStokes equations can be simplified. For example, lowspeed flow about a thin airfoil outside of the boundary layer can be assumed to be incompressible, ρ = constant, and inviscid, μ = 0, if the airfoil is at a conservative angle of attack and large Reynolds Numbers. With these two assumptions, Eq. (2.21) simplifies to the Euler equation. ∂ ∂t qj + ∂ ∂xi (qjqi) = fj + 1 ρ ∂p ∂xj (2.22) 6 2.2 Potential Flow The potential flow assumption is of interest here because it describes the flow regime examined in the current investigation. 2.2.1 Velocity Potential If a flowfield can be considered incompressible, then the continuity equation, Eq. (2.9), simplifies to ∂qi/∂xi = 0. If the flowfield is also inviscid, μ = 0, then vorticity in the flowfield must remain constant with respect to time, ∂ζ/∂t = 0. Given these assumptions, a scaler potential function Φ exists that is a solution to the Laplace equation describing the flowfield ∂2 ∂x2j Φ = 0 (2.23) The potential function, Φ, is often denoted the velocity potential because the velocity field is equal to the gradient of Φ. qj = ∂ ∂xj Φ (2.24) Inversely, the potential at any point, P, in the flowfield can be calculated from any arbitrary reference point, P0, by integrating the velocity field along any path between P0 and P Φ(x1, x2, x3) = P P0 x1 dx1 + x2 dx2 + x3 dx3 = P P0 ∂Φ ∂x1 dx1 + ∂Φ ∂x2 dx2 + ∂Φ ∂x3 dx3 (2.25) Note that with the assumptions of irrotationality and incompressibility, the integrand of Eq. (2.25) is an exact differential, and as such the potential is independent of the integration path. [8] 2.2.2 Superposition Because the velocity potential describes the potential flowfield, and is the solution to the Laplace equation, it holds that [1]: 7 1. Any irrotational incompressible flow has a velocity potential and stream function (for twodimensional flow) that both satisfy Laplace’s equation. 2. Conversely, any solution of Laplace’s equation represents the velocity potential or stream function (twodimensional) for an irrotational, incompressible flow. Since the Laplace equation is a secondorder, linear, partial differential equation, it holds that the sum of two or more particular solutions is also a valid solution. Thus, a complex flowfield with a total potential Φ can be modeled as the superposition of multiple potential solutions, Φk, giving Φ = Φk (2.26) 2.2.3 Boundary Conditions Since solving the Laplace equation is a boundary value problem, applying the correct boundary conditions is essential. The two physical phenomon considered here are the noflow boundary condition at the fluidbody interface, and the farfield condition forcing bodyinduced disturbances to decay to zero strength far from the body. There are two types of boundary condition formulations, the “direct” Neumann boundary condition, and the “indirect” Dirichlet boundary condition. The Dirichlet boundary condition is not explained here because it is not employed in this investigation. See References [1] and [8] for a full explanation. The Neumann boundary condition specifies the normal velocity on the fluidbody boundary must equal zero, ∂Φ ∂n = 0 (2.27) and the potential field due to the presence of the body must be negligible in the farfield (r→∞). lim r→∞ (Φbody) = 0 (2.28) 8 2.3 Theorems And Relations 2.3.1 Bernoulii To compute pressure in a potential flow, the relation between potential and velocity, qj = ∂Φ/∂xj , and the assumption of a conservative body force with potential E, fj = −∂E/∂xj , are substituted in to the Euler equation, Eq. (2.22). ∂ ∂xj E + p ρ + qj 2 2 + ∂Φ ∂t = 0 (2.29) Thus, upon spatial integration E + p ρ + qj 2 2 + ∂Φ ∂t = C(t) (2.30) where C(t) is a spatially independent constant over the entire flowfield, but is a function of time. This is the Bernoulli equation [8]. Because the left hand side of Eq. (2.30) is constant over the entire flowfield at a given point in time, pressure and velocity can be compared at different points in the flow if the potential is known. 2.3.2 Coefficient of Pressure The pressure coefficient is a nondimensional parameter relating pressure between two different locations in a flowfield. Cp = p∞ − p 1 2ρqj 2 ∞ (2.31) Using the Bernoulli equation, Eq. (2.30), the pressure difference in Eq. (2.31) becomes p∞ − p = ρ E + qj 2 2 + ∂Φ ∂t ∞ − ρ E + qj 2 2 + ∂Φ ∂t p (2.32) If care is taken with the choice of reference point, denoted by ∞, such that it exists at a location in the farfield not influenced by any bodyinduced disturbances, then the change in potential with time can be neglected at the reference point. 9 ∂Φ∞ ∂t = 0 (2.33) If the reference point is also chosen such that the difference in the body forces is negligible, E∞ = Ep (2.34) then Eq. (2.32) reduces to p∞ − p = ρ qj 2 ∞ 2 − qj 2 p 2 − ∂Φp ∂t (2.35) Dividing the pressure difference by the free stream dynamic pressure, 1 2ρqj 2 ∞, Eq. (2.35) becomes Cp = 1− qj 2 p qj 2 ∞ − 2 qj 2 ∞ ∂Φp ∂t (2.36) 2.4 Angular Velocity, Vorticity, and Circulation 2.4.1 Motion of a Fluid Element Motion of a fluid element is comprised of translation, rotation, and deformation, where each type of motion is usually caused by different phenomena in the flowfield. Translation is caused by a uniform velocity, where all parts of the element move at the same velocity, disallowing deformation and rotation. Rotation and deformation occur when velocity gradients exist across the element, as can be the case when viscous effects are not negligible. 2.4.2 Angular Velocity and Vorticity The angular velocity of a fluid element relates to the element deformation caused by a velocity gradient. Generally these velocity gradients are caused by shear stresses. The angular velocity of a fluid element, ωi, is defined as the curl of the velocity vector, or ωi = −1 2 εijk ∂qj ∂xk (2.37) 10 Another measure of fluid angular velocity, used to simplify several equations, is vorticity, defined as twice the angular velocity. ζi = 2ωi = −εijk ∂qj ∂xk (2.38) 2.4.3 Circulation Circulation, Γ, is a measure of the vorticity in a fluid region, and equals the integral of the vorticity normal to a surface, S. Γ = S ζini dS (2.39) By substituting the definition of vorticity, Eq. (2.38), into Equation (2.39) and using Stokes Theorm [8], S −εijk ∂qj ∂xk ni dS = C qi dxi (2.40) circulation can be defined as Γ = C qi dxi (2.41) 2.4.4 Kelvin’s Theorem Kelvin’s theorem relates the time rate of change of circulation in a potential flow inside a closed region C. Simply stated, it states that the time rate of change of circulation in a closed fluid region must equal zero, DΓ Dt = 0 (2.42) or the total circulation of a closed fluid region is constant with time. In the case of a lifting body, the body carries some bound circulation related to the body lift. If the body is at steady state, then lift and circulation are constant with time, and Eq. (2.42) is satisfied. For a body that is not at steady state, lift and circulation are functions of time. Therfore, to satisfy Eq. (2.42), another source of equal and opposite circulation must 11 exist in the closed region. From physical observations, the additional circulation is known to be confined to a wake behind the body, Γwake, giving DΓ Dt = (Γbody +Γwake) Δt = 0 (2.43) 12 CHAPTER 3 PANEL CODES This chapter describes the solution of twodimensional potential flowfields using the SmithHess panel method [5]. The description starts with the solution of the flowfield about a nonlifting body, incorporates the Kutta condition to account for bound circulation, and then incorporates timedependent effects to solve for timedependent flowfields using the method of Basu and Hancock [3] as modified by Ardonceau [2]. The solution of the inviscid and incompressible flowfield about a nonlifting body represents the fundamental case to which a panel method can be applied. It also provides a starting point to describe the basic implementation of the panel method which will be expanded upon for the later liftingbody and timedependent solutions. 3.1 Nonlifting Body As described earlier, the inviscid and incompressible flowfield about a nonlifting body can be described by a potential field, which is the combination of the body and freestream potentials. Φ = Φbody +Φ∞ (3.1) To model the body potential, a distributed strength source sheet (of strength λ(s)) is placed along the fluidbody interface, s, as illustrated in Figure 3.1. This allows the body potential at an arbitrary point in the flow, P (x1, x2), to be computed in terms of the potential due to the source sheet. 13 Figure 3.1. Airfoil modeled with a continuous source sheet. Φbody(P) = s λ(s) 2π lnr ds (3.2) Correspondingly, of freestream flow is uniform, parallel to the x1axis, and the origin is a reference point where Φ∞(0, 0) = 0, the potential due to the freestream at point P is Φ∞(P) = qj∞xj (3.3) Substituting Eqs. (3.2) and (3.3) into Eq. (3.1), gives the total potential at point P due to both the body and freestream. Φ(P) = s λ(s) 2π lnr ds + qj∞xj (3.4) The only unknown parameter in Eq. (3.4) is the body source distribution, λ(s). However, by applying the noflow Neumann boundary condition from Eq. (2.27), qjnj = nj ∂Φ ∂xj = 0 (3.5) 14 Figure 3.2. Airfoil discretized into constant strength source elements. to the total potential at the fluidbody interface, nj ∂Φ ∂xj = s nj ∂ ∂xj λ(s) 2π lnr ds + qj∞nj = 0 (3.6) the unknown source distribution can be determined. Unfortunately, solving Eq. (3.6) for the source distribution is a nontrivial exercise for all but the simplest geometries. However, by applying geometric simplifications, determining the body source distribution as a function of body geometry and freestream conditions can be reduced to solving a set of linear equations. 3.1.1 Discretization By discretizing the continuous source distribution, shown in Figure 3.1, into a series of straight segments, or panels, as shown in Figure 3.2, Eq. (3.6) may be reduced to a set of dependent linear equations. For this discussion, each panel represents a unique distributed source element having a constant source strength along the length of the element. A further simplification is made in that the noflow boundary condition is not enforced at all locations 15 on the body. Rather, the noflow boundary conditions are applied to a single location, or collocation point, at the midpoint of each panel, as shown in Figure 3.2. By discretizing the body into N panels, numbered clockwise from panel 1 at the lower body trailingedge to panel N at the upper body trailing edge, the potential at the collocation point of any panel, panel α, can be determined as a function of freestream potential, body geometry, and panel strength distribution along the body. In this manner, the potential on panel α due to a source element on panel β and the freestream is Φαβ = λβ 2π β ln rαβ dsβ + qj∞xjα (3.7) The potential on panel α due to the entire body can be calculated using superposition. Thus, the potential on panel α due to the entire body is the sum of the potential due to the N panels on the body. Φα = N β=1 λβ 2π β ln rαβ dsβ + qj∞xjα (3.8) Applying the noflow boundary condition, Eq. (2.27), to Eq. (3.8) gives the normal velocity on panel α due to the body and freestream. qnα = N β=1 λβ 2π β ∂ ∂nβ ln rαβ dsβ + qj∞njα = 0 (3.9) As in Eq. (3.6), the source strengths in Eq. (3.9) are the unknown. However, because the parameters in the integrand of Eq. (3.9) are based strictly on body geometry, the integral can be replaced by a geometric influence coefficient, aαβ, which represents the geometric influence of panel β on panel α. aαβ = 1 2π β ∂ ∂nα ln rαβ dsβ (3.10) Using the influence coefficient method, the noflow normal condition on panel α, given previously in Eq. 3.9 becomes N β=1 (λβ aαβ) = −qj∞njα (3.11) 16 Figure 3.3. Constant Strength Panel Discretization Equation (3.11) is the basis for a set of linear equations relating the unknown panel source strengths λβ to the noflow boundary condition. This system of equations begins with the influence matrix, Aαβ, which is made up of the influence coefficients, aαβ, based only on the body geometry. Aαβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N a21 a22 ... a2N ... ... ... aN1 aN2 ... aNN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.12) The element strengths for each panel are stored in the column vector xβ. xβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.13) 17 Finally, the column vector Bα represents normal velocity components at the collocation point not induced by the body, such as the freestream normal velocity. Bα = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.14) Combined, these matrices and vectors form a system of equations Aαβxβ = Bα, or ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N a21 a22 ... a2N ... ... ... aN1 aN2 ... aNN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.15) the solution of which is trivial, or nonunique, describing the potentialflow about the nonlifting body. In physical terms, the trivial solution does not include the effects of bound circulation about the body, and therefore does not model lift. 3.2 Lifting Body To model the effects of lift and bound circulation about a body, additional constraints must be considered. To model bound circulation on the body, a set of constant strength vortex panels, each of the same strength, are added to the existing source panel discretization. Since each vortex panel has the same strength, only a single variable must be added to the set of linear equations modeling the nonlifting solution. The additional variable necessitates an additional constraint to solve for the vortex panel strength. This additional constraint is provided by the Kutta condition, which is based on observations of physical flow phenomena about a lifting body, or airfoil, with a sharp trailingedge. 18 3.2.1 Kutta Condition The Kutta condition is a means to relate possible potentialflow solutions about a body to observed physical flow characteristics, thereby generating a unique solution for the flowfield. The general definition of the Kutta condition specifies that the flow must detach from the airfoil at the airfoil trailingedge and that the trailingedge has zero loading [8]. The aforementioned potentialflow solution described for the nonlifting body possess a trailingedge singularity, thus the Kutta condition as specified can not be satisfied at the airfoil trailingedge. For a lifting body, a commonly used first approximation is employed in which the zero loading condition is enforced on the panels adjacent to the airfoil trailingedge. For a discretized airfoil, the condition of zero trailingedge loading is approximately satisfied by specifying equal pressure on the airfoil upper and lower trailingedge panels. The unsteady Bernoulli equation, Eq. (2.30), is used to relate the fluid flow on the upper, u, and lower, l, panels, giving p ρ + q2 j 2 + ∂Φ ∂t l = p ρ + q2 j 2 + ∂Φ ∂t u (3.16) For a steadystate flow, the timedependent potential terms can be neglected in Eq. (3.16), and the condition of equal pressure simplifies to the specification of equal flow velocity on the airfoil upper and lower trailingedge panels. qjl = qju (3.17) Thus, Eq. (3.17) provides the additional constraint necessary to solve for the unique panel strengths on the airfoil in a steadystate flow. 3.2.2 Equations Placing vortex panels along the airfoil does not change the noflow boundary condition described in Eq. (2.27), but the potential at panel α due to panel β and the freestream must now include the influence of the discreyized vortex sheet. 19 Φαβ = λβ 2π β ln rαβ dsβ − γ 2π β θαβ dsβ (3.18) Accordingly, the potential at panel α due to the N source panels, N vortex panels, and the freestream influence along the body becomes Φα = N β=1 λβ 2π β ln rαβ dsβ − γ N β=1 1 2π β θαβ dsβ + qj∞xjα (3.19) Hence, the normal velocity on panel α due to the N panels and freestream can be writen in terms of the potential gradient normal to the body at panel α qnα = N β=1 λβ 2π β ∂ ∂nα ln rαβ dsβ − γ N β=1 1 2π β ∂ ∂nα θαβ dsβ + qj∞njα = 0 (3.20) Again, the integrand for the circulatory term in Eq. (3.20) is based solely on body geometry and therefore may be calculated as a geometric influence coefficient, bα, representing the influence of the N discretized vortex panels on panel α. bα = − N β=1 1 2π β ∂ ∂nα θαβ dsβ (3.21) Substituting the influence coefficients, Eq. (3.10) and Eq. (3.21), the noflow boundary condition gives N β=1 (λβ aαβ) + γ bα = −qj∞njα (3.22) Equation (3.22) still provides N equations, but there are now N + 1 variables (N source strengths, λα, and one vortex strength, γ) describing the potential field about the lifting body. The Kutta condition, Eq. (3.17), provides the N + 1’th condition needed to solve the linear system of equations for the source and vortex strengths. Using the noflow boundary condition to simplify the Kutta condition (i.e. all flow on the trailingedge panels must be tangential) the tangential flow velocity on panel α can calculated 20 in terms of the potential gradient along the body. The tangential flow velocity on panel α due to panel β and the freestream is therefore qsαβ = λβ 2π β ∂ ∂sβ ln rαβ dsβ + γ 2π β ∂ ∂sα θαβ dsβ + qj∞sj (3.23) giving a tangential flow velocity on panel α due to all N body panels of qsα = N β=1 λβ 2π β ∂ ∂sα ln rαβ dsβ + γ N β=1 1 2π β ∂ ∂sα θαβ dsβ + qj∞sj (3.24) Examining Eq. (3.24), two new influence coefficients are introduced, cαβ, the tangential flow component along panel α due to source panel β cαβ = 1 2π β ∂ ∂sα ln rαβ dsβ (3.25) and dα, the tangential flow component along panel α due to the N body vortex panels. dα = N β=1 1 2π β ∂ ∂sα θαβ dsβ (3.26) Rewriting the steadystate Kutta condition, Eq. (3.17), in terms of the geometric influence coefficients, qsl = N β=1 (λβ c1β) + γ d1 = N β=1 (λβ cNβ) + γ dN = qsu (3.27) and rearranging to position the terms on the left hand side, n j=1 (λj (cnj − c1j)) + γ (dn − d1) = 0 (3.28) gives the Kutta condition in a suitable form to incorporate into the system of linear equations, Eq. (3.15). 21 Rewriting Eq. (3.15) to include the vortex influence and the Kutta condition, the Aαβ matrix becomes, Aαβ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N b1 a21 a22 ... a2N b2 ... ... ... ... aN1 aN2 ... aNN bN (cN1 − c11) (cN2 − c12) ... (cNN − c1N) (dN − d1) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.29) the xβ vector becomes, xj = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN γ ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.30) and the Bα vector becomes, Bi = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ −qj∞nj1 −qj∞nj2 ... −qj∞njN 0 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.31) Combining Eqs. (3.29), through (3.31) gives the linear system of equations which model the flowfield about the lifting body, ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1N b1 a21 a22 ... a2N b2 ... ... ... ... aN1 aN2 ... aNN bN (cN1 − c11) (cN2 − c12) ... (cNN − c1N) (dN − d1) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λN γ ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ −qi∞ni1 −qi∞ni2 ... −qi∞niN 0 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.32) 22 providing a unique solution to the flowfield which includes the effects of lift and bound circulation in the solution. 3.3 TimeDependent Solutions The nonlifting and lifting body solutions described in Sections 3.1 and 3.2 provide methods to model steadystate flowfields about a body in a uniform freestream flow. If the body is in motion relative to the freestream, or if the freestream includes perturbations about its timeaverage mean, assumptions neglecting timedependent terms are no longer valid and a timedependent solution methodology must be found. The basic formulation of a timedependent solution is similar to that of the lifting body solution; i.e. the body is discretized using source and vortex panel discretization, the noflow boundary condition provides N linear equations, and the Kutta condition provides the one additional relation necessary to formulate a unique liftingbody solution. The difference between the timedependent and steadystate solutions is the application of the Kutta condition and the incorporation of a model to account for the airfoil wake. The following method describes a solution for the timedependent flowfield about an airfoil in motion relative to the influence of an otherwise uniform freestream flow. 3.3.1 Frame of Reference The choice of coordinate system and reference frame determine the complexity of the mathematical model. For this discussion, the noflow boundary condition is calculated in a bodyfixed coordinate system that is allowed to translate and pitch in the global reference frame with velocity qjrel and pitch rate Ω. 3.3.2 Wake In a viscous solution for attached flow over an airfoil, a low energy boundary layer along the airfoil is shed into the freestream flow from the airfoil trailing edge to form the airfoil wake. The wake represents a pertubation of the freestream flow aft of the airfoil due to the 23 flow about the airfoil. The wake is significant because it can have a profound influence on the flow about the airfoil, even as the wake convects with the freestream. Since viscous effects are neglected in a potentialflow solution, an airfoil wake must be modeled in a way representing the effect of shed bound circulation to satisfy the Kelvin theorm. As such, the wake is often called a “time history” because it represents the change in bound circulation on the airfoil with time. Because the inviscid wake represents the change in bound circulation on the airfoil with time, or the shed circulation, it is possible to model the wake using a set of discrete vortex elements. Basu and Hancock [3] use a set of point vortices and a single constant strength vortex panel to model the inviscid timedependent wake. The strength and orientation of the wake vortex panel, or wake panel, play a key roll in satisfying the timedependent Kutta condition (explained in detail later in this chapter). The the strength of the wake panel is dependent on the amount of circulation shed by the airfoil between time steps, and wake panel orientation is determined by the Kutta condition. After the Kutta condition has been satisfied, calculations necessary to determine the timedependent flowfield solution have been performed, and any necessary post solution calculations have been completed, all wake vortices are convected with the local flow in preparation for the next time step. The wake panel is not convected with the local flow, however, rather the wake panel is replaced by a single point vortex of strength equal to the shed circulation from the previous timestep. This new point vortex is then allowed to convect with the local flow. In this manner, the wak panel and point vortices model shed circulation from the airfoil, which in turn can influence the flow about the airfoil. 3.3.3 Unsteady Kutta Condtion The specification of the timedependent Kutta condition is similar to the steadystate specification described in Section 3.2.1. The difference is in the application of the timedependent Kutta condition, which can no longer neglect timedependent terms in the un 24 steady Bernoulli equation relating pressures on the upper and lower airfoil trailingedge panels. The inclusion of timedependent terms means that the timedependent Kutta condition is a quadratic equation which must be solved iteratively. Two applications of the timedependent Kutta condition are described below, the method of Basu and Hancock [3], and the modification of that method by Ardonceau [2] used in the current investigation. 3.3.3.1 BasuHancock Basu and Hancock propose that “...there is no definitive statement of the Kutta condition for a steady airfoil, each mathematical model requiring its own consistent ‘Kutta’ condition to ensure a unique solution...” [3] Based on that statement, the assumption that the flow separates from the airfoil at the airfoil trailingedge, zero loading exists across the shed vorticity at the trailingedge, and zero loading occurs across the trailingedge elements of the airfoil, Basu and Hancock propose the folowing mathematical model for the Kutta condition. This model determines the orientation, θk, length, Δk, and strength, (γw)k, of the wake panel at time tk. Beginning with the unsteady Bernoulli equation applied at the airfoil trailingedge panels p ρ + q2 j 2 + ∂Φ ∂t l = p ρ + q2 j 2 + ∂Φ ∂t u (3.33) and specifying equal pressure at the trailingedge, pl − pu ρ = q2 ju 2 − q2 jl 2 + ∂Φu ∂t − ∂Φl ∂t = 0 (3.34) a quadratic relation develops for the flow velocity on the upper and lower airfoil trailingedge panels. Because the velocity relation is not linear, an iterative solution is necessary to determine the orientation and strenth of the wake panel which satisfies the Kutta condition. 25 Using the Kelvin theorem, Eq (2.42), the rate of change of circulation about the airfoil must be balanced by the rate of change of the shed circulation in its wake, Δk (γw)k Δt = −∂Γ ∂t = −Γk−1 − Γk Δt (3.35) or the change in circulation about the airfoil from tk−1 to tk must be balanced by an equal and opposite circulation about the wake panel. Δk (γw)k = Γk − Γk−1 (3.36) The rate of change of potential across the airfoil trailingedge is related to the rate of change in circulation by ∂ (Φl − Φu) ∂t = ∂Γ ∂t (3.37) Therefore, substituting Eq. (3.37) into Eq. (3.34) relates the upper and lower trailingedge velocities to the rate of change of circulation about the airfoil. q2 ju − q2 jl 2 + Γk − Γk−1 tk − tk−1 = 0 (3.38) Substituting Eq. (3.36) into Eq. (3.38) gives the circulation strength about the wake panel in terms of trailingedge panels velocities and wake panel length. (γw)k = (tk − tk−1) (q2 l − q2u ) 2Δk (3.39) Wake panel orientation is determined by local velocity on the wake panel, neglecting the effect of the wake panel on itself, tan θk = (q1w)k (q2w)k (3.40) and wake panel length is proportional to the magnitude of the local velocity and the time step. Δk = (qjw)k (tk − tk−1) (3.41) 26 3.3.3.2 Ardonceau Ardonceau proposed a modification to Basu and Hancock’s Kutta condition based on experimental studies [2]. The modified solution method is nearly identical to that of Basu and Hancock, but the wake panel geometry is altered. Instead of allowing the wake panel to change both orientation and length, the wake panel orientation is fixed along the bisector between the upper and lower trailingedge panels. θk = θu + θl 2 (3.42) The length of the Ardonceau wake panel then equals the average of the trailingedge panel velocities porportional to the time step. Δk = 1 2 (qju + qjl)k (tk − tk−1) (3.43) The calculation of wake panel strength is the same as Eq. (3.39). 3.3.4 Method of Solution Regardless of the mathematical formulation of the unsteady Kutta condition, the solution methods are the same. As in the steadystate solutions, the N source strengths, one vortex strength, and freestream along the body are related through the noflow boundary condition which gives a system of N linear equations. As outlined above, however, the Kutta condition becomes a quadratic relation in an unsteady flow which must be solved using an iterative techique. The noflow boundary condition for the timedependent solution also includes induced velocity terms due to body motion relative to the freestream and induced velocity terms due to the airfoil wake. Modifying Eq. (3.22) to include the effects of body rotation, qjrotation = Ω× rα (3.44) 27 body translation, qjtranslation = qjrel (3.45) and the influence of the wake panel and point vortices, qjwake = γw bαN+1 + k−1 β=1 Γβ ∂ ∂nα θαβ 2π (3.46) the time dependent noflow relation becomes N β=1 (λβ aαβ) + γ N β=1 bαβ + γw bαN+1 + k−1 β=1 Γβ ∂ ∂nα θαβ 2π + (qj∞ +Ω× rα + qjrel) njα = 0 (3.47) Rearranging to place all nonsource terms on the right hand side gives N β=1 (λβ aαβ) = −γ N β=1 bαβ−γw bαN+1− k−1 β=1 Γβ ∂ ∂nα θαβ 2π −(qj∞ +Ω× rα + qjrel) njα (3.48) Note that Eq. (3.48) is very similar to Eq. (3.11) but with extra terms on the right hand side. Therefore, rewriting Eq. (3.14) to include the new terms of Eq. (3.48) gives Bi = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −γ N β=1 b1β − γw b1N+1 − k−1 β=1 Γβ ∂ ∂n1 θ1β 2π − (qj∞ +Ω× r1 + qjrel) nj1 −γ N β=1 b2β − γw b2N+1 − k−1 β=1 Γβ ∂ ∂n2 θ2β 2π − (qj∞ +Ω× r2 + qjrel) nj2 ... −γ N β=1 bNβ − γw bNN+1 − k−1 β=1 Γβ ∂ ∂nN θNβ 2π − (qj∞ +Ω× rN + qjrel) njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.49) Substituting Eq. (3.49) for Eqs. (3.14) in Eq. (3.15) gives a linear system of equations for the timedependent solution. 28 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ a11 a12 ... a1n a21 a22 ... a2n ... ... ... an1 an2 ... ann ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ λ1 λ2 ... λn ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ −γ N β=1 b1β − γw b1N+1 − k−1 β=1 Γβ ∂ ∂n1 θ1β 2π − (qj∞ +Ω× r1 + qjrel) nj1 −γ N β=1 b2β − γw b2N+1 − k−1 β=1 Γβ ∂ ∂n2 θ2β 2π − (qj∞ +Ω× r2 + qjrel) nj2 ... −γ N β=1 bNβ − γw bNN+1 − k−1 β=1 Γβ ∂ ∂nN θNβ 2π − (qj∞ +Ω× rN + qjrel) njN ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ (3.50) An iterative solution scheme is used to find the unique solution satisfying both the system of N linear equations and one quadratic relation. The iterative Kutta condition assumes initial values for the wake panel orientation, length, and circulation strength at the initialization of the simulation. The inital values are used in the solution of N linear equations to determine the strength of the body source elements. The calculated body source strengths are then used to recalculate the orientation, length, and circulation strength of the wake panel, and the process is repeated until the orientation, length, and circulation strength of the wake panel meet a given convergence criteria. 29 CHAPTER 4 CODE DESCRIPTION To facilitate investigations into the interaction between elastic airfoil response and arbitrary aerodynamic forcing, two components are added to the timedependent panel code described in Section 3.3. The first component is a gust model originally proposed by Basu and Hancock [3] which uses singularity elements to model the influence of a sharp edge gust. The second component is a structural model which, when coupled with the aerodynamic model, determines the airfoil responce to aerodynamic forcing. This chapter describes the gust and structural models, as well as their implementation and integration into the unsteady panel code. 4.1 Frame of Reference The source and vortex elements modeling the airfoil, wake, and freestream perturbations are tracked in a Lagrangian reference frame. The origin of this reference frame is located at the leading edge of the undisturbed airfoil, with the airfoil trailing edge lying on the x1 axis. 4.2 Gust Model The gust model uses singularity elements to induce velocity perturbations about an otherwise uniform freestream flow. Because this investigation is initially interested in the effect of transverse velocity perturbations, or perturbations perpendicular to the timeaveraged freestream flow, the gust is modeled by a set of vortex sheets. For this discussion, a vortex sheet will be defined as a collection of vortex elements, each sharing at least one end point with a neighboring vortex element. Collectively, these vortex elements produce a continuous 30 Figure 4.1. Influence of a Vortex sheet located in the Freestream flow compared to the Freestream influence. vorticity “sheet” that convects with the freestream flow. Each vortex sheet possess a finite amount of bound circulation that remains constant with time, and is initially distributed evenly along the length of the sheet. The influence of a single gust sheet in the freestream flow is shown in Figure 4.1. By placing gust sheets into the flowfield prior to initialization of the simulation, and specifying a timeinvariant total circulation about each gust sheet, Kelvin’s Theorm is implicitly satisfied. Thus, the Kutta condition discussed in Section 3.3 remains valid. 4.2.1 Deformation The key to properly modeling the transverse gust, such that the gust responds to the influence of the airfoil and its wake, lies in modeling gust convection. To allow the gust sheet to deform and react to the influence of the airfoil, wake, and other elements in the flowfield, each gust sheet is discretized into a finite number of panels, or elements. At the end of each computational time step, the gust sheet is convected by propagating the endpoints of each gust element with the local fluid flow. In this manner, the gust sheet is alowed to deform 31 Figure 4.2. Deformation of a vortex sheet approaching the airfoil leading edge. due to local velocity gradients in the flow. Figure 4.2 illustrates the deformation of a gust sheet as it approaches the airfoil leading edge. Because the total bound circulation about each gust element is timeinvariant, the influence of each gust element on the surrounding fluid flow is a function of element length. 4.2.2 AirfoilGust Interaction Each gust sheet initialized upstream of the airfoil eventually encounters the airfoil as it convects with the freestream flow. Since the gust sheet provides aerodynamic forcing for forcedresponse simulations, proper modeling of airfoilgust interaction is a critical aspect of the overall gust model. To properly model gust sheet influence on the airfoil, the gust sheet must propagate around the airfoil, not propagate through the airfoil. Therefore, the continous gust sheet must be “split” when the gust sheet encounters the forwardmost edge of the airfoil, allowing one section of the sheet to convect across the upper airfoil surface and the other to convect 32 Figure 4.3. Case One: Gust element straddling the airfoil. along the lower airfoil surface. Techniques used to determine if and when a gust sheet must be split, and methods used to split each gust sheet are described below. Two distinct cases can arise when a gust sheet reaches the airfoil. Case One involves a gust element straddling the airfoil leading edge. In this case, one endpoint attempts to convect above the airfoil while the other endpoint convects below, as illustrated in Figure 4.3. To maintain the noflow boundary condition, the gust element straddling the airfoil must be split into two separate elements, one ending on the upper airfoil surface and one ending on the lower airfoil surface. The two “new” elements must also posses a combined bound circulation equal to the bound circulation of the original gust element to satisfy Kelvin’s theorem. Case Two involves a gust element, or pair of elements sharing an endpoint, where the endpoint attempts to convect into the airfoil interior, as illustrated in Figure 4.4. To maintain the noflow boundary condition, the erroneous endpoint must be to the either the upper or lower airfoil surface. In this case, no new gust elements are created, and each affected elements retains its original bound circulation. 33 Figure 4.4. Case Two: Gust element endpoint convected into the airfoil. In each case, once the erroneous elements have been relocated to the airfoil surface, the element endpoints lying on the airfoil surface are convected along the airfoil at the surface tangential velocity, until the gust element propagates past the airfoil trailing edge. 4.2.2.1 Determining Gust Element Condition To ensure a gust sheet does not breach the airfoil interior, the following conditions are checked for each gust element after it is convected in preparation for the next time step. 1. Does the gust element currently terminate on the airfoil surface? 2. Is either element endpoint located between the airfoil leading and trailing edges in the x1 direction? 3. Is one element endpoint located above the airfoil while the other is located below the airfoil in the x2direction? 4. Is either gust element endpoint located inside the airfoil surface? 34 Based on the four conditions, the state of the gust element with respect to the airfoil can be determined. If condition 1 is true, the gust element must be convected along the airfoil at the surface tangential velocity instead of the local flow velocity. If condition 1 is false and conditions 2 and 3 are true, the gust element is an example of Case One and must be split. If condition 1 is false and conditions 2 and 4 are true, the gust element is an example of Case Two and the element endpoint must be relocated to the airfoil surface. If conditions 1 through 4 are false, the gust element is located in the freestream and no gust sheet modifications are required. 4.2.2.2 Case One Case One involves splitting a gust element straddling the airfoil, as illustrated in Figure 4.3, and determining the bound circulation about the split gust elements. Because the airfoil may be at some arbitrary orientation relative to the timeaveraged freestream flow, the airfoil leadingedge node may not be the airfoil node the gust sheet first encounters, therefore the term “forward most” edge, or node, will be defined as the node closest to the gust when the gust impacts the airfoil. In addition, depending on airfoil orientation and the influence of the freestream (including the gust), the upstream stagnation point on the airfoil may not correspond to either the forwardmost airfoil node or the airfoil leading edge. This distinction is subtle, as illustrated in Figure 4.5. The variance in x1 location between the leading edge and stagnation point may only be a few hundredths of a chord length, but the influence of this variance on resulting airfoil forcing can be significant. For example, consider the case where a gust element is split about the airfoil leading edge at time tk+1, but the upstream stagnation point on the airfoil does not correspond to the airfoil leading edge. If the upstream stagnation point is located on the lower airfoil surface, the gust element ending on the lower surface between the airfoil leading edge and the upstream stagnation point will convect towards the airfoil leading edge at time tk+2 instead of towards the airfoil trailing edge, as desired. This process is depicted in Figures 4.6 and 35 Figure 4.5. Airfoil leading edge vs. the airfoil forwardmost node. 4.7. In fact, the lower gust element will eventually propagate around the leading edge and convect towards the airfoil trailing edge along the upper surface. This will stretch the gust element through the airfoil, invalidating the noflow boundary condition. A similar circumstance occurs if the gust element is simply split about the upstream stagnation point. For example, if the upstream stagnation point does not correspond to the leading edge, but rather lies on the lower airfoil surface, the gust element propagating above the airfoil will stretch through the airfoil and end on the lower airfoil surface, as illustrated in Figure 4.8. The upper gust element will eventually propagate around the airfoil leading edge before it propagates towards the trailing edge along the upper airfoil surface, as desired, but this gives the gust element an undue influence on the airfoil as it propagates around the airfoil leading edge and will invalidate the noflow boundary condition. Because the upstream stagnation point and the forwardmost node both exhibit large influences on the gust element, gust elements straddling the airfoil leading edge are split about both the forwardmost airfoil node and the upstream stagnation point. In this manner, if the upstream stagnation point is located on the lower airfoil surface, the lower gust element will 36 Figure 4.6. Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+1. Figure 4.7. Gust element split about the leading edge with the upstream stagnation point on the lower airfoil surface at time tk+2. 37 Figure 4.8. Gust element split about the upstream stagnation point. convect along the airfoil from the upstream stagnation point while the upper gust element will convect along the airfoil from the forwardmost airfoil node, and visa versa for an upstream stagnation point located on the upper airfoil surface. In most cases, this distinction is negligible, but the method ensures that split gust elements will not convect towards the airfoil leading edge, or stretch through the airfoil surface. 4.2.2.3 Implementation As mentioned, once it has been determined that a gust element straddles the airfoil, the element must be split into two “new” elements, one convecting above the airfoil and one convecting below the airfoil. Because the unsteady panel code models the flowfield using discrete time steps, it is unlikely that the instant a gust element impacts the forwardmost airfoil node will correspond exactly to a panel code time step. Therefore, an interpolation routine is employed to accurately determine the instant in time a gust element impacts the forwardmost airfoil node, as illustrated in Figure 4.9. It is necessary to know this time because, for example, if a gust element impact occurs at midway between timesteps, the 38 Figure 4.9. Interpolation to determine the time of GustAirfoil impact. element should be convected along the airfoil during the remaining amount of the time step after being split. In addition to accurately determining the instant in time that a gust element impacts the airfoil, the interpolation routine also provides information regarding what percentage of the original gust element should convect above and below the airfoil. Knowing these percentages is necessary so proper fractions of the original bound circulation can be assigned to each “new” gust element, thereby maintaining a constant total circulation in the flow. 4.2.2.4 Case Two Case Two involves a gust element, or pair of elements, possessing an endpoint that convects into the closed airfoil surface, as illustrated in Figure 4.4. This is the less common of the two cases, and for a simulation with a suitably small time step only occurs if the initial gust sheet contains an element possessing an endpoint close to the x1 axis. Therfore, in an effort to simplify the panel code, this case is controled through well considered initial discretization of the gust sheet. 39 Figure 4.10. Gust element convection along the upper airfoil surface. 4.2.3 Convection For a gust element ending on the airfoil surface, the endpoint on the surface is convected at the surface tangential velocity instead of the local flow velocity. Because the airfoil itself is discretized into a set of discrete panels and the noflow boundary condition is enforced only at each panel midpoint, a gust element endpoint convected at the local flow velocity could convect into the airfoil surface, or off the airfoil surface into the freestream flow. Basu and Hancock [3] calculated the surface tangential velocity at the gust element endpoint by interpolating tangential velocities across adjacent airfoil panels. The interpolated surface tangential velocity value was then multiplied by the local time step to find the distance the element endpoint should convect along the airfoil surface. This method provides a good first aproximation for coarse airfoil discretizations, but fails for finely discretized airfoils in locations where a large velocity gradient exists between adjacent panels, such as at the airfoil leading edge. To acount for large tangential velocity gradients, an alternate method of convecting a gust element along the airfoil surface has been developed. This alternate method estemates the 40 amount time nessisary to convect the gust endpoint along a surface panel based on the length of the surface panel and the surface tangential velocity at the panel midpoint. The estimated time to convect the gust element endpoint to the end of the surface panel is compared to the amount of time remaining in the computational timestep. Based on whether the estimated time is greater than the remaning time step, a decision is made to convect the endpoint a fractional distance along the surface panel, based on the surface tangential velocity and the remaining time step, or to convect the endpoint to the end of the current surface panel, and repeat the time estimation on the next surface panel. For example, to convect the gust element endpoint initially located at some location along Panel a, as depicted in Figure 4.10, the distance between the gust endpoint and the downstream node of Panel a is used with the surface tangential velocity at the midpoint of Panel a to estimate the amount of time necessary to convect the gust endpoint to the downstream node of Panel a. If the estimated time to convect the gust endpoint to the end of Panel a is less than the local time step, Δt, or for convenience, the time remaining, tr, then the gust endpoint is relocated to the downstream airfoil node shared by Panels a and b, and the estimated time is subtracted from tr. In this manner, using the lengths of Panels b, c, and d, along with their respective tangential velocities, the time necessary to convect the gust endpoint across Panels b, c, and d is estimated to be greater than tr, but the time necessary to convect the gust endpoint across only Panels b and c is less than tr. Thus, the gust endpoint is relocated to the shared airfoil node between Panels c and d, and the estimated time to convect the gust endpoint across Panels b and c is subtracted from tr. Because the time necessary to convect across Panel d is greater than tr, the gust endpoint is relocated a fractional distance along Panel d, as determined using the surface tangential velocity at the midpoint of Panel d and tr. 41 Figure 4.11. Gust element convection along the upper airfoil surface. 4.2.4 Gust Influence on the Airfoil Since the gust sheet is composed of singularity elements, each with an influence proportional to 1/r, the influence of a gust element ending on the airfoil depends on the proximity of the element endpoint to an airfoil collocation point. If a gust element ends on a collocation point, r approaches zero and the influence of that gust element becomes infinite. This skews the flowfield solution in a nonphysical manner. Basu and Hancock [3] prevented this possibility by replacing the each gust element ending on the airfoil with a pair of “imaginary” elements, illustrated in Figure 4.11. The two imaginary gust elements share the freestream endpoint with the original gust element, but instead of terminating at some location along an airfoil panel, airfoil panel a, with the original gust element, the imaginary element pair terminate at corresponding endpoints of airfoil panel a. The imaginary elements share the bound circulation of the original gust element in a manner dependent on the location of the original element endpoint on panel a. As such, the influence of the gust continues to propagate across the airfoil surface but the possibility of discontinuities arrising due to a gust element coexisting with an airfoil colocation point is eliminated. 42 Figure 4.12. Pitching and Plunging Airfoil 4.3 Free Response The inclusion of an airfoil structural model enables the unsteady panel code to model timedependent airfoil response due to arbitrary and selfinduced aerodynamic forcing. 4.3.1 Model The airfoil structural model is a two degree of freedom (TDOF) springmass system allowing coupled airfoil motion in rotation and translation, or pitch and plunge. Figure 4.12 shows a basic schematic detailing parameters important to the model. The equations governing twodimensional body motion in terms of sectional characteristic and generalized external forces are m¨h + Sα ¨α + mω2h h =Qh (4.1a) Sα¨h + Iα ¨α + Iαω2α α =Qα (4.1b) For a thin airfoil, the generalized external forces correspond to aerodynamic lift and moment about the elastic axis. 43 Qh = − L (4.2a) Qα =My (4.2b) For compatibility with the developed unsteady panel code, which calculates nondimensional forces and moments through integration of instantanious surfacepressure coefficents, Eq. (4.1) is nondimensionalized with respect to chord, freestream velocity, time, and mass. The resulting nondimensional equations of motion are then rewritten in terms of the nondimensional sectional characteristics, such as density ratio, μ, radius of gyration, rα, static imbalance, xα, reduced bending frequency, kh, reduced pitching frequency, kα, normalized plunge, ˆh , normalized pitch, ˆα, and nondimensional time, ˆt. μˆh + xαμ 2 ˆα + μk2h ˆh = 2 π Cl (4.3a) xαμ 2 ˆh + r2α μ 4 ˆα + r2α μk2α 4 ˆα = 2 π Cmy (4.3b) Expressing Eq. (4.3) in matrix notation, M X + K X = F (4.4) where 44 M = ⎡ ⎢⎣ μ xαμ 2 xαμ 2 r2α μ 4 ⎤ ⎥⎦ (4.5a) K = ⎡ ⎢⎣ μk2h 0 0 r2α μk2α 4 ⎤ ⎥⎦ (4.5b) X = ⎧⎪⎨ ⎪⎩ ˆh ˆα ⎫⎪⎬ ⎪⎭ (4.5c) F = 2 π ⎧⎪⎨ ⎪⎩ Cl Cmy ⎫⎪⎬ ⎪⎭ (4.5d) The second derivative in Eq. (4.4) can be isolated on the LHS, X = M −1 F − M −1 K X (4.6) allowing the equations of motion to be writen as a set of coupled first order ordinary differential equations. X = Y (4.7a) Y = M −1 F − M −1 K X (4.7b) Airfoil orientation and position at time tk+1 is determined by solving Eq. (4.7) with a fourthorder RungeKutta method using nondimensional aerodynamic forces computed at time tk. 4.3.2 Solution The aerodynamic solution and TDOF structural model are coupled directly in the developed unsteady panel method to calculate free and forced response of an arbitrary thin airfoil. The nondimensional forces and moments calulated at each time step are used as inputs to the structural model, predicting airfoil orientation and position at the next timestep. The 45 new airfoil position and orientation are used to calculate new nondimensional aerodynamic forces, and the process is repeated. 4.4 Forced Response By coupling the gust model described in Section 4.2 and the structural model described in Section 4.3, airfoil responce to arbitrary gust induced forcing can be modeled. As will be shown in Chapter 5, the influence of multiple gust sheets can be superimposed to model periodic freestream disturbance having arbitrary shapes, frequencies, and amplitudes. Thus, airfoil responce due to external forcing can be systematically studied by varying the characteristics of the freestream gust. 46 CHAPTER 5 CODE VERIFICATION To verify the accuracy and applicability of the developed panel code, a set of test cases were examined. These test cases compare unsteady panel code simulations with fundamental problems in unsteady aerodynamics having known analytical or computational solutions. In this manner, the accuracy and applicability of the panel code is established prior to its extension to problems of interest not having known solutions. 5.1 Wagner The Wagner problem, one of the fundamental problems in unsteady aerodynamics, explores the lift response of a flat plate to a flowfield which is instantaneously accelerated from one equilibrium state to another. The problem demonstrates the effect of body wake development on lift and moment during transition between equilibrium states. 5.1.1 Description Consider a stationary flat plate, or airfoil of infinitesimal thickness, at some angular orientation relative to a freestream flow, α0, illustrated in Figure 5.1. At time t < 0, the magnitude of the freestream relative to the flat plate is zero, q∞ = 0. Since the noflow boundary condition is implicitly satisfied, the body produces zero lift, and perhaps more importantly, carries zero bound circulation. At time t = 0, the freestream instantaneously accelerates to a finite nonzero velocity, q∞ = c. By applying the unsteady Kutta condition and noflow boundary condition discussed in Section 3.3, a lifting solution can be found for the flat plate. 47 Figure 5.1. Flat plate at time t = 0 It should be recalled from Section 3.3 that the body wake for an inviscid solution represents shed bound vorticity from the body which is necessary to satisfy Kelvin’s theorem. As such, the shed vorticity magnitude in the wake at time t = 0 equals the magnitude of the bound circulation change about the body, but in the opposite direction. The shed circulation caused by the flowfield transition between equilibrium states is often called a “starting vortex” because the magnitude of this vortex is significantly greater then the rest of the wake. Shed vorticity in the wake produces an aerodynamic downwash on the body, influencing the noflow boundary condition. Wake influence on lift is normally of a small magnitude relative to the freestream and the relative body motion, but in the case of a starting vortex where the shed circulation magnitude is on the same order as the bound circulation about the body, the wakeinduced downwash suppresses lift generation on the body. As such, the starting vortex significantly influences lift development on the body until the starting vortex propagates into the far field. 48 5.1.2 Solution By modeling the induced body wake as a continuous vortex sheet of varying strength, originating at the body trailing edge and oriented parallel to the freestream flow, Wagner [4] developed a timeaccurate solution for lift on an instantaneously accelerated flat plate. L = 2πbρq2α0φ (s) (5.1) This solution depends on a modified Bessel function known as the Wagner function, φ (s). φ (s) = 1 2πi ∞ −∞ C (k) k eiks dk (5.2) An approximate representation [4] of the Wagner function has been computed as, φ (s) ≈ 1 − 0.165e −0.0455s − 0.335e −0.3s (5.3) the solution of which is shown in Figure 5.2, along with the solution to the approximate Kussner function described in Section 5.4. 5.1.3 Comparison To assist verification of the developed panel code, lift solutions for thin symmetric airfoils computed using the panel code are compared the Wagner lift solution for a flat plate. Panel code solutions were obtained for instantaneously accelerated NACA 0006, 00010, and 0014 airfoils oriented at α0 = 1, 2, and 4 deg relative to a uniform freestream in the x1 direction. Solutions were computed using nondimensionalized time steps of 0.005, 0.075, and 0.010, corresponding to 4000, 3000, and 2000 computational iterations, respectively. Calculated lift coefficients for each simulation were normalized by corresponding steadystate lift values, allowing a comparison to the Wagner function, Eq. (5.3). Note that differing fundamental assumptions between the panel code and the Wagner solution affect direct comparison of the results. For example, the Wagner solution assumes 49 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 1.2 s φ(s) Approximate Wagner Approximate Kussner Figure 5.2. Solutions for the approximate Wagner function, Eq. (5.3), and the approximate Kussner function, Eq. (5.19) the body wake is a continous vortex sheet convecting at the mean freestream velocity, while the panel code discretizes the wake into a set of discrete vortices convecting at the local velocity. Also, the Wagner solution models a flat plate with negligible thickness, while the panel code models a thin symmetric airfoil. Figure 5.3 compares panel code solutions for airfoils of different thicknesses to the Wagner function. The panel code solutions in Figure 5.3 are computed for NACA 0006, 0010, and 0014 airfoils at α0 = 1 deg using a normalized time step of 0.005. As airfoil thickness decreases, the panel solutions approach the Wagner function. Figure 5.4 compares panel code solutions for a single airfoil at several orientation angles to the Wagner function. Panel code solutions in Figure 5.4 are computed for a NACA 0010 airfoil at α0 = 1, 2, and 4 deg using a normalized time step of 0.010. As Figure 5.4 shows, airfoil orientation does not have a discernable effect on normalized lift. Figure 5.5 compares panel code simulations for a single airfoil thickness and orientation but at varying normalized time steps. Panel code solutions in Figure 5.5 are computed for a 50 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t NACA 0006 NACA 0010 NACA 0014 Wagner Figure 5.3. Normalized lift for the NACA 0006, 0010, and 0014 airfoils at α0 = 1 deg using a normalized time step of 0.005 compared to Eq. (5.3) 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t α = 1 α = 2 α = 4 Wagner Figure 5.4. Normalized lift on a NACA 0010 at α0 = 1, 2, and 4 deg using a normalized time step of 0.010 compared to Eq. (5.3) 51 0 1 2 3 4 5 6 7 8 9 10 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 Cl/Clss t dt = 0.005 dt = 0.075 dt = 0.010 Wagner Figure 5.5. Normalized lift on a NACA 0010 at α0 = 2 deg computed using nondimensionalized time steps of 0.005, 0.075, and 0.010 compared to Eq. (5.3) NACA 0010 airfoil at α0 = 2 deg using nondimensionalized time steps of 0.005, 0.075, and 0.010. The panel code solutions show no significant dependence on the selected normalized time steps. Since neither time step nor orientation significantly affects the panel code solutions, differences between the panel code and Wagner solutions can be attributed primarily to airfoil thickness effects. Despite their differences, however, good overall agreement exists between the two lift solutions. 52 Figure 5.6. Notation used to describe the Theodorsen pitching and plunging flat plate 5.2 Theodorsen The Theodorsen problem, or the problem of a periodically pitching and plunging airfoil in an otherwise steady uniform freestream flow, demonstrates the effect of body motion and timedependent wake on unsteady lift and moments. 5.2.1 Description For an airfoil translating and rotating relative to an otherwise uniform freestream flow, induced flow perturbations near the airfoil surface due to its relative motion can be significant. For an airfoil undergoing periodic translational and rotational relative motion, the influence of the induced surface flow perturbation is a function of motion frequency and amplitude. In addition to motion induced flow perturbations, wake circulation will induce velocity perturbations which also influence the unsteady airfoil lift and moment as a function of motion frequency and amplitude. Figure 5.6 illustrates common notation used to describe the Theodorsen problem. 53 5.2.2 Solution Using conformal mapping techniques and a wake model similar to that employed in the Wagner solution, Theodorsen developed an analytic solution for lift and moment on a flat plate undergoing periodic pitching and plunging. The solution relates induced lift and moment on a flat plate to reduced frequency of the relative motion. L = LNC + LC C (k) (5.4) M = MNC +MC C (k) (5.5) Note that the Theodorsen lift and moment solutions separate noncirculatory terms, the irrotational component of lift and moment due to body motion, LNC = πρb2 ! ¨h + Uα˙ − baα¨ " (5.6) MNC = πρb2 ba¨h − Ub 1 2 − a α˙ − b2 1 8 − a2 ¨α (5.7) from circulatory terms, the rotational component of lift and moment necessitated by the Kutta condition. LC = 2πρUb ˙h + Uα + b 1 2 − a α˙ (5.8) MC = 2πρUb2 a + 1 2 ˙h + Uα + b 1 2 − a α˙ (5.9) The distinction between noncirculatory and circulatory terms is of importance because cirulatory terms depend on motion reduced frequency, as related through the Theodorsen function. C (k) = H2 1 (k) H2 1 (k) + iH2 0 (k) (5.10) 54 Combining Eqs. (5.6) through (5.8) with Eqs. (5.4) and (5.5) produces timedependent Theodorsen lift and moment equations for a flat plate undergoing periodic pitching and plunging relative to the freestream flow. L = πρb2 ! ¨h + Uα˙ − baα¨ " + 2πρUb ˙h + Uα + b 1 2 − a α˙ C (k) (5.11) My =πρb2 ba¨h − Ub 1 2 − a α˙ − b2 1 8 − a2 ¨α + 2πρUb2 a + 1 2 ˙h + Uα + b 1 2 − a α˙ C (k) (5.12) 5.2.3 Comparison To further assist verification of the unsteady panel code, namely the effects of relative body motion, computed solutions for thin symetric airfoils undergoing periodic pitching, periodic plunging, and periodic pitching and plunging are compared to the corresponding Theodorsen solution for a flat plate. As with comparisons to the Wagner function, the effects of thickness and wake model limit direct comparison between the panel code and analytic solutions. 5.2.3.1 Pure Pitching Panel code solutions for NACA 0006, 0010, and 0014 airfoils pitching relative to the freestream flow were computed for reduced frequencies of k = 0.25 and 0.75 and amplitudes of ¯α = 1, 2, and 4 deg about the airfoil quarterchord location. Timedependent lift and moment for a flat plate undergoing similar motion were also computed using Eqs. (5.11) and (5.12). Figures 5.7 through 5.9 demonstrate the effect of airfoil thickness on panel code lift and moment solutions, as compared to the Theodorsen solution. Panel code solutions in Figures 5.7 through 5.9 were computed for NACA 0006, 0010, and 0014 airfoils pitching at a reduced frequency of k = 0.25 and an amplitude of ¯α = 2deg. 55 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.7. Cl vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg For pure pitching at small reduced frequencies, airfoil thickness exhibits a small influence on the phase between the panel code and Theodorsen lift solutions as well as the amplitude ratio of the two solutions. However, both amplitude and phase of the panel code solution approach the Theodorsen solution as airfoil thickness decreases. Figures 5.10 through 5.12 shows the relative agreement between the panel code lift and moment to the Theodorsen solution, for a range of pitching amplitudes. Panel code solutions for Figures 5.10 through 5.12 were computed for a NACA 0010 airfoil pitching at a reduced frequency of k = 0.25, and pitching amplitudes of ¯α = 1, 2, and 4 deg. For pure pitching at a constant reduced frequency, pitching amplitude does not appear to exhibit a significant influence on either the phase between the panel code and Theodorsen lift solutions or the lift ratio between the two solutions. The lift ratio, computed as the maximum panel code lift coefficent divided by the maximum Theodorsen lift coefficient, remains constant around 1.08 for the pitching amplitudes computed. 56 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmle t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.8. Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmea t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.9. Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitude of ¯α = 2deg 57 0 2 4 6 8 10 12 14 16 18 20 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Cl t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.10. Cl vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 Cmle t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.11. Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg 58 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen α = 1 UVPM α = 1 Theodorsen α = 2 UVPM α = 2 Theodorsen α = 4 UVPM α = 4 Figure 5.12. Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at a reduced frequency of k = 0.25 and amplitudes of ¯α = 1, 2, and 4 deg Figures 5.13 through 5.15 show the relative agreement of the panel code lift and moment solutions to the Theodorsen solution, for a range of reduced frequencies. Panel code solutions in Figures 5.13 through 5.15 are a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with a pitching amplitude of ¯α = 2 deg. For pure pitching at a constant amplitude, reduced frequency does not appear to exhibit an influence on the phase between the panel code and Theodorsen lift solutions, but does appear to influence the amplitude ratio between the two solutions. It appears that the pase between the panel code and Theodorsen lift solutions remains constant as reduced frequency varies, however, the amplitude ratio between the panel code and Theodorsen lift solution decreases as reduced frequency increases. 59 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.13. Cl vs. Time for a NACA 0010 airfoil pitching at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cmle t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.14. Cmle vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 60 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.15. Cmea vs. Time for a NACA 0010 airfoil pitching about the quarterchord at reduced frequencies of k = 0.25 and 0.75 with an amplitude of ¯α = 2deg 61 0 2 4 6 8 10 12 14 16 18 20 −0.025 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 Cl t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.16. Cl vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 5.2.3.2 Pure Plunging The unsteady panel code was used to generate solutions for NACA 0006, 0010, and 0014 airfoils plunging at reduced frequencies of k = 0.25 and 0.75 with plunging amplitudes of ¯h = 0.010, 0.025, and 0.050 relative to the mean freestream flow. Timedependent lift and moment for a flat plate undergoing similar motion were also computed using Eqs. (5.11) and (5.12). The halfchord was used to calculate moments about the elastic axis. Figures 5.16 through 5.18 demonstrate the effect of airfoil thickness on panel code lift and moment solutions, as compared to the Theodorsen solution. Panel code solutions in Figures 5.16 through 5.18 were computed for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and amplitude of ¯h = 0.025. For pure plunging at small reduced frequencies, as in the case of pure pitching, airfoil thickness exhibits a small influence on the phase between the panel code and Theodorsen lift solutions as well as the amplitude ratio of the two solutions. However, both amplitude and phase of the panel code solution approach the Theodorsen solution as airfoil thickness decreases. 62 0 2 4 6 8 10 12 14 16 18 20 −0.01 −0.005 0 0.005 0.01 Cmle t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.17. Cmle vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 0 2 4 6 8 10 12 14 16 18 20 −6 −4 −2 0 2 4 6 x 10−3 Cmea t Theodorsen UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0014 Figure 5.18. Cmea vs. Time for NACA 0006, 0010, and 0014 airfoils plunging at a reduced frequency of k = 0.25 and an amplitude of ¯h = 0.025 63 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cl t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.19. Cl vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 Figures 5.19 through 5.21 show the relative agreement between the panel code lift and moment solution to the Theodorsen solution for different plunging amplitudes. Panel code solutions in Figures 5.19 through 5.21 were computed for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050. As is the case for pure pitching, for pure plunging at a constant reduced frequency, pitching amplitude does not appear to exhibit a significant influence on either the phase between the panel code and Theodorsen lift solutions or the lift ratio between the two solutions. The amplitude ratio remains constant around 0.99 for the plunging amplitudes computed. Figures 5.22 though 5.24 show relative agreement between the panel code lift and moment solution to the Theodorsens solution at different reduced frequencies. Panel code solutions in Figures 5.22 though 5.24 were computed for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and amplitude of ¯h = 0.025. As is the case for pure pitching, for pure plunging at a constant amplitude, reduced frequency does not appear to exhibit an influence on the phase between the panel code and Theodorsen lift solutions, but does appear 64 0 2 4 6 8 10 12 14 16 18 20 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 Cmle t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.20. Cmle vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 0 2 4 6 8 10 12 14 16 18 20 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04 0.05 Cmea t Theodorsen h = 0.010 UVPM h = 0.010 Theodorsen h = 0.025 UVPM h = 0.025 Theodorsen h = 0.050 UVPM h = 0.050 Figure 5.21. Cmea vs. Time for a NACA 0010 airfoil plunging at a reduced frequency of k = 0.25 and amplitudes of ¯h = 0.010, 0.025, and 0.050 65 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.22. Cl vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 to influence the amplitude ratio between the two solutions. It appears that the pase between the panel code and Theodorsen lift solutions remains constant as reduced frequency varies, however, the amplitude ratio between the panel code and Theodorsen lift solution decreases as reduced frequency increases. 66 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cmle t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.23. Cmle vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen k = 0.25 UVPM k = 0.25 Theodorsen k = 0.75 UVPM k = 0.75 Figure 5.24. Cmea vs. Time for a NACA 0010 airfoil plunging at reduced frequencies of k = 0.25 and 0.75 and an amplitude of ¯h = 0.025 67 0 2 4 6 8 10 12 14 16 18 20 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Cl t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.25. Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 5.2.3.3 Combined Pitching and Plunging Because solutions for an airfoil undergoing combined pitching and plunging motion represent a superposition of solutions for pure pitching and pure plunging, which have already been examined, this section will use a subset of the previously examined test cases to demonstrate that the panel code properly models combined pitching and plunging. It is expected that the same observations on the effects of amplitude and reduced frequency for the pure pitching and pure plunging cases will hold for the combined pitching and plunging. Figures 5.25 through 5.27 demonstrate the relative agreement between the panel code lift and moment solutions and the Theodorsen solution. Panel code solutions in Figures 5.25 through 5.27 are for a NACA 0010 airfoil pitching and plunging about the quarterchord at a reduced frequency of k = 0.25 with amplitudes of ¯α = 1, 2, and 4 deg, and ¯h = 0.025. Figures 5.28 through 5.30 demonstrate the relative agreement between the panel code lift and moment solutions and the Theodrsen solution. Panel code solutions in Figures 5.28 through 5.30 are for a NACA 0010 airfoil pitching and plunging about the quarterchord at 68 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 Cmle t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.26. Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmea t Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 2 h = 0.025 UVPM α = 2 h = 0.025 Theodorsen α = 4 h = 0.025 UVPM α = 4 h = 0.025 Figure 5.27. Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 1, 2, and 4 deg, and ¯h = 0.025. 69 0 2 4 6 8 10 12 14 16 18 20 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 Cl t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.28. Cl vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. a reduced frequency of k = 0.25 with a amplitudes of ¯α = 1 deg and ¯h = 0.010, 0.025, and 0.050. 5.2.3.4 Discussion As observed in the pure pitching and pure plunging examples, the panel code solution showed small variations in both phase and amplitude as compared to the Theodorsen solution. It was also shown that these small variations were dependent only on airfoil thickness and reduced frequency. Since the variations do not do not show a dependence on motion amplitude, the variation between the two solutions may be attributed to differences inherent in the wake models. As described earlier, the Theodorsen solution models the shed bound vorticity as a continuous vortex sheet of variable strength, released from the undisturbed airfoil trailing edge. The vortex sheet then convects with the timeaveraged freestream flow, essentially confining the vortex sheet to the x1 axis. The panel code models the shed circulation as a set of discrete vortex elements, released from the airfoil trailing edge location at each time step. The discrete 70 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cmle t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.29. Cmle vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmea t Theodorsen α = 1 h = 0.010 UVPM α = 1 h = 0.010 Theodorsen α = 1 h = 0.025 UVPM α = 1 h = 0.025 Theodorsen α = 1 h = 0.050 UVPM α = 1 h = 0.050 Figure 5.30. Cmea vs. Time for a NACA 0010 airfoil pitching and plunging about x = c/4 at k = 0.25, ¯α = 2 deg, and ¯h = 0.010, 0.025, and 0.050. 71 vortex elements are allowed to convect with the instantaneous local flowfield, and thus the wake is allowed to deform in time. By convecting the wake at the local velocity and allowing the wake to deform, the panel code wake model induces a different influence on the airfoil than the Theodorsen wake model, which would be dependent on reduced frequency. 72 Figure 5.31. Stationary plate of infinitesimal thickness with periodic transverse gust 5.3 Sears Periodic Gust The Sears periodic gust problem examines timedependent lift and moment generated on a stationary airfoil under the influence of a timeaveraged uniform freestream flow with sinusoidal transverse velocity perturbations, or transverse gusts. 5.3.1 Description Consider a stationary airfoil immersed in a freestream flow where the timeaveraged flow is aligned with the airfoil chord. If the airfoil is symmetric, it will not generate lift. However, if a transverse velocity perturbation is introduced into the freestream, the velocity perturbation will induce a local timedependent angle of attack change on the airfoil that will induce unsteady lift. The Sears gust problem investigates the influence of periodic transverse velocity perturbations on unsteady lift and moment generated on a flat plate. 73 5.3.2 Solution By assuming the transverse gust is not influenced by the presence of the airfoil, i.e. the “frozen gust” approximation, the downwash induced by the gust on the airfoil as a function of time can be writen as w = ¯weiω(t− x U ) = ¯weik(s−x∗) (5.13) Using a wake model similar to that of Wagner and Theodorsen, the relation for wake induced downwash on the airfoil as a function of reduced gust frequency is given by the Sears function. S (k) = 2 πk [H2 0 (k) − iH2 1 (k)] (5.14) Modifying the noflow boundary condition to include both gust and wake induced downwash, Sears lift and moment solutions for a flat plate under the influence of a sinusoidal transverse gust become L = 2πρU ¯ wbeiωt S (k) (5.15) My = L 1 2 + a b (5.16) 5.3.3 Comparison The panel code can implement two separate methods to model a periodic transverse gust. The first method, refered to later as the modified panel code, accounts for the influence of the periodic gust directly by modifying the implementation of the noflow boundary condition on the airfoil surface. The modified boundary condition includes the influence of the velocity pertubation by replacing the constant freestream velocity term with a time and position dependent function. This time and position dependent function must then be included in every other calculation which depends on the freestream velocity, such as the computation of unsteady surface pressures and the convection routines used to convect wake elements with the local flowfield. The application of this method to other problems, such as forced responce, is limited because gust influence on the airfoil is directly modeled as a function of 74 time and location in the flowfield, and as such does not allow for gust deformation due to body or wake influences. The second method employs the gust model descibed in Section 4.2. To use this discrete gust model, the continuous periodic gust is discretized into a set of gust sheets which propagate across the airfoil at the local flow velocity. This method does not require modifications to the original noflow boundary condition since the gust sheet is composed of constant strength vortex elements whose influence was included in the original noflow boundary condition. Since the gust sheets convect with the local flowfield, this method allows for gust deformation due to body and wake influences. The comparison of this method to the Sears solution is only limited by the discretization of the continuous periodic gust into a corresponding gust sheet representation. As such, care must be taken in choosing the method of discretization, since different representations of the same continuous gust will result in different lift and moment solutions. 5.3.3.1 Modified NoFlow Boundary Condition Figures 5.32 through 5.37 demonstrate the effect of airfoil thickness on panel code lift and moment solutions as compared to the Sears lift and moment solutions for reduced frequencies of k = 0.25, 1.0, and 4.0. Panel code solutions in Figures 5.32 through 5.37 were computed for NACA 0006, 0010, 0012, and 00014 airfoils under the influence of a continuous sinusoidal gust having an amplitude of ¯ w = 0.01 using the modified noflow boundary condition. Airfoil thickness exhibits a small influence on the phase between the modified panel code and Sears lift solutions as well as the amplitude ratio of the two solutions. The difference in amplitude of the lift solution computed by the panel code is attributed to the effects of airfoil thickness, because the solutions approach, but do not reach the Sears solution as airfoil thickness decreases. It is interesting to note that the amplitude ratio, computed as the maximum panel code lift coefficent divided by the maximum Sears lift coefficent, remains constant at roughly 0.7 for a NACA 0010 airfoil regardless of reduced frequency. The 75 0 2 4 6 8 10 12 14 16 18 20 −0.1 −0.05 0 0.05 0.1 0.15 Cl t Sears k = 0.25 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.32. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 influence of reduced frequency on the phase of the solutions is more dramatic. The phase of the modified panel code solution lags the Sears solution at small reduced frequencies, and shifts such that it leads at higher reduced frequencies. Because the amplitude ratio does not appear to be influenced by reduced frequency, it is assumed that differences in wake models, as discussed in Section 5.2.3.4, are responsible for the phase shift with reduced frequency. 5.3.3.2 Vortex Sheet Gust Model Figures 5.38 and 5.39 compares lift and moment coefficents calculated by the panel code utilizing the freestream gust model to a corresponding Sears solution. Panel code solutios in Figures 5.38 and 5.39 were computed for a NACA 0010 airfoil oriented at α0 = 0.0 deg to the timeaveraged freestream. The gust, having a reduced frequency of k = 1.0 and amplitude of ¯ w = 0.01, was modeled using a set of six gust sheets per gust period for five and a half periods upstream of the airfoil. The strength of each gust sheet is based on the velocity pertubation at the gust sheet’s initial x1 location in the flowfield. Figure 5.40 shows 76 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.01 0 0.01 0.02 0.03 Cmle t Sears k = 0.25 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.33. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 0.25 and a gust amplitude of ¯ w = 0.01 0 2 4 6 8 10 12 14 16 18 20 −0.04 −0.02 0 0.02 0.04 0.06 Cl t Sears k = 1.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.34. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 77 0 2 4 6 8 10 12 14 16 18 20 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 Cmle t Sears k = 1.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.35. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 1.0 and a gust amplitude of ¯ w = 0.01 0 1 2 3 4 5 6 7 8 9 10 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 Cl t Sears k = 4.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.36. Cl vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 78 0 1 2 3 4 5 6 7 8 9 10 −0.01 −0.005 0 0.005 0.01 0.015 Cmle t Sears k = 4.00 w = 0.01 UVPM NACA 0006 UVPM NACA 0010 UVPM NACA 0012 UVPM NACA 0014 Figure 5.37. Cmle vs. Time for the Sears solution compared to the alternate panel code solution for NACA 0006, 0010, 0012, and 0014 airfoils under the influence of a sinusoidal gust with a reduced frequency of k = 4.0 and a gust amplitude of ¯ w = 0.01 circulation strength per unit length about each gust sheet relative to the initial x1 location of the sheet. One drawback to the use of the discrete gust model is that the gust sheets do not produce a sinusoidal velocity pertubation. The induced pertubation due to the gust sheets closely resembles a set of sharp edge gusts, as shown in Figure 5.41. Figure 5.41 combines a visualization of the location of the airfoil, wake, and gust sheets in the top panel, the instantaneous pressure coefficients along the airfoil in the lower 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


