THE INFLUENCE OF CHORDWISE FLEXIBILITY ON THE FLOW STRUCTURE AND
STREAMWISE FORCE OF A SINUSOIDALLY PITCHING AIRFOIL
By
David Arthur Olson

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering – Doctor of Philosophy
2017

ABSTRACT
THE INFLUENCE OF CHORDWISE FLEXIBILITY ON THE FLOW STRUCTURE AND
STREAMWISE FORCE OF A SINUSOIDALLY PITCHING AIRFOIL
By
David Arthur Olson
Many natural flyers and swimmers need to exploit unsteady mechanisms in order to generate
sufficient aerodynamic forces for sustained flight and propulsion. This is, in part, due to the low
speed and length scales at which they typically operate. In this low Reynolds number regime,
there are many unanswered questions on how existing aerodynamic theory for both steady and
unsteady flows can be applied. Additionally, most of these natural flyers and swimmers have
deformable wing/fin structures, three dimensional wing planforms, and exhibit complex kinematics
during motion. While some biologically-inspired studies seek to replicate these complex structures
and kinematics in the laboratory or in numerical simulations, it becomes difficult to draw explicit
connections to the existing knowledge base of classical unsteady aerodynamic theory due to the
complexity of the problems. In this experimental study, wing kinematics, structure, and planform
are greatly simplified to investigate the effect of chordwise flexibility on the streamwise force
(thrust) and wake behavior of a sinusoidally pitching airfoil.
The study of flexibility in the literature has typically utilized flat plates with varying thicknesses
or lengths to change their chordwise flexibility. This choice introduces additional complexities
when comparing to the wealth of knowledge originally developed on streamlined aerodynamic
shapes. The current study capitalizes on the recent developments in 3D printer technology to
create accurate shapes out of materials with varying degrees of flexibility by creating two standard
NACA 0009 airfoils: one rigid and one flexible.
Each of the two airfoils are sinusoidally pitched about the quarter chord over a range of oscillation amplitudes and frequencies while monitoring the deformation of the airfoil. The oscillation
amplitude is selected to be small enough such that leading edge vortices do not form, and the vortical structures in the wake are formed from the trailing edge. Two-component Molecular Tagging

Velocimetry (MTV) is employed to measure the vortical flowfield over the first chord length behind
the airfoil. A control volume method is used to estimate the mean thrust of the airfoil based on the
mean and fluctuating velocity profiles from the MTV results.
The mean thrust results show chordwise flexibility increases the thrust produced by the airfoil
over the range of motion parameters and the flexibility considered in this study. The flexible airfoil
is also seen to experience the drag-to-thrust crossover at a lower oscillation frequency than its rigid
counterpart. The relative change in thrust due to flexibility decreases with increasing amplitude.
The increase in thrust can, however, be captured as an amplitude effect when the Strouhal number
based on the actual trailing edge displacement, Stte , is used for scaling. Scaling based strictly on
the prescribed motion, typically employed in the literature, is not sufficient for the data to collapse.
Motion trajectories which produced a classical von KĂĄrmĂĄn vortex street or a reverse von
KĂĄrmĂĄn vortex street (depending on the arrangement of the vortices), are considered for further
study. The vortices in the wake are characterized in terms of their strength, size, and spacing using
phase-averaged MTV results. The circulation of the vortices are shown to collapse for both rigid
and flexible airfoils when plotted against Stte . The actual trailing edge displacement is used as a
length scale to normalize the transverse and streamwise spacing, and the vortex core size. These
measurements also now collapse when plotted against Stte across oscillation amplitude for both the
rigid and flexible airfoils.

Copyright by
DAVID ARTHUR OLSON
2017

This thesis is dedicated to science.

v

ACKNOWLEDGEMENTS

Few things in life are accomplished individually. There is nearly always someone or something
that has significantly contributed or enabled one to succeed. I will make an attempt to express my
gratitude to those individuals and entities. First and foremost, I’d like to humbly thank my wife,
Olivia, for all of the love and support she gives me. Secondly, I’d like to thank all of my friends
and family (near and far), for all their support over the many years.
I’ve had the opportunity to have two advisors during my entire time here at MSU. Dr. Ahmed
Naguib has been an invaluable co-advisor; his ability to translate graduate student ramblings into
coherent messages is remarkable, especially when I almost got the point across, but it needed a
little tweaking to actually make sense. His work ethic and still-in-the-trenches style of advising is
aspirational. Lastly, I need to thank one of the best advisors, Dr. Manoochehr Koochesfahani, I
have had the opportunity to work with, meet, or even hear about. Sometimes, I would feel left out
when graduate students get together and complain about their advisors, but then I would remember
that I had it pretty good. I sincerely thank you for all of the guidance and support you have given
me over the years.
Visitors to TMUAL commonly ask how is it that Dr. Patrick Hammer and I get along, with him
being a computationalist and me an experimentalist working on the same projects. My response is
always: Patrick has thick skin; and for that, I am grateful. Patrick’s high fidelity simulations have
been an invaluable resource, especially the numerous cases he ran specifically to compare with my
experiments. For all of that, and more, Thanks! I’d like to thank all of the rest of the members (and
escapes) of TMUAL over the years. I’ll name a few: Dr. Rohit Nehe (for many years of friendship
and camaraderie), Dr. Bruno Monnier (who designed and build the 3-DOF motion apparatus while
I got to be the first one to use it), and Dr. Shahram Pouya (whose kind demeanor is only surpassed
by his knowledge of TMUAL).
Lastly, this work was sponsored by AFOSR award numbers FA9550-10-1-0342 and FA955015-1-0224.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

CHAPTER 1 INTRODUCTION
1.1 Motivation . . . . . . . .
1.2 Background . . . . . . .
1.3 Flexibility . . . . . . . .
1.4 Current Study . . . . . .
1.5 Outline . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

1
1
2
4
6
8

CHAPTER 2 EXPERIMENTAL METHODS . . . . . . . . . . . . .
2.1 Flow Facility . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Description . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Flow Boundary Conditions . . . . . . . . . . . . . . .
2.2 Airfoil Motion Trajectory . . . . . . . . . . . . . . . . . . . .
2.2.1 Equipment Description . . . . . . . . . . . . . . . . .
2.2.2 Prescribed Motion . . . . . . . . . . . . . . . . . . .
2.2.3 Measured Trajectory . . . . . . . . . . . . . . . . . .
2.3 Airfoil Properties . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 3D Printer Description . . . . . . . . . . . . . . . . .
2.3.2 Material Testing . . . . . . . . . . . . . . . . . . . . .
2.3.3 Design and Shape Validation . . . . . . . . . . . . . .
2.3.4 Estimated Flexibility . . . . . . . . . . . . . . . . . .
2.3.5 Deformation Tracking . . . . . . . . . . . . . . . . . .
2.4 Molecular Tagging Velocimetry . . . . . . . . . . . . . . . .
2.4.1 Background . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Current Implementation . . . . . . . . . . . . . . . .
2.4.3 MTV Imaging System . . . . . . . . . . . . . . . . .
2.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Phase Averaged Data . . . . . . . . . . . . . . . . . .
2.5.3 Sample Measurements . . . . . . . . . . . . . . . . .
2.5.4 Control volume method for estimating streamwise force
2.5.5 Wake Vortex Properties . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

10
10
10
11
14
14
15
17
19
20
21
26
30
32
36
36
37
38
40
40
44
45
50
51

CHAPTER 3 STREAMWISE FORCE
3.1 Rigid Airfoil Results . . . . .
3.1.1 Velocity Profiles . . .
3.1.2 Streamwise forces . . .
3.2 Effect of Flexibility . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

54
54
54
60
64

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

vii

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

3.3

3.2.1 Velocity Profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2.2 Streamwise force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

CHAPTER 4 WAKE VORTEX PROPERTIES . . . . . . . . . . . . . . . . . . . .
4.1 Validation of Rigid Airfoil Results . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Comparison with 2D Simulations of the NACA 0009 airfoil . . . . . .
4.1.2 Effect of NACA 0009 on vortex properties . . . . . . . . . . . . . . .
4.1.3 Comparison with Experiments and 2D Simulations of the NACA 0012
4.2 Effect of flexibility on vorticity field . . . . . . . . . . . . . . . . . . . . . .
4.3 Wake vortex properties of rigid and flexible airfoils . . . . . . . . . . . . . .
4.3.1 Peak Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Vortex core size . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 Convective Speed and Wavelength . . . . . . . . . . . . . . . . . . .
4.3.5 Transverse Spacing . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
CHAPTER 5

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

74
74
74
79
82
84
86
86
90
92
96
100
105

CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
APPENDIX A
TWO CAMERA CALIBRATION . . . . . . . . . . . . . . . . .
APPENDIX B
UNCERTAINTY ESTIMATES . . . . . . . . . . . . . . . . . . .
APPENDIX C
EFFECT OF VORTICITY CUTOFF LEVEL . . . . . . . . . . .
APPENDIX D
STREAMWISE FORCES . . . . . . . . . . . . . . . . . . . . . .
APPENDIX E
VORTICITY FIELDS FOR RIGID AND FLEXIBLE AIRFOILS .
APPENDIX F
POSITIVE AND NEGATIVE VORTEX PROPERTIES . . . . . .
APPENDIX G
SPATIAL RESOLUTION AFFECT ON VORTEX PROPERTIES

.
.
.
.
.
.
.
.

108
109
114
120
123
129
133
136

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

viii

LIST OF TABLES

Table 2.1: Non-dimensional bending stiffness parameters of rigid and flexible airfoils . . . . 31
Table 3.1: Deviations of rigid airfoil results and numerical simulations of NACA 0009 . . . 55
Table 4.1: Vortex parameter nondimensionalizations and scaling variable . . . . . . . . . . 105
Table C.1: Effect of vorticity cutoff level on rc and Γ . . . . . . . . . . . . . . . . . . . . . 122
Table G.1: Change in the vortex properties due to reduced spatial resolution . . . . . . . . . 137

ix

LIST OF FIGURES

Figure 1.1: Schematic of a reverse von KĂĄrmĂĄn vortex street . . . . . . . . . . . . . . . . .

3

Figure 2.1: Scaled drawing of the entire flow facility . . . . . . . . . . . . . . . . . . . . . 11
Figure 2.2: Airfoil and end walls assembly, as installed in the water tunnel . . . . . . . . . 13
Figure 2.3: The three degree of freedom apparatus above the test section . . . . . . . . . . 15
Figure 2.4: Electronics chassis and local network topology . . . . . . . . . . . . . . . . . . 16
Figure 2.5: The average maximum deviation from the best fit sine wave of the motion
trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Figure 2.6: The average RMS deviation from the best fit sine wave of the motion trajectory . 19
Figure 2.7: Spanwise bending measured at the tip of the airfoil . . . . . . . . . . . . . . . . 20
Figure 2.8: Test specimen of 3D printer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Figure 2.9: Photochemical properties of 3D printer materials after exposure to MTV triplex

23

Figure 2.10: Trailing edge deflection as a function of the submersion time . . . . . . . . . . 24
Figure 2.11: Change in the trailing edge deflection from the first sFOV of the day to the
second sFOV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Figure 2.12: The change in the trailing edge deflection across all FOVs of the flexible airfoil . 27
Figure 2.13: Cross-sectional views of the 3D printed airfoils . . . . . . . . . . . . . . . . . . 27
Figure 2.14: Measured thickness distribution of the rigid and flexible airfoils . . . . . . . . . 29
Figure 2.15: Measured airfoil profile of the rigid and flexible airfoils . . . . . . . . . . . . . 29
Figure 2.16: Shape monitoring camera setup and sample image . . . . . . . . . . . . . . . . 33
Figure 2.17: Trailing edge displacement for all motion trajectories of the flexible airfoil . . . 34
Figure 2.18: Relative change in trailing edge displacement . . . . . . . . . . . . . . . . . . . 35
Figure 2.19: Example 2c-MTV image pair and resulting 2-component velocity vector field . 36
x

Figure 2.20: Top view of optical beam path of excimer laser . . . . . . . . . . . . . . . . . . 37
Figure 2.21: Photograph of dual camera MTV imaging system . . . . . . . . . . . . . . . . 39
Figure 2.22: Fields of view for the rigid airfoil . . . . . . . . . . . . . . . . . . . . . . . . . 41
Figure 2.23: Fields of view for the flexible airfoils . . . . . . . . . . . . . . . . . . . . . . . 42
Figure 2.24: Summary of all motion parameters . . . . . . . . . . . . . . . . . . . . . . . . 43
Figure 2.25: Representative mean and rms streamwise velocity field . . . . . . . . . . . . . 47
Figure 2.26: Representative mean and rms transverse velocity field . . . . . . . . . . . . . . 48
Figure 2.27: Representative vorticity field and its uncertainty . . . . . . . . . . . . . . . . . 49
Figure 2.28: Schematic of typical vortex pattern in wake of oscillating airfoil . . . . . . . . . 51
Figure 3.1: Velocity profiles of the rigid NACA 0009 airfoil compared to the numerical
simulations of Hammer (2017) . . . . . . . . . . . . . . . . . . . . . . . . . . 56
Figure 3.2: Velocity profiles of the rigid airfoil compared to the experiments of Bohl &
Koochesfahani (2009) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Figure 3.3: Mean Ct of the rigid airfoil and numerical simulations at α0 = 2◦ . . . . . . . . 61
Figure 3.4: Mean thrust coefficient (Ct ) plotted against reduced frequency (k) for various
different studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Figure 3.5: Mean thrust coefficient (Ct ) plotted against Strouhal number (St) for various
different studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Figure 3.6: Selected mean velocity profiles at α0 = 2◦ . . . . . . . . . . . . . . . . . . . . 66
Figure 3.7: Selected fluctuating streamwise velocity profiles at α0 = 2◦ . . . . . . . . . . . 67
Figure 3.8: Selected mean streamwise velocity profiles at α0 = 4◦ and α0 = 6◦ . . . . . . . 69
Figure 3.9: Mean thrust plotted against reduced frequency for rigid and flexible airfoils . . . 70
Figure 3.10: Mean thrust plotted against Strouhal number for rigid and flexible airfoils . . . . 71
Figure 3.11: Mean thrust coefficient plotted against Strouhal number based on the actual
trailing edge displacement for rigid and flexible airfoils . . . . . . . . . . . . . 72

xi

Figure 4.1: Comparison of experiment and simulations of phase-averaged vorticity fields . . 76
Figure 4.2: Vortex Properties of Rigid Airfoil and Simulations of NACA 0009 . . . . . . . 78
Figure 4.3: Vortex Properties of the rigid airfoil compared to numerical simulations of
the NACA 0009 and NACA 0012 airfoils . . . . . . . . . . . . . . . . . . . . . 80
Figure 4.4: Comparisons of vortex properties of the rigid airfoil and other experiments
and simulations of a NACA 0012 airfoil . . . . . . . . . . . . . . . . . . . . . 83
Figure 4.5: Phase-averaged vorticity field of rigid (left) and flexible (right) airfoils . . . . . 85
Figure 4.6: Non-dimensional peak vorticity for different oscillation amplitudes of the
rigid and flexible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Figure 4.7: Scaling of the peak vorticity by h p /U∞ . . . . . . . . . . . . . . . . . . . . . . 88
Figure 4.8: Peak Vorticity scaling with trailing edge acceleration . . . . . . . . . . . . . . 89
Figure 4.9: Nondimensional circulation plotted against k and St, Stte . . . . . . . . . . . . 91
Figure 4.10: Circulation normalized following the scaling of Buchholz et al. (2011) . . . . . 93
Figure 4.11: Vortex core size plotted against k and St . . . . . . . . . . . . . . . . . . . . . 94
Figure 4.12: Core size of rigid and flexible airfoils plotted against Stte . . . . . . . . . . . . 95
Figure 4.13: Vortex convective speed of the rigid and flexible airfoils . . . . . . . . . . . . . 97
Figure 4.14: Vortex convective speed of the rigid and flexible airfoils compared to selected
literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Figure 4.15: Wavelength scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Figure 4.16: Wavelength normalized by the trailing edge displacement. . . . . . . . . . . . . 100
Figure 4.17: Transverse spacing normalized by trailing edge displacement of the rigid and
flexible airfoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Figure 4.18: Transverse spacing normalized by trailing edge displacement compared with
selected literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Figure 4.19: Transverse spacing normalized by wavelength for the rigid and flexible airfoils . 103
Figure 4.20: Transverse spacing normalized by wavelength compared with selected literature 104
xii

Figure A.1: MTV Imaging System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
Figure A.2: Targets of each camera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Figure A.3: Uncorrected displacement field between the two cameras . . . . . . . . . . . . 112
Figure A.4: Corrected displacement field using affine transform . . . . . . . . . . . . . . . 113
Figure B.1: Independent measurements of Ct plotted against reduced frequency (k) . . . . . 119
Figure C.1: Effect of vorticity cutoff on Γ and rc . . . . . . . . . . . . . . . . . . . . . . . 121
Figure D.1: Mean thrust coefficient plotted against k, full range . . . . . . . . . . . . . . . 123
Figure D.2: Mean thrust coefficient plotted against k, selected range . . . . . . . . . . . . . 124
Figure D.3: Mean thrust plotted against St, full range . . . . . . . . . . . . . . . . . . . . . 125
Figure D.4: Mean thrust coefficient plotted against St, selected range . . . . . . . . . . . . . 126
Figure D.5: Mean thrust coefficient plotted against Stte , full range . . . . . . . . . . . . . . 127
Figure D.6: Mean thrust coefficient plotted against Stte , selected range . . . . . . . . . . . . 128
Figure E.1: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 2◦ . . . . . . 130
Figure E.2: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 4◦ . . . . . . 131
Figure E.3: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 6◦ . . . . . . 132
Figure F.1: Peak vorticity of both positive and negative sign vortices . . . . . . . . . . . . . 134
Figure F.2: Circulation of both positive and negative sign vortices . . . . . . . . . . . . . . 135

xiii

CHAPTER 1
INTRODUCTION

1.1

Motivation

There are numerous engineering applications for unmanned aerial vehicles such as remote hazardous material identification, search and rescue, and reconnaissance. Many of these applications
require vehicles that are small and lightweight, maneuverable, and operate in complex environments. Collectively, these vehicles are typically referred to as Micro Air Vehicles (MAVs), and
engineers have turned to nature for inspiration as these vehicles are the same physical size of small
birds and/or fish (Spedding & Lissaman, 1998). At the low speeds and length scales in which
these vehicles operate, viscous effects become increasingly important (Lissaman, 1983); necessitating many of the natural flyers and swimmers, and the biomimetic MAVs, to exploit unsteady
mechanisms in order to generate sufficient aerodynamic forces for sustained flight and propulsion.
One challenge of investigating the performance and aerodynamics of complex biomimetic
wings is generalizing the results to enhance our understanding of the underlying mechanisms of
lift and thrust generation. The wings (or fins) of the natural flyers (or swimmers) are structurally
intricate and undergo complex kinematics; one approach to studying the problem, from a direct
biological perspective, utilizes replicas or the actual insect, bird, or bat wing and are experimentally
or numerically studied with either realistic or simplified kinematics (Spedding & Lissaman, 1998;
Shyy et al., 2010). While these types of studies have yielded tremendous insight into the flight
and aerodynamics of the natural flyers, it is difficult to generalize the results to understand the
fundamental unsteady aerodynamic and aeroelastic effects. Alternatively, another approach is to
start with more traditional airfoils under 2D boundary conditions and simplified kinematics in order
to parametrically investigate the lift and thrust generation. Systematically increasing the complexity
of the airfoils, kinematics, boundary conditions, or a combination thereof, is seen as a path for
improving our understanding of the fundamental unsteady aerodynamics.

1

1.2

Background

The initial theoretical explanations of the generation of thrust of flapping flight are credited to
Knoller (1909) and Betz (1912), with Katzmayr (1922) being the first to experimentally observe the
effect. The classical unsteady aerodynamic theory was developed typically under the assumptions
of inviscid flow, low amplitude oscillations, thin and high aspect ratio (AR) wings, as motivated
by the modeling of the forces associated with flutter (Theodorsen, 1934; von KĂĄrmĂĄn & Burgers,
1935; Garrick, 1937; von KĂĄrmĂĄn & Sears, 1938). The complex kinematics of flapping wings
were reduced to several canonical motion types: pitching, plunging or heaving, and combined
pitch-plunge motion. It is known from these classical unsteady theories that the generation of
thrust can be accomplished by simple oscillatory pitching motion. While the current study focuses
on pure harmonic pitching motion to simplify the kinematics to the minimum necessary to generate
thrust, some background on harmonic plunging motion is required as both motion kinematics are
studied throughout the relevant literature.
The generation of thrust from an oscillating airfoil was connected to an array of alternating-sign
vortices shed from the airfoil by von KĂĄrmĂĄn & Burgers (1935) and Lighthill (1969). Koochesfahani
(1989) showed experimentally for a pitching airfoil that the shed vorticity form these coherent
structures as a von KĂĄrmĂĄn vortex street, indicative of a momentum deficit, and at a high enough
frequency, a so-called reverse von KĂĄrmĂĄn street, indicative of a momentum excess or jet-like
velocity profile. Ramamurti & Sandberg (2001) computed similar results for the pitching airfoil,
while similar patterns were seen for a harmonically plunging airfoil experimentally by Freymuth
(1988) and numerically by Platzer et al. (1993), who also later saw it experimentally in Jones et al.
(1998) and Lai & Platzer (1999).
An array of alternating-sign vortices is shown schematically in Figure 1.1 where a pitching airfoil
about its quarter chord (with a trailing edge displacement of h p ) generates a reverse von KĂĄrmĂĄn
vortex street. The vortex pattern is quantified by its transverse spacing, b, and its wavelength,
a. The transverse spacing (b) is the displacement from the center of the negative vortex (yc,n ) to
the positive vortex (yc,p ). The induced velocity from the vortices, which may be visualized by
2

y

a

U∞

x
yc,p
hp
yc,n

+
–

+
–

b

Figure 1.1: Schematic of a reverse von KĂĄrmĂĄn vortex street. An airfoil pitching about the quarter
chord is depicted with the transverse spacing (b), wavelength (a), and the trailing edge displacement
(h p ) indicated. The sign of the circulation of the vortex is indicated by the arrow and by the (+) or
(-) sign.

the their sign in circulation (Γ) rotation direction shown in Figure 1.1, creates a velocity excess
in the streamwise direction when b is positive. When the vortices are arranged in the opposite
orientation (with the positive sign vortices below the negative sign vortices) the induced velocities
of the vortices would create a velocity deficit in the streamwise direction.
Bohl & Koochesfahani (2009) found the switch from drag to thrust does not correspond with
the frequency at which the vortex pattern changes, but they are close to one another; the switch
in the vortex pattern occurs at a lower frequency than the drag-to-thrust crossover frequency.
Bohl & Koochesfahani (2009) also studied the details of the vortex properties of a harmonically
pitching NACA 0012 airfoil at a chord Reynolds number of Rec = 1.25 × 104 for an amplitude of
α0 = 2◦ , and α0 = 4◦ is included in Bohl (2002). The peak vorticity (ω pk ) and circulation (Γ)
were used to quantify the strength, and the radius of gyration rc , was used to quantify the size of
the vortices. The wavelength (a) and transverse spacing (b), as depicted in Figure 1.1, were also
quantified. Computationally, Hammer (2016) was able to closely replicate the experimental results
of Bohl & Koochesfahani (2009) with high fidelity 2D simulations. Hammer (2016) investigated
several additional pitching amplitudes, and Reynolds numbers while investigating the effect of these
parameters on the wake vortex properties.

3

1.3

Flexibility

Many natural flyers and swimmers have extremely complex three dimensional deformable wing
structures and undergo equally complex kinematics. Combes (2003) found that spanwise flexural
stiffness was typically 1-2 orders of magnitude higher than chordwise flexural stiffness for many
different insect species. A starting point for considering flexibility is to isolate the chordwise
direction as it appears to be the more flexible than the spanwise direction in many insects in nature.
An inviscid model developed by Katz & Weihs (1978, 1979) was used to investigate passive
chordwise flexibility for harmonically oscillating wings undergoing large amplitude motions along
curved paths. The oscillating airfoil in Katz & Weihs (1978) is taken to be moving at high Reynolds
number in an incompressible fluid such that potential theory developed earlier (Wu, 1961; Lighthill,
1975; James, 1973; Chopra, 1976) could be used. The airfoil motion along a curved path is such
that while the transverse oscillation amplitude is large, the angle of attack (Îą) variation is low. The
interaction between the chordwise flexible structure and the flow field was achieved by an elastic
model iteratively interacting with the pressure field described by the potential flow model. The
model utilizes the Kutta condition for shedding a discrete into the wake for the wake behavior;
an assumption that is now known to be invalid for unsteady low amplitude oscillations (Bohl &
Koochesfahani, 2009). Overall, Katz & Weihs (1978) observe that flexibility tends to improve
efficiency but at a sacrifice to thrust production. In their subsequent paper, Katz & Weihs (1979)
incorporated lifting line theory to account for the downwash of a 3D airfoil which resulted in
slightly modified, but otherwise the same conclusions as their earlier work.
More recently, Moore (2014), used a torsional flexibility model where the leading edge is
modeled as torsional spring, allowing the otherwise rigid airfoil to pitch about its leading edge.
Moore (2014) extends the model of Wu (1961) using a torque balance between the fluid flow and
the torsional spring that determines the pitch angle. Moore (2014) was able to capture the observed
“ideal flexibilities” seen previously in different experiments (Heathcote & Gursul, 2007; Alben,
2008; Spagnolie et al., 2010; Dewey et al., 2013; Monnier et al., 2015). Moore (2014) shows in
the limit of a massless compliant airfoil, thrust production is hindered compared to the rigid airfoil
4

but moderate flexibilities can increase thrust production. While this is not a new observation, the
simplicity of the model enabling these features to be seen for a passive pitching (or flexing) airfoil
is useful for further investigating these flows.
Moore (2015) extended his previous model by considering the wing as a bending beam with
chordwise dependent rigidity. They investigated the optimal thrust production from airfoils modeled
with a torsional spring at the leading edge, together with linear, quadratic, or cubic distributions of
chordwise rigidity. All of the ideal conditions for stiffness distribution resulted in the airfoil closely
resembling a rigid wing pitching about its leading edge. This is reminiscent of the “feathering
effect” as described by Katz & Weihs (1978) and the the idea mentioned by Zhu (2007), that while
complete feathering will produce almost no thrust, partial feathering would lead to a small amount
of thrust reduction, but increased efficiency. Moore (2014, 2015) does not address efficiency in his
models, and so, the efficiency of the “optimal” chordwise flexible airfoil remains unknown based
on their model.
There is little consensus in the literature regarding the mechanism that drives the efficiency
of passive chordwise flexible airfoils undergoing flapping motion (pitching, plunging or combined
motion). Two competing theories are: (1) the motion oscillation frequency is tuned to the natural
resonance frequency of the airfoil, and (2) the bending of the airfoil is timed such that its release of
stored elastic energy is beneficial to the flow structures produced.
There have been many studies considering chordwise flexible wings undergoing plunging
motion. Collectively, these studies (Heathcote & Gursul, 2007; Shyy et al., 2010; Kang et al., 2011;
Cleaver et al., 2014, 2016) have shown there exists an optimal amount of chordwise flexibility to
improve thrust production and propulsive efficiency. Cleaver et al. (2016) notes that the thrust
improvement is observed in conjunction with stronger trailing edge vortices with a larger transverse
spacing (which is known, based on vortex models such as Naguib et al. (2011), to correspond to an
increase in thrust). Ramananarivo et al. (2011), based on the plunging, self-propelled flexible wing
of Thiria & Godoy-Diana (2010), hypothesize that optimization of efficiency in unsteady flexible
flyers is not achieved through structural resonance but through an optimal timing of the release

5

of the stored elastic energy into creating flow structures that are beneficial for thrust. Thiria &
Godoy-Diana (2010) substantiate their hypothesis with a nonlinear oscillator model which shows
that energy transfer may be achieved at significantly lower frequency than structural resonance.
This is in contrast to the conclusions of Kang et al. (2011) which state that for a plunging airfoil, the
maximum thrust is achieved near resonance while maximum efficiency is reached at approximately
half of the natural frequency, which they found from numerical simulations and non-dimensional
modeling.
Some closure on the debate regarding chordwise flexibility and efficiency/thrust production was
presented by Dewey et al. (2013), albeit for a pitching panel compared to the previously mentioned
plunging airfoils. Dewey et al. (2013) observed that chordwise flexibility can increase or decrease
the thrust, depending on the amount of flexibility; consistent with other pitching airfoil studies
of Prempraneerach et al. (2003), Marais et al. (2012), and Monnier et al. (2013, 2015). Dewey
et al. (2013) also observed efficiency peaks below, at, and above the airfoil’s structural natural
frequency, depending on the flexibility. More importantly, Dewey et al. (2013) found that the peak
in efficiency for each airfoil typically occurred when the Strouhal number was between 0.25 < St <
0.35, which was observed by Triantafyllou (1993) in a survey of Strouhal numbers seen in various
fish species and through a stability analysis of the jet-like velocity profiles that flapping wings
generate. Furthermore, Dewey et al. (2013) observed that when this oscillation frequency (falling
within the 0.25 < St < 0.35 range) is also tuned to the first natural frequency, there exists a global
maximum in efficiency for the optimally-flexible airfoil.

1.4

Current Study

In this experimental study, wing kinematics, structure, and planform are greatly simplified to
investigate the effect of chordwise flexibility on the streamwise force (thrust) and vortical wake
structure of a sinusoidally pitching airfoil.
Most of the previous studies discussed thus far on both pitching and plunging chordwise flexible
airfoils focused primarily on thrust production and propulsive efficiency. Quantifying the vortical

6

flow structures shed into the wake have been given less attention in the literature, even for rigid
airfoils (Bohl & Koochesfahani, 2009; Hammer, 2016). The current study focuses on quantifying
the vortex parameters and spacing in the wake of a harmonically-pitching chordwise flexible airfoil
at a chord Reynolds number of nominally Rec = 1.2 × 104 , which is chosen to allow direct
comparisons to previous experiments (Bohl, 2002; Bohl & Koochesfahani, 2009) and simulations
(Hammer, 2016) on wake vortex properties. Consistent with these studies, the pitching amplitudes
considered at α0 = 2◦ , 4◦ , and 6◦ . The motivation to study wake vortex properties is two-fold. The
first is the potential use of low order discrete vortex models as an engineering tool for thrust/lift
predictions, as discussed in Naguib et al. (2011) and Monnier et al. (2015). The second is the
aforementioned general dearth of detailed studies into these flow fields. Detailed experimental
studies on such dynamic and intricate flowfields are important for validating high fidelity numerical
simulations.
Previous studies focusing on the vorticity in the wake of harmonically pitching flexible propulsors is limited. Cleaver et al. (2014) studied chordwise flexibility by considering a plunging
NACA 0012 airfoil with an attached trailing edge flap of different lengths and thicknesses. While
Cleaver et al. (2014) did not systematically study the wake vortex properties, they did observe, for
a few selected cases, increased circulation and increased transverse spacing of the trailing edge
vortices for the optimal flexible trailing edge flap compared to the rigid flap, albeit, for a plunging
airfoil.
Monnier et al. (2013, 2015) considered the vortex properties behind rigid and flexible propulsors,
but at a lower Reynolds number of Rec = 2.3 × 103 . The model of Monnier et al. (2013, 2015) was
not a standard airfoil shape, but a geometry similar to Heathcote & Gursul (2007) where the leading
edge was the shape of a standard airfoil but had a long trailing edge attached to it. This geometry is
a strong departure from standard airfoils where there is significant leading edge curvature that helps
delay leading edge separation. The Reynolds number was also low enough to affect the results
compared to Rec = 1.2 × 104 as discussed by Hammer (2016).
Prempraneerach et al. (2003) and Egan et al. (2016) created standard symmetric airfoils (by

7

creating urethane casts of different rigidities) to study chordwise flexibility. Their pitching amplitudes, α0 = 10◦ -30◦ , chord Reynolds numbers, Rec = 3 × 104 − 3 × 105 , are both considerably
higher than the current study.
The current research utilizes a slightly different approach than most of the other studies looking
at flexibility by studying one aerodynamic shape, specifically a NACA 0009, that is either rigid or
flexible in the chordwise direction. Uniformity in the span, with regards to flexibility and shape,
simplifies the boundary conditions and facilitates the direct comparisons to previous information
that is known for rigid airfoils; this is important for comparing to existing literature as few other
studies have investigated the vortex properties and their arrangement in the wake of the airfoil.
The oscillation amplitudes are kept low enough to prevent leading edge vortices from forming and
interacting with trailing edge vortices to simply the flowfield and subsequent discussion.

1.5

Outline

The remainder of this dissertation is broken into four chapters. Chapter 2 describes the
experimental methods used to measure and quantify the airfoil’s shape, shape deformation, and
trajectory. The technique and data processing applied to quantify the velocity field is discussed.
The method used to calculate the mean streamwise force (thrust) acting on the airfoil from the
mean and fluctuating velocity field is outlined. Further, the parameters used to quantify the vortex
properties and arrangement are described.
The primary results and discussion is broken into two chapters. The first, Chapter 3, begins with
a comparison of the mean velocity profiles with other studies under the same conditions. Chapter 3
then compares the mean velocity profiles of the rigid and flexible airfoils, and lastly, compares the
mean thrust of the airfoil based on a control volume method. The second results chapter, Chapter 4,
initially compares the vorticity field characteristics of the rigid airfoil in this study to those in
the literature. The influence of chordwise flexibility on the various vortex properties and spacing
parameters is then compared and discussed. A primary focus is on the non-dimensional parameters
that are useful for scaling the properties for overall collapse between the rigid and flexible airfoil.

8

Lastly, Chapter 5 concludes the dissertation with a brief summary of the important results and
outlines future areas of research on the topic.

9

CHAPTER 2
EXPERIMENTAL METHODS

2.1
2.1.1

Flow Facility
Description

Experiments are conducted in a water tunnel located in the Turbulent Mixing and Unsteady Aerodynamics Laboratory (TMUAL) at Michigan State University (MSU). The water tunnel is a 10,000
liter closed-return facility with a 61 cm × 61 cm × 244 cm test section (shown in Figure 2.1). The
flow management consists of two perforated plates in the initial diffuser followed by a 0.25" cell
diameter honeycomb and fine mesh screen at the entrance of the settling chamber. After the settling
chamber, a 6:1 contraction leads to the entrance of the test section where a 0.125" cell diameter
honeycomb and fine mesh screen are located.
The test section has full optical access in the visible spectrum on both sides, the bottom, and
a full height end-view window in the return plenum. Two quartz windows (41 cm × 84 cm each
which can be seen in the top of Figure 2.1) were installed on one side of the test section to extend the
transmission range of the test section into the deep UV (ultraviolet) and near IR (infrared) ranges.
The impeller of the water tunnel is powered by a Toshiba 20 hp motor and controlled with
a Toshiba VF-AS1 drive which ensures day-to-day repeatability of the same freestream velocity;
a necessity since the experiments were conducted over a six month time period. The maximum
freestream velocity (U∞ ) of the facility is 1 m/s but is set to 10 cm/s in all of the current experiments.
The water level was maintained at 63.5 cm in the test section to ensure all air pockets have been
removed from the sealed sections upstream and downstream of the test section. This helps reduce
some of the low frequency oscillations that are naturally present in these types of facilities. The
freestream turbulence intensity of the facility (FST I = urms /U∞ ) consists of nearly spatially
uniform low frequency (<0.2 Hz) oscillations in addition to broadband turbulence; the FST I of the

10

Flow

→

Figure 2.1: Scaled drawing of the entire flow facility. The facility is 8.6 m in length, 1.9 m in height,
and 2.5 m in width. The test section, translucent in the side profile view (top), has an internal cross
section of 61 cm × 61 cm and is 244 cm in length. The quartz windows are each 41 cm × 84 cm
installed on the back side of the test section, as shown by darker translucent rectangular sections in
top drawing.

facility, most recently measured in Olson (2011), is 1.8% from all contributions and 0.5% from
broadband turbulence.

2.1.2

Flow Boundary Conditions

The scope of the current work has defined the desired boundary conditions to be those that
approximate 2D flow conditions, i.e. large aspect ratio (AR = span/c). The span of the airfoil
is limited to 33.6 cm due to manufacturing limitations described in Section 2.3, which yields an
aspect ratio of AR = 2.8, since the chord was kept at 12 cm for consistency with previous studies

11

in the facility. With this AR, the current study will be effected by three dimensionality more than
the previous studies in this facility, notably Bohl & Koochesfahani (2009) and Bohl (2002) who
had an AR = 4. The impact of three dimensionality on the conclusions drawn from this study will
be discussed in Chapter 5.
A false wall on the bottom of the test section was designed, fabricated, and installed in the test
section to compensate for the span of the airfoil being less than the height of the test section. The
assembly of the false walls and airfoil is shown in Figure 2.2. The walls are 91.4 cm × 59.7 cm
(7.6c × 4.9c) in length and width, respectively, and are made from scratch resistant polycarbonate.
The airfoil’s leading edge is located in the center of the test section’s width and 2c downstream
from the leading edge of the plates. The airfoil is also centered vertically between the two false
walls with a gap distance of < 2 mm between the wall and the airfoil. The lower plate’s leading
edge is rounded and the entire plate is offset 27.6 cm from the bottom of the test section. A 7.6 cm
long flap is attached to the plate with a stainless steel piano hinge and was adjusted to an angle of
15◦ ± 1◦ on the top (airfoil side) of the false wall. The upper wall has a single radius on the leading
edge and is suspended approximately 15 cm from the top of the test section such that the water level
was at the middle of the plate. The center portion of the top plate has two removable sections that
allow the airfoil shaft assembly (a 3.81 cm diameter stainless steel cylinder) to pass through and
connect to the motors.
The airfoil is considered to have stationary wall boundary conditions located at Âąy=2.54c and
at Âą z=1.4c with a small rotating wall of diameter 0.31c centered around the pitch axis (0.25c) at
the upper wall.

12

z
U∞

y
x

Figure 2.2: Airfoil and end walls assembly, as installed in the water tunnel. The translucent parts
are made from chemical resistant polycarbonate and use aluminum L-brackets for added support
on the top plate. The flap is attached using a stainless-steel piano hinge. The upper surface of the
lower plate is positioned 27.6 cm above the bottom of the test section. The airfoil (black, centered
between the plates) is positioned two chord lengths downstream of the leading edge of the plates.
The plates extend 2c upstream from the leading edge and an additional 4.6c beyond the trailing
edge of the airfoil before the flap begins. The shaded green area at the midspan of the airfoil,
represents the wake region of interest. The freestream (U∞ ) is shown in addition to the directions
of the the (x, y, z) coordinate system. The origin is at the mid-span trailing edge of the stationary
airfoil at α = 0◦ .

13

2.2
2.2.1

Airfoil Motion Trajectory
Equipment Description

The airfoil is mounted to a custom1three degrees of freedom electromechanical system that is
capable of providing arbitrary motion trajectories (e.g. pitch, plunge, pitch-plunge, and perching
maneuvers) in the water tunnel.

All of the components are from the Parker Hannifin Corpora-

tion and include: two Trilogy Ironless Positioners with 1 Âľm resolution incremental encoders; a
T4DB0576NPAMA4 and a T4DB0436NPAMA4 for streamwise and cross-stream directions respectively. A MPP Rotary Servo Motor (MPP1154A9D-KPSN) with a single turn, 217 count,
absolute encoder is used for the pitching axis. The entire assembly is mounted above the water
tunnel (Figure 2.3) and raised and lowered into position using scissor jacks at each corner of the
support structure. The current study requires only the use of the rotary motor, so once the airfoil was
centered (to within ±0.3 cm) in the test section, the two linear motors’ carriages were immobilized
with clamps.
Each positioner/motor is driven using a 1300 Watt Aries Drive (AR-13AE) and controlled with
an ACR9000 Controller (9000P1U4M1). A portable chassis, shown in Figure 2.4a, houses the three
drives, controller, and interfacing electronics used for synchronization. Communication with the
controller is over TCP/IP on a private Fast Ethernet connection, whose topology is also indicated in
Figure 2.4b. The controller’s Binary Host Interface is utilized for configuring the motion trajectory
and additional program functionalities in addition to transferring the captured encoder positions
back to the computer for analysis. The software on the central computer (shown in the figure) is a
suite of custom MATLAB programs used for interacting with the controller.
The encoder positions of each axis may be either externally (up to 50 Hz) or internally triggered
at the servo loop period and stored in the onboard memory with a maximum of over 140,000 values.
The memory is read off of the controller and stored on a computer for further analysis. The position
1 The

design, assembly, and initial testing of the system was done by Dr. Bruno Monnier
concurrently with the work of Monnier et al. (2013).

14

Figure 2.3: The three degree of freedom apparatus above the test section. The entire system is
attached to linear bearings on top of a separate (blue) support structure from the test section and
may be positioned at different locations in the streamwise direction. The motor assembly is raised
and lowered to move the airfoil into the test section using scissor jacks installed at each corner of
the frame (not shown).

acquisition program typically records two sets of position data for each data set: (1) an externally
triggered set synchronized with each MTV image pair, and (2) an internally triggered set recorded
at 2 kHz for 5 seconds at the completion of the first set. The two data sets allow proper phase
ordering of the MTV image pair data relative to the motion and for validation of motion trajectory
of each data set. The accuracy of the motion trajectory is discussed in Section 2.2.3.

2.2.2

Prescribed Motion

The motion trajectory considered in the current work is pure sinusoidal pitching motion about the
quarter chord as described by
α = α0 sin (2π f t + φ0 ) + αm,

15

(2.1)

(a)

(b)

Figure 2.4: Electronics chassis and local network topology. (a) Portable chassis that houses
various components: high-powered drives (lower compartment); power switches, indicators, and
emergency stop button (lower front panel); controller, DC power supplies, and I/O signal conditioning (middle compartment); I/O connectors (middle front panel). In the top compartment and
on the top of the chassis two remotely accessed PCs used for image acquisition and controller
configuration (not shown). (b) Simplified diagram indicating various interfacing done with the
controller and computer on a local network.

where α is the instantaneous angle of attack, α0 is the amplitude, f is the frequency, φ0 is a
phase angle, Îąm is the mean angle of attack, and t is the time. The mean angle of attack is zero
for all cases considered, and the phase angle is 180◦ (since the motion starts with al pha going
negative). The prescribed motion is described by two non-dimensional parameters: the reduced
frequency, k = 2π f c/U∞ , and the Strouhal number, St = f h p /U∞ . The trailing edge deflection,
h p , is the peak-to-peak trailing edge amplitude and may be further classified between the actual
and prescribed: h p,m and h p,r , where the former must be measured experimentally and the latter is
computed based on the simple geometric relation for a perfectly rigid airfoil: h p,r = 2¡0.75 sin (ι0 ).

16

For simplicity, h p will be used from here on in place of h p,m .
The single period of the motion trajectory is calculated and stored in a 1000-element-array on the
controller. The time between each element is set based on desired period of the motion. The period
of the motion is always an integer number of milliseconds, with no long-term deviation which the
controller is able to repeat the motion cycle with no long-term deviations. The controller linearly
interpolates the position array as needed at each servo loop period (0.5 ms in these experiments) and
repeats the motion cycle until commanded to stop. Starting and stopping of the motion is ramped
up and down by linearly adjusting the amplitude over 3 motion cycles. The motion, as indicated
by the phase angle of 180◦ starts with the airfoil’s trailing edge moving upwards (negative angle of
attack).

2.2.3

Measured Trajectory

The deviation from sinusoidal motion has been characterized for all data sets used in this study.
The instantaneous angle of attack, measured by the rotary motor encoder, was recorded for each
image pair (synchronized to the first image of the pair) used in measuring the flow field. The dataset
was phase ordered based on the known period of the motion and the image acquisition frequency
(23.93 Hz). The phase ordered position data were fit, using the linear least squares method, to
the sinusoidal function given by Equation 2.1 with α0 , φ0 , and αm as unknowns. Each dataset
was analyzed to determine the deviations from the best fit sine function in terms of the maximum
deviation, amplitude, mean angle of attack, and the RMS deviation from the fitted trajectory as
overall metrics of the accuracy of the motion.
The maximum deviation was less than 0.08◦ and the RMS deviation less than 0.02◦ for 95% of
all motion trajectories considered. The maximum deviation is shown in Figure 2.5 and the RMS
deviation is shown in Figure 2.6 for both the rigid and flexible airfoils. The average deviation from
each of the datasets is plotted with the error bar indicating the 95% confidence interval of the entire
data set which is comprised of approximately 36 sub-fields of view for the rigid airfoil and 40
sub-fields of view for the flexible airfoil. Overall, the prescribed motion trajectory was faithfully

17

Figure 2.5: The average maximum deviation from the best fit sine wave of the motion trajectory.
Abbreviations for the rigid airfoil (R) and flexible airfoil (F) are used in the legend. The error bars
indicate the 95% confidence interval of all of the independent datasets of each motion trajectory.

created experimentally.
In addition to the measured trajectory from the rotary motor’s encoder, the shape of the airfoil
at the bottom tip of the airfoil, closest to the false wall in the middle of the tunnel, was also
recorded. The primary motivation of these measurements was to document the flexible airfoil’s
deformation, as described in Section 2.3.5, but these measurements were also used to measure the
spanwise bending of the airfoil. The spanwise bending of the airfoil needs to be checked to verify
the prescribed 2D motion satisfactory describes the actual experimental conditions since the airfoil
is not being held at the other end. The cross-stream displacement of the quarter chord position was
determined from the shape deformation images and displayed in Figure 2.7. The data show that
with increasing reduced frequency, the amount of non-ideal deflection of the airfoil’s tip increases,
however, only up to a maximum of 0.4% of the chord (0.48 mm).

18

Figure 2.6: The average RMS deviation from the best fit sine wave of the motion trajectory.
Abbreviations for the rigid airfoil (R) and flexible airfoil (F) are used in the legend. The error bars
indicate the 95% confidence interval of all of the independent datasets of each motion trajectory.

2.3

Airfoil Properties

The manufacturing of chordwise flexible airfoil shapes was undertaken using a state of the art
multi-material 3D printer. This technology enables single piece construction of parts with varying
levels of flexibility. In addition to a baseline fully rigid airfoil, regions of the substructure can
be modified to have a flexible leading or trailing edge (or any proportion) while maintaining the
same geometric shape. This is a fundamentally different approach than previous studies, such as:
Monnier et al. (2015), Heathcote & Gursul (2007), and Thiria & Godoy-Diana (2010) that relied on
shape variations (typically by varying the thickness of a material) in order to modify the flexibility,
or in the case of Marais et al. (2012) a rudimentary approximation of an airfoil profile. Utilizing the
3D printer is also simpler than the cast urethane models used by Prempraneerach et al. (2003) and

19

Figure 2.7: Spanwise bending measured at the tip of the airfoil. The cross-stream displacement
of the quarter-chord (pitch axis) is plotted against the reduced frequency. The ordinate is in % of
chord; for reference, 0.4%c is 0.48 mm. The error bars indicate the 95% confidence interval of the
measurement.

Egan et al. (2016) and provides more flexibility in the varying the model’s chordwise flexibility.
The state of the art multi-material 3D printers are not without fault, as there are many nuances that
need to be considered, that will be discussed in the rest of this section, to ensure the final output is
of sufficient quality.

2.3.1

3D Printer Description

The printer, an Objet Connex350, is an inkjet-based 3D printer that uses photopolymers that cure
to form materials ranging from rubber-like materials (with a selectable Shore Hardness values
in the range of 27-95) and plastic-like material similar to PVC. While there are a wide array of
different materials available for the printer (and the materials are under continual improvement by
the company), the two materials used in this study are VeroWhitePlus (RG835), an opaque rigid
material and Objet TangoBlackPlus (FLX980), a rubber-like material. The printer has a resolution

20

of 30 Âľm in the vertical build direction while utilizing these materials and a resolution of 42 Âľm
in the horizontal directions. The maximum print size is 340 mm × 340 mm × 200 mm. In order to
achieve uniform flexibility in the spanwise direction only a single part construction is considered;
which limits the span of the airfoil to 33.7 cm.

2.3.2

Material Testing

The two materials utilized by the Objet Connex 350 printer in the current study are the rigid
opaque Objet VeroWhitePlus (RG835) and the rubber-like Objet TangoBlackPlus (FLX980). While
some material properties are documented and specified from the manufacturer, many of them
are estimated and can vary greatly depending on printer configuration. Additionally, with no
institutional knowledge with regards to their photochemical properties and material compatibility
with our MTV triplex some basic material tests were required. These tests aimed to assess the
feasibility of using these materials in this project by determining their compatibility with our MTV
triplex, the effect of prolonged fluid submersion, and their ability to withstand exposure to high
intensity UV laser light.
Initial tests were performed on a NACA 0012 airfoil profile whose first half of the chord was the
rigid material and the latter half was the pure flexible material. The model had a 12 cm chord and
a span of only 5.1 cm as shown in Figure 2.8. Immediately, several less than ideal ‘features’ of the
printing process became apparent. The first was delamination at the rubber and plastic interfaces
as shown in the image. The second was a considerable “drooping” of the rubber portion near the
trailing edge (not visible in the image) in the downward direction relative to when it was printed;
e.g. gravity appears to pull the flexible material down slightly during curing.
Initial structural response tests involved pitching the airfoil over a large range of frequencies
and amplitudes to address its deformation relative to a rigid airfoil. These initial tests indicated
that the NACA 0012 model’s trailing edge only marginally (10-15%) deformed more than a rigid
airfoil. Increasing the aeroelastic effects would require either a more flexible material, which is not
available using this 3D printer, or a thinner airfoil; to this end, a NACA 0009 profile was designed

21

Figure 2.8: Test specimen of 3D printer. The airfoil has a chord of 12 cm and a span of 5.1 cm. A
tapped 1/4-20 hole at the 0.25c location was used to pitch the airfoil. The dark rectangular marks
on the plastic (white) surface are burn marks from an unfocused 308 nm laser beam. The small
chuck of rubber (black) is detached and the interface plane between the rubber and plastic parts can
be seen to have delaminated.

for use in this study.
The photochemical properties and the material’s interaction with the MTV triplex was studied.
The surface of the part was imaged as a pulsed Excimer laser beam (308 nm) was directed at the
surface while it was in air. As shown in the image in Figure 2.8, the plastic material readily left
burn marks after very few shots. The rubber material did not readily burn, although the surface of
the material did heat up over time. No fluorescence or phosphorescent emission from the materials
were visually observed. The model was submerged in the MTV triplex (see Section 2.4 for details
on the composition of the solution) for a several minutes, wiped dry, and re-tested in air. The plastic
material readily emitted a green phosphorescence, the rubber material did not have any noticeable
emission.
To further quantify the emission from the material after its exposure to the MTV solution, the
test was repeated but the results were quantified with a simple camera setup, recording images
at different time delays relative to the laser pulse. The average intensity of the spatial region in
which the laser beam was impacting the surface was used to estimate an exponential lifetime of the
observed emission. The TangoBlackPlus material shows an exponential decay lifetime of τ=0.3 ms,
whereas the VeroWhitePlus material had a lifetime of τ=3.5 ms (see Figure 2.9). The intensity
between two tests were directly comparable, indicating that not only is there a considerable lifetime,

22

Figure 2.9: Photochemical properties of 3D printer materials after exposure to MTV triplex. The
spatially averaged intensity of a region of the 3D printed airfoil excited by a 308 nm pulsed laser at
different time delays relative to the laser pulse. Data has been normalized by the RGD835 (rigid
opaque) material.

but the intensity is initially over 10x higher with the rigid material. A reasonable explanation to the
behavior observed is that the tracer, 1-Bromonathalene, has been absorbed by the plastic material
but not the rubber material. The implications are that for near-surface measurements, the airfoils
should utilize a rubber coating on the surface.
An early observation in the material also indicated that the rubber material would increase
in flexibility over time as it was continually submerged. The flexibility of the material tends to
increase several percent (measured by comparing the trailing edge deflection for a particular motion
trajectory) every few hours. This trend continues for the first couple of days of submersion and
only begins to stabilize after a week or longer of continual submersion. This effect is shown in
Figure 2.10, which was characterized early in the material testing phase of the NACA 0012 profile.
The maximum change in h p at k = 12 being only 12%, motivated the future use of the NACA 0009
profile to increase the effect of flexibility. Several of the test conditions (k=22, and 32) far exceed

23

the conditions that will be used in this study, but are included as they were the only conditions to
be tested at the 21 day submersion test. The lower k cases all appear to behave the same, with
regards to the h p increasing with submersion time at approximately the same rate up until the 7
day submersion test, so the conclusions reached with this test should be transferable to the lower k
regime that will be considered in this study.

Figure 2.10: Trailing edge deflection as a function of the submersion time. The ordinate is the
percent change in the trailing edge deflection of the flexible airfoil (h p− f lex ) compared to the
computed trailing edge deflection of rigid airfoil (h p−rigid ). These experiments were conducted on
the NACA 0012 airfoil shown in Figure 2.8

Within a couple days of submersion, the airfoil’s rubber coating would begin to blister and
delamination between the rubber and plastic parts would start to develop. With the duration required
for achieving a constant level flexibility is over a week, an alternative method for conditioning the
material to achieve a constant level of flexibility is required. The process developed started with
precondition the airfoil in a saturated air environment for several weeks prior to the experiments,
which softened the material without causing blistering. The airfoil was submerged in the water
tunnel prior to the experiments for two hours before flowfield measurements were performed. The

24

flowfield measurements on each day consisted of two sub-fields-of-view (sFOV) for which the
entire parameter space was repeated twice to increase the overall data resolution (see Section 2.4.3
for complete details). The airfoil was raised out of the water and wrapped in plastic to keep the
airfoil wet at the end of the second sFOV approximately after 6 hours of continuously running
experiments or 8 hours of total submersion time. The process was repeated daily until the entire
flowfield was acquired.
Verification that this process produced a constant level of flexibility throughout the experiments
is quantified by measuring the changes in the trailing edge deflection, h p , which was measured
throughout the experiments on the flexible airfoil. Complete details of the shape tracking methods
and results are presented in Section 2.3.5. The first metric that will be considered is the change in
flexibility between the two sFOVs acquired on a single day. Figure 2.11 shows the percent change
in h p from the first to the second sFOV. The average change between the two sFOVs is plotted
against k with error bars indicating the 95% confidence interval of all of the FOVs. The average
was computed from the 21 pairs of sFOVs used in mapping the entire flowfield. The second sFOV
produces, on average across all k and Îą0 , a 0.12% higher h p than the first sFOV of the day. Overall,
the deviations in h p between the two sFOVs acquired on the same day is less than 1.5%. Figure 2.12
depicts the percent change in the average h p measured for all (42) sFOVs as error bars representing
the 95% confidence interval. The variation in h p increases with increasing k but remains within 4%
for all motion trajectories over the entire duration of the flexible airfoil experiments. The α0 = 2◦
cases vary more than α0 = 4◦ and α0 = 6◦ , both of which see less than 3% variation for all cases.
In physical units, the deviation in h p by 3% for the largest h p observed at k = 6 and α0 = 6◦ is
0.7 mm, while the maximum deviations of 4% for the highest h p observed for α0 = 2◦ at k = 12
is 0.4 mm. In summary, the change in flexibility throughout the course of the measurements has
been well controlled with the procedure outlined and the change in flexibility, as measured by h p
has been documented. Further analysis and the full presentation of the measured h p is located in
Section 2.3.5.

25

Figure 2.11: Change in the trailing edge deflection from the first sFOV of the day to the second
sFOV. Error bars indicate a 95% confidence interval of all of the FOVs in the flexible airfoil results
(21 FOVs total).

2.3.3

Design and Shape Validation

The internal structure of the airfoil was designed to address several of the issues discovered during
early testing of the 3D printer’s materials. The design provisions included: (1) use of a rubber
coating to prevent surface emission, and (2) a remedy for the observed drooping of the trailing
edges. A cross sectional view of the rigid and flexible airfoils can be found in Figure 2.13. The
chord length, c, of the airfoil is nominally 12 cm in length.
The first of the provisions listed above was implemented by including a 0.51 mm thick shell
of a rubber coating on the entire outer surface of the airfoil; which was the minimum thickness
that could cure properly to the plastic material without delamination issues. The second provision

26

Figure 2.12: The change in the trailing edge deflection across all FOVs of the flexible airfoil.
The error bars indicate a 95% confidence interval of h p from the average h p of all FOVs (42 total
sFOVs).

Figure 2.13: Cross-sectional views of the 3D printed airfoils. The rigid (top) and flexible (bottom)
NACA 0009 airfoils printed on the 3D Printer. The flexible TangoBlackPlus (FLX980) material is
shown in black and the rigid VeroWhitePlus (RG835) material is gray. The internal void in the part
is where a stainless steel shaft (1.000 in × 0.250 in with 0.080 in chamfers) is centered at 0.25x/c
and is used throughout the span to interface the airfoil with the pitch motor.

27

entailed the introduction of a piece of solid plastic material in the last 14.2% of the chord. This
added structure helps preserve the trailing edge shape, and since it runs the entire span, the structure
also improves the two dimensionality of the chordwise flexing. However, the very tip of the trailing
edge (0.01c or 1.2 mm) was still deformed in both the rigid and flexible airfoils and was removed,
leaving the airfoil with a slightly blunt trailing edge of 0.4 mm, compared to the theoretical trailing
edge thickness of 0.23 mm.
The airfoil cross sections shown in Figure 2.13 are colored to indicate the material used in each
of the different regions. The void in the airfoils is centered at the quarter chord to make room for
a 1.000 in×0.250 in stainless steel shaft with 0.080 in chamfers which runs through the entire span
of the airfoil and connects, via several connectors and adapters to the pitching motor. The flexible
material starts to make up the airfoil at 0.4c, which can be seen most easily in the top of Figure
2.14.
The actual shape of the airfoils (when not moving) was measured by imaging the surface of the
airfoil at midspan. A thin (2 mm) laser sheet illuminated the surface so it could be imaged and
manually located with an accuracy of Âą0.08 mm (Âą0.07% of the chord). The top and bottom surface
were imaged in this fashion and combined together to estimate the airfoil thickness distribution
shown in Figure 2.14. This is the thickness distribution of the airfoil relative to the chord line.
The thickness distribution of the rigid and flexible airfoils agree very well up until x/c=0.35
where the flexible airfoil starts to lose thickness compared to the rigid airfoil as the internal material
starts to rapidly transition from the rigid to rubber material. The thickness of the two airfoils also
agree very well again after x/c=0.75. Between x/c=0.35 and x/c=0.75 the difference in thickness
between the rigid and flexible airfoil has a maximum value of 0.005c or 0.65 mm compared to
the airfoil’s maximum thickness of 10.8 mm. The overall shape of the both the rigid and flexible
airfoils along with the theoretical profile, in an exaggerated scale, in Figure 2.15. Overall, the
rigid and flexible airfoils compare well against the NACA 0009 airfoil, especially towards the
leading edge, but both have slight deviations compared to the theoretical profile and to each other.
Direct comparisons will be made in Section 3.1.1 between numerical simulations of a NACA 0009

28

0.1
0.08
0.06
0.04
0.02
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 2.14: Measured thickness distribution of the rigid and flexible airfoils. The thickness
distribution measured at the mid-span. The internal structure of the flexible airfoil is shown above
the plot for reference.
airfoil and the rigid airfoil measurements that will provide additional evidence that the measured
differences in airfoils static shapes are not significant.

Figure 2.15: Measured airfoil profile of the rigid and flexible airfoils. The measured airfoil cross
section (at mid-span) of the rigid and flexible airfoils compared to the theoretical NACA 0009
airfoil section

29

2.3.4

Estimated Flexibility

Non-dimensional analysis applied towards the effect of flexibility has been well discussed in the
literature; with the typically two non-dimensional parameters surfacing, one related to the ratio of
the densities between the airfoil and its working fluid, and the other relates the bending rigidity
of the airfoil to the aerodynamic forces. Heathcote & Gursul (2007) define the bending stiffness
parameter using a method they attribute to being similar to that of Murray (2000). Their bending
stiffness parameter is
Îťx =

(E I) x
2 Sc2
0.5ρU∞

,

(2.2)

where (E I) x is the flexural stiffness in the chordwise direction, S is the planform area, and c is
the chord. This is similar, as they discuss in Gursul et al. (2014) to the “effective stiffness”, Π1 ,
used by Shyy (Shyy et al., 2010; Kang et al., 2011). The primary difference between the two is Π1
written explicitly for an isotropic plate using the classical thin plate bending theory for the flexural
stiffness. The Π1 term is given by
Π1 =

Et 3 /(12(1 − ν 2 ))
2 c3
ρU∞

,

(2.3)

where E is again Young’s modulus, ν is Poisson’s ratio, and t is the thickness of the plate. The
density of the fluid (ρ), the freestream velocity (U∞ ), and the chord length (c) are also used in the
Π1 term. In the current study, both terms will be estimated for comparison, but the first term, Ν x , is
more directly applicable as the current study uses relatively thicker airfoils than the flat plates Π1
was written for.
The flexible material used in this study, TangoBlackPlus (FLX980), has a manufacturer supplied
data on the Shore Hardness value of 26-28 on the A scale (Stratasys, 2017). Converting between
Shore Hardness and Young’s Modulus (E) is not straightforward. A recent method for converting
between Shore Hardness scales and E was developed in Mitchell et al. (2011), the details of
which are excluded for brevity, leads to Young’s Modulus estimates of E=1.3±0.06 MPa for the
FLX980 material. The uncertainty range is based strictly on the Shore Hardness range specified
by the manufacturer. The increase in flexibility of the material that was observed due to prolonged
30

Îťx
Rigid
1.22×103
Flexible
1.04

Π1
1.45 × 103
2.30

Table 2.1: Non-dimensional bending stiffness parameters of rigid and flexible airfoils. Estimated
non-dimensional bending stiffness of the rigid and flexible airfoils.
submersion (see Figure 2.10) most likely decreases this value of E by 10%-20%, but no further
material testing were performed to improve the estimates of this value. Mitchell et al. (2011) state
that the Poisson’s ratio for rubber is well estimated by ν=0.5.
The manufacturer (Stratasys, 2017) provides the Young’s Modulus of the rigid material used
in this study, VeroWhitePlus (RGD835), as E=2 GPa-3 GPa. A value of 0.4 is used for Poisson’s
ratio, which is a common value for nylon. The average thickness of airfoil, 0.0411c or 4.93 mm, is
used as the thickness of the plate for computing Π1 . The chord length used for the flexible airfoil is
0.6c, which is the flexible portion of the airfoil (as shown in Section 2.3.3). When computing for
the rigid airfoil, the entire chord length is used for c in the equations of Ν x and Π1 . The average
thickness of only the flexible portion is used when computing the area moment of inertia for the
flexible airfoil, while the average thickness of the entire airfoil is used for the rigid airfoil. The
average value in the range for E will be used in estimating Ν x and Π1 .
The computed values for non-dimensional bending stiffness are summarized in Table 2.1. The
Π1 parameter for the rigid airfoil is estimated to be Π1 = 1.45 × 103 for the rigid airfoil, and
Π1 = 2.3 for the flexible airfoil. The rigid airfoil in the study by Monnier et al. (2015) was Π1 = 97
and their flexible airfoils were in the ranged of 0.6 < Π1 < 5.3. In the study by Heathcote & Gursul
(2007), λ x varied from 0.07 to 260, compared to λ x = 1.04 for the flexible and λ x = 1.22 × 103
for the rigid airfoil of the current study. This puts the current study within the range of “moderate”
flexibility where it would be expected, based on these and other studies, to experience increased
thrust and efficiency compared to the rigid airfoil.

31

2.3.5

Deformation Tracking

The shape of the flexible airfoil is monitored throughout the experiments by imaging the tip of
the airfoil nearest to the bottom of the tunnel. This was accomplished by installing a Point Grey
Grasshopper3 GS3-U3-28S4M camera with a 1928 x 1448 pixel CCD underneath the tunnel as
pictured in Figure 2.16. This camera setup utilizes a mirror mounted at 45◦ relative to the bottom
wall of the test section such that it imaged the end of the airfoil. The airfoil’s tip itself was painted
with a white rubberized coating (Rust-Oleum LeakSeal) to create contrast between the airfoil and
the background as shown in the raw example image in Figure 2.16. The tip of the airfoil was
illuminated by an LED flashlight that was externally focused and the beam was blocked such that
it only illuminated the end of the airfoil and did not interfere with the velocimetry measurements.
The camera was synchronized with the rest of the imaging system during the experiments.
The trailing edge displacement, h p , was measured from the shape deformation images for
each dataset to ensure the day to day repeatability of the flexible airfoil’s flexibility, as discussed
previously. The actual trailing edge deflection measured is compared to that based on the prescribed
motion of the rigid airfoil in Figure 2.17. The trailing edge deflection normalized by the calculated
(rigid) trailing edge deflection based on the prescribed motion (h p,r ) is shown in Figure 2.18. For all
amplitudes the ratio of (h p − h p,r )/h p,r increases with reduced frequency at the same rate until k=4,
where the rate of α0 = 6◦ begins to decrease, while the cases with α0 = 4◦ stay with the α0 = 2◦
cases until k=6, where its rate of increase too begins to decrease. For the α0 = 2◦ trajectories, there
appears to be a peak in the trailing edge deflection at k=12, after which the trailing edge deflection
is reduced with increasing k.

32

Figure 2.16: Shape monitoring camera setup and sample image. The shape monitoring camera
setup (top) is installed underneath the test section. The airfoil is not shown in the within the test
section. The sample image (bottom) is for the flexible airfoil oscillating at k=8 with α0 = 6◦ . The
shadows and bright spots in the bottom image is primarily due to non-uniform illumination and the
applied coating of Rust-Oleum LeakSeal bleeding slightly around the edges of the airfoil.

33

Figure 2.17: Trailing edge displacement for all motion trajectories of the flexible airfoil. The solid
lines indicate the trailing edge displacement of the rigid airfoil, with the color indicating the Îą0 .
The error bars indicate a 95% confidence interval of h p from all FOVs.

34

Figure 2.18: Relative change in trailing edge displacement. The percent change in the trailing
edge deflection for all motions of the flexible airfoil. The error bars indicate a 95% confidence
interval of h p all FOVs.

35

2.4
2.4.1

Molecular Tagging Velocimetry
Background

The current experiments utilize molecular tagging velocimetry for measurements of the velocity
field of the oscillating airfoils. Molecular Tagging Velocimetry (MTV) is a whole-field optical
technique in which naturally occurring or premixed molecules (or atoms) in a flowing medium
are turned into long lifetime tracers through photon excitation of an appropriate wavelength.
Typically, a pulsed laser is used to excite or ‘tag’ regions of interest and the resulting emission
is interrogated twice with a prescribed time delay to form an image pair. Spatial gradients in
intensity form regions of interest in the image pair that are correlated with one another to yield a
measurement of the Lagrangian displacement. Combined with the prescribed time delay between
the images, the velocity field is determined as demonstrated in Figure 2.19. The reader is referred to
Gendrich & Koochesfahani (1996) for details on the image processing procedures, Gendrich et al.
(1997) for information on the phosphorescent supramolecule used in the current application, and
Koochesfahani & Nocera (2007) for a comprehensive review of the development, photochemistry,
and several applications of MTV.

Figure 2.19: Example 2c-MTV image pair and resulting 2-component velocity vector field. (a)
The fluid flow of a vortex ring impinging on a flat plate with normal incidence is imaged 1 Âľs
after the laser pulse. (b) The same tagged fluid is interrogated again 8 ms after the first image. (c)
The resulting 2-component velocity vector field is derived from correlating the undelayed (a) and
delayed (b) images. (Gendrich & Koochesfahani (1996)) With permission of Springer

36

2.4.2

Current Implementation

In the current implementation, an excimer laser, a Coherent COMPexPro 205 XeCl (308 nm) is used
as the excitation source. The optics and beam path used to create the tagging pattern are depicted in
Figure 2.20. The optical assembly enables the tagging pattern to be placed in the test section along
the streamwise, cross-stream and vertical directions independently. The current measurements are
limited to a single horizontal plane at the midspan of the airfoils.

Test Section
Airfoil
(not to scale)

Beam Blocker

Beam
Expander

Mirror
50/50 Beam Splitter

Translating
Optics Platform

From Laser

Sheet Forming Optics

Figure 2.20: Top view of optical beam path of excimer laser. At the location of the red circle, indicated as the sheet forming optics, the beam is also directed vertically upwards from approximately
41 cm off of the floor, to the midspan of the airfoil, approximately 163 cm from the floor. As the
beam travels vertically, it passes through the sheet forming optics. The airfoil depicted in the test
section is 2.5 times larger than its actual size relative to the width of the test section width.

37

Referring to Figure 2.20, the beam path travels parallel to the test section before being directed
vertically by a mirror that moves with the assembly only in the streamwise direction. The beam
travels vertically to a set of sheet forming optics which move both in the streamwise and vertical
direction with the rest of the optics platform. The sheet forming optics consist of a pair of lenses
that form a very long focal length lens, creating a sheet thickness of <1 mm in the test section.
The vertically traveling laser sheet passes through another mirror that turns the beam into the
cross-stream direction on the measurement plane where the rest of the optics are located on a
large wooden platform. The wooden platform moves with the overall system in the streamwise
and vertical directions but is also moveable in the cross-stream direction to aid in repositioning the
tagging pattern. While traveling on this platform the sheet first intersects a 50/50 beam splitter each
of the two sheets are directed through 3:2 beam expanders in the horizontal plane (making the laser
sheets ~51 mm × 1 mm in cross-section). Lastly, the sheets are directed into the test section and
pass through brass beam blockers that have vertical slots to form the desired tagging pattern whose
lines are approximately 0.6 mm in width and less than 1 mm thick in the spanwise (z) direction of
the airfoil. Typically, there are 30 lines spaced approximately 1.4 mm apart from each laser sheet,
which intersect to form approximately 800 intersections in each sFOV.

2.4.3

MTV Imaging System

The imaging system used for the current work utilizes two 1392×1024 pixel PCO.Pixelfly cameras
which are mounted underneath the water tunnel as shown in Figure 2.21. The camera assembly is
connected to the streamwise and vertical traverses such that it also moves with the entire optical
assembly. It is on a separate cross-stream traverse since it must move at a different rate than the
optical system due to the refraction of the UV beam passing through the quartz windows at an
angle. A single 105 mm f/2.0 Nikon Nikkor Lens, a shortened F-mount to C-mount adaptor, and
a standard 1” C-mount cube beam splitter from Edmund Optics completes the setup. The compact
two-camera system relies on a calibration routine based on the Affine transform to compensate for
large (greater than 10 pixel) sensor misalignments, which is described in detail in Appendix A, but

38

yields an uncertainty of less than 0.02 pixels RMS in calculating the displacement.

Figure 2.21: Photograph of dual camera MTV imaging system. Two PCO.Pixelfly cameras (blue)
are flush mounted to a backing plate and a cube beam splitter (highlighted in red). A single lens is
mounted to a shortened c-mount to Nikon F-mount adapter (highlighted green). The trigger boxes
can be seen in the upper right of the image. The cameras are located underneath the water tunnel
looking up towards the test section. The gray pipe in the bottom of the image is the return pipe of
the water tunnel.

39

2.5
2.5.1

Experiments
Overview

The current study focuses on the wake structures of an oscillating rigid and flexible airfoil. The
experimental study of Bohl & Koochesfahani (2009) will be used as a baseline comparison, and
as such, the overall region of interest and required measurement resolution can be used as a
guide for the current study. The measurement domain required to study the wake is nominally
−0.1 < x/c < 1.5 in the streamwise direction, with the origin located at the location of the trailing
edge of the stationary airfoil at α = 0◦ . In the transverse direction, the domain of interest is
Âą0.3c over the entire streamwise extent of the domain but a larger extent of Âą0.6c is required at
x/c = 1 where a control volume method for estimating the streamwise force (drag/thrust) will be
applied. This region was shown previously in Figure 2.2. The required spatial resolution, based
on the expected flowfield as measured by Bohl & Koochesfahani (2009), is approximately 0.01c,
to accurately measure the expected vorticity field.
Multiple fields of view (FOVs) are required to perform the measurements at the desired resolution over the entire domain. The FOVs positions for the rigid airfoil measurements are shown
in Figure 2.22, and Figure 2.23 for the flexible airfoil measurements. At each FOV displayed in
the figures two sub-FOVs (sFOVs) are required to achieve the required spatial resolution. At these
sFOVS the cameras do not move, but the laser tagging pattern is shifted by nominally one half
of the grid spacing. The final velocity vector spacing is nominally 0.008c in both the streamwise
and transverse directions. During the flexible airfoil experiments, a second set of FOVs centered
at x/c=1 were acquired to have a complete redundant set of data for applying the control volume method (not shown on Figure 2.23). The camera’s pixels were 2×2 binned on the sensor to
improve image quality and acquisition rate; which was 23.93 Hz. Each FOV was approximately
5.5 cm × 4.1 cm (width × height) with a scale factor of 125 pixels/cm.
The experiments were carried out one field of view at a time. A calibration image was acquired
for each field of view and is used to compute the Affine transformation for the two camera calibration,

40

Figure 2.22: Fields of view for the rigid airfoil. The dashed yellow lines are located at Âą0.6c, the
dashed red lines are located at Âą0.3c and the dashed vertical line is located at x/c = 1c. The airfoil
is shown as black at an α = 0◦ , with its trailing edge located at (x, y) = (0, 0).

the image scale factor, and the sub-pixel tracking of the relative movement of one FOV to the next.
To conduct the experiment, the airfoil starts from rest at α = 0◦ and begins its motion. At least
50 convective times are allowed to pass after the start of each motion before measurements begin.
After the conclusion of all measurements (described in the following paragraph) the motion is
stopped and an additional 50 convective times are allowed to pass before the controller gains and
settings are updated for the next set of motion parameters. The stop/start procedure was automated,
but had to be monitored to ensure synchronization between the four separate acquisition systems.
The laser grid pattern is manually shifted by approximately one half the grid spacing after all cases
have been completed and the process is immediately repeated for the second sFOV.
The non-dimensional motion parameters for each of the airfoils are outlined in Figure 2.24,

41

Figure 2.23: Fields of view for the flexible airfoils. The dashed yellow lines are located at Âą0.6c,
the dashed red lines are located at Âą0.3c and the dashed vertical line is located at x/c = 1c. The
airfoil is shown as black at an α = 0◦ , with its trailing edge located at (x, y) = (0, 0).

and represent 43 distinct trajectories for the rigid airfoil and 33 distinct trajectories for the flexible
airfoil. The oscillation frequency ranged from 0.2 Hz to 4.8 Hz, and with amplitudes of α0 = 2◦ ,
4◦ , and 6◦ . This corresponds to a reduced frequency range of 0 < k < 18 and Strouhal number
range of 0 < St < 0.6. The entire parameter space considered in this study is displayed in Figure
2.24, plotting k against the St of the prescribed motion. The mean angle of attack (Îąm ) was found
to within ±0.1◦ by considering the symmetry of the velocity profile in the wake of the rigid airfoil.
The freestream velocity of the tunnel was set to 10 cm/s, as verified by independent 1c-MTV
and 2c-MTV measurements at various parts of the test section (both with and without the airfoil
present). The temperature of the water tunnel is monitored daily throughout all of the experiments
and was in the range of 20.3◦ C and 23.3◦ C and is used in the computation of the viscosity based on

42

Figure 2.24: Summary of all motion parameters. The entire motion parameter space considered
for the rigid and flexible airfoils.

the temperature dependence, which was measured to be the same as deionized water. Combined
with a chord length of 12 cm yields a chord Reynolds number range of Rec = [11,900, 12,800] due
to the viscosity change as a function of temperature.
The data set acquired for each case at each sFOV comprises of 2048 image pairs consisting of
an undelayed image acquired < 5 Âľs after the laser has fired and a delayed image acquired between
0.5 ms and 4 ms after the laser pulse (depending on the FOV and motion parameters). The exposure
of the two cameras were identical, typically 1/10th of the delay time (select cases with very short
delay times required longer exposures, but never exceeding an exposure to delay time ratio of 1/5).
The airfoil’s pitch angle is also recorded at the same time of each undelayed image. For the flexible
airfoil experiments, the airfoil shape monitoring camera also records an image synchronized with
each undelayed image.
The images are pre-processed by applying the Affine transform to the undelayed image se-

43

quence. The images are then processed using in-house cross-correlation program (Gendrich &
Koochesfahani, 1996), the residual displacement field from the Affine transformation is subtracted,
bad vectors are removed, and finally, the velocity field is phase-averaged relative to the airfoil
oscillation cycle into 64 phase bins, such that there are 32 velocity samples per phase bin. The
phase averaged velocity fields are combined from each sFOV and are remapped onto a regular
grid and the vorticity is calculated using a central 4th-order finite difference method (Cohn &
Koochesfahani, 2000). The processing requirements are as follows for a single motion trajectory
at each sFOV: (1) pre-processing: 180 s, (2) cross-correlation: 2100 s, (3) two camera correction:
300 s, (4) data filtering: 120 s, and (5) phase averaging: 15 s, for a combined processing time of
approximately 45 minutes. Combining all of the FOVs together and remapping onto a regular grid
takes an additional 30 minutes or processing time per motion parameter.

2.5.2

Phase Averaged Data

Phase averaging is a useful measurement technique well suited to oscillatory flows where a single
frequency (or its harmonics) dominate the flow. This is evident by considering the measured
instantaneous velocity as
u(x, y, t) = u(x, y) + u0φ (x, y, φ) + u01 (x, y, t) + u02 (x, y, t) + u03 (x, y, t)

(2.4)

where u is the time average, u0φ is the oscillatory velocity fluctuation, u01 is the cycle to cycle
variation, u02 is the random fluctuation from turbulence, and u03 is the random fluctuation from
measurement error. This is the same as the triple decomposition of Hussain & Reynolds (1970), but
the third random term has been split into several distinct contributions. When u0φ is significantly
larger than u01 and u02 then phase averaging will well represent u(x, y, t). Phase averaging will
thereby also reduce the rest of the random velocity fluctuations; which includes the real cycle-tocycle variations, turbulence, and measurement noise. Therefore, in oscillating airfoil experiments,
phase averaging is commonly employed.

44

The phase of a particular velocity measurement may be computed by
φi = i ×

f
facq

, where i = 0, 1, 2, ... (M − 1)

(2.5)

where M is the number of samples in the acquisition (2048) and f is the frequency of the motion
and facq is the acquisition rate of the data (23.93 Hz). The phase is reduced to the range of [0,1)
by simply subtracting the integer number of cycles that have already occurred, i.e. φ = φi − bφi c.
Phase averaging is the process of averaging the u(x, y, t) over a finite range of the the phase
(φ1, φ2 ) where φ1 = φ − ∆φ/2 and φ2 = φ + ∆φ/2. Care must be taken with selecting facq and f
to ensure each phase bin is equally populated. The size of the phase bin, ∆φ is 1/64 for the current
measurements, leading to having 64 phase bins whose center is φn = n(∆φ) − (∆φ)/2 for n=1 to
64. The phase averaged streamwise velocity can be considered
hui (x, y, φn ) =

1
∆φ

φn∫+∆φ/2




u(x, y) + u0φ (x, y, φ) dφ,

(2.6)

φn −∆φ/2

in its integral form. The phase-bin-rms velocity, huirms is defined as RMS of the instantaneous
velocities that contributed to a single phase bin’s hui (x, y, φn ). It is a measurement of the overall
fluctuations measured in the experiment, including cycle-to-cycle variation, turbulence, and measurement noise. If the phase bins are too large, the natural variation in the underlying value would
contribute significantly to the phase-bin-rms as well.
The phase averaged variables are used to compute the overall mean velocity profiles hui and
the fluctuating velocity profiles hui0. While hui is equal to u, hui0 will be less than overall temporal
RMS since the phase averaging process will remove some random noise (and some real components
of the flow as discussed earlier).

2.5.3

Sample Measurements

A sample set of measurements of the rigid airfoil oscillating at a reduced frequency of k = 8,
and a pitching amplitude of α0 = 2◦ (St = 0.13) are shown in Figures 2.25-2.27. The mean and
fluctuating streamwise velocity are shown in Figure 2.25, while the transverse velocity component
45

and its fluctuation are shown in Figure 2.26. A single vorticity field, extracted from φ = 0.5 is
shown in Figure 2.27 along with the uncertainty of the vorticity field. Further discussions of the
uncertainties of the measured or calculated variables are included in the next sections and Appendix
B.

46

(a) Mean Streamwise Velocity

(b) Streamwise Velocity R.M.S.

Figure 2.25: Representative mean and rms streamwise velocity field. The calculated mean streamwise velocity (top) and its R.M.S fluctuation (bottom) for the rigid airfoil oscillating at k = 8 and
α0 = 2◦ (St = 0.13). The black region in the lower left and lower center portion of the plots are
regions where there was no velocimetry data.

47

(a) Mean Transverse Velocity

(b) Transverse Velocity R.M.S.

Figure 2.26: Representative mean and rms transverse velocity field. The calculated mean transverse
velocity (top) and its R.M.S fluctuation (bottom) for the rigid airfoil oscillating at k = 8 and α0 = 2◦
(St = 0.13). The black region in the lower left and lower center portion of the plots are regions
where there was no velocimetry data.

48

(a) Spanwise Vorticity

(b) Spanwise Vorticity Uncertainty

Figure 2.27: Representative vorticity field and its uncertainty. The calculated vorticity (top) and its
corresponding uncertainty (bottom) for the rigid airfoil oscillating at k = 8 and α0 = 2◦ (St = 0.13)
extracted at φ = 0.5. The cuttoff level in vorticity (top) is −3 ≤ hωz i ≤ 3 and there is no cutoff
level used in the uncertainty contour (bottom). The black region in the lower left and lower center
portion of the plots are regions where there was no velocimetry data.

49

2.5.4

Control volume method for estimating streamwise force

The control volume method for computing the streamwise force in the current study was described
by Bohl & Koochesfahani (2009), whose inclusion of the fluctuating velocity profiles in the analysis
demonstrated a marked improvement compared to the traditional analysis of using only the mean
velocity profile. They found that only using the mean velocity profile would tend to overestimate
the thrust produced by an oscillating airfoil. The control volume method is
!
!
(
∍
hui
2 +H hui hui
−1 +ε
−1
CF =
c −H U∞ U∞
U∞
!)
 0 2  0 2
Uo2
hui
hvi
1
+
−
+
1−
δy,
2
U∞
U∞
2
U∞

(2.7)

where


1
U0
Îľ=
1−
,
2
U∞

(2.8)

where U0 is the average streamwise velocity at the extrema of the domain. In the current experiments
U0 is estimated by averaging 0.05c at the extrema of the streamwise velocity profile. The integral
form of the equation is approximated using the numerical quadrature method with the current
measurements at 1c downstream of the airfoil, with the the integration performed over the transverse
region of ÂąH = Âą0.6c from the center of the velocity profile.
A detailed analysis of the experimental uncertainty in CT is relegated to Appendix B.4, where
standard uncertainty analysis is used to propagate the experimental uncertainties of all of the
measured quantities in Equation 2.7. Also included in this Appendix is a comparison of the CT
computed from each of the two independent sets of measurements acquired for the flexible airfoil.
The uncertainty was, as expected, a function of the motion parameters, but the averaged 95%
confidence interval for all cases was CT Âą 0.01, with the uncertainty increasing with increasing
Stte . From the repeated measurements for the flexible airfoil, all but 3 cases had differences in CT
less than CT = 0.01 .

50

y

a

U∞

x
yc,p
hp
yc,n

+
–

+
–

b

Figure 2.28: Schematic of typical vortex pattern in wake of oscillating airfoil. A reverse von
KĂĄrmĂĄn vortex street with the transverse spacing (b), wavelength (a) indicated. The trailing edge
displacement (h p ) is also depicted of the airfoil pitching about its quarter chord. The sign of the
vortex is indicated by the arrow and by the (+) or (-) sign.

2.5.5

Wake Vortex Properties

The properties of the vortices in the wake will be investigated in Chapter 4. The vortices of the
rigid and flexible airfoils will be compared in terms of their strength, size, and spacing. In this
section, these parameters will be defined and how they are calculated from the phase averaged
experimental data will be described. These properties are extracted at x/c=0.5 downstream from
the trailing edge of the airfoil. Only motion parameters producing a single pair of vortices, one of
each sign, will be considered for analysis. An example of such a case was shown in Figure 2.27a
and a schematic of this vortex configuration was shown previously in Figure 1.1, but is repeated in
Figure 2.28 for convenience.
The vortices are located in the wake by monitoring the vorticity profile at extraction location,
x/c = 0.5, throughout the entire phase of the motion. The maximum and minimum vorticity are
found, recording their transverse location (yc,p, yc,n ) and the phase (φ p , and φn ) at which these



events occurr. The peak vorticity, ω pk , are taken directly as these maximum and minimum
vorticity values, respectively. The transverse location of the vortices are refined using a second
order polynomial fit to the vorticity profile in the y-direction using Âą2 grid points around the initial
estimate of the peak location. The transverse spacing, b, is now computed as b=yc,p − yc,n . A
51

conservative estimate of the uncertainty in b is Âą0.005c, which is half the grid spacing.
The convective speed of the positive and negative vortices are computed by repeating the
previous procedure at 0.25 ≤ x/c ≤ 0.75. The phase at which the maximum (or minimum)
vorticity arrives at the monitoring station is plotted against the monitoring station. A line is fitted to
the extracted data and the inverse of the slope, ∆x/∆φ, is multiplied by the frequency of the motion,
f , to yield the convective velocity, Uc . The streamwise spacing, or wavelength (a) is computed by
a = Uc / f .

(2.9)

Further analysis of each (+) and (-) vortex is now performed. A vorticity cutoff is applied



to the entire phase averaged vorticity field such that  ω pk  < 2s−1 is set to 0. This helps
remove background noise in the vorticity field from effecting the integrated terms of the core
size and circulation (to be defined below). The same vorticity cutoff limit was used by Bohl &
Koochesfahani (2009) and Hammer (2016).
The circulation and core size, based on the radius of gyration, are computed from
Γ0 =
and

∏
hωi dA,

(2.10)

r 2 hωi dA
∏
,
hωi dA

(2.11)

v
t∏
rc =

respectively. The integrals, estimated by numerical quadrature, are evaluated only within a radial
distance of r=0.0833c from the center of the vortex to facilitate direct comparisons with previous
studies (Bohl & Koochesfahani, 2009; Hammer, 2016). Additionally, only the same sign of vorticity
as the vortex currently being analyzed is included in the integration. The effect of the radius of
integration and vorticity cutoff on Γ0 and rc are discussed in Appendix C.
The core size, rc , as estimated by the radius of gyration is more sensitive to the integration
domain size and the vorticity cutoff level. It became necessary to include a secondary core size
estimator to account for the wide range of parameters considered in this study compared to those of

52

Bohl & Koochesfahani (2009). This alternative core size is the estimate of the gaussian core size,
rÎł , computed from the best fit (in a least squares sense) to



hωi (x, y) = ω pk e

−



 
 
x−xc 2 + y−yc 2
rÎł
rÎł

,

(2.12)

where the peak vorticity, vortex center (xc ,yc ) are also left as fit parameters. The same vorticity
cutoff level is used but all vorticity within half a wavelength is included for determining the fit.
The positive and negative vortices are processed separately yielding individual measurements
for Γ0 , ω pk , rc , rγ , a, Uc . The average value of the positive and negative vortices will be reported to
focus the discussion of these properties and their dependence on k, St, and Stte . A brief comparison
between the values obtained from the positive and negative vortices is included in Appendix F.

53

CHAPTER 3
STREAMWISE FORCE

In this chapter the phase-averaged results will be used to in conjunction with a control volume
method to compute the mean streamwise force (drag/thrust) acting on oscillating airfoils. The
chapter starts with a comparison of the rigid airfoil results to both previous computational and
experimental data to establish a baseline for this study. The effect of using the thinner airfoil
profile, a NACA 0009, rather than the NACA 0012 used in much of the literature on rigid airfoils is
investigated. After validation of the baseline rigid airfoil, the effect of flexibility on the mean flow
in the wake of the airfoilis considered. Velocity profiles of the rigid and flexible airfoil are extracted
at the same downstream location where the control volume method is applied for comparisons.
The mean thrust of the rigid and flexible airfoils is compared and scaling based on k, St and Stte is
considered.

3.1

Rigid Airfoil Results

The rigid airfoil serves as the baseline of comparison in this study. To this end, the rigid
airfoil results will be compared against similar experiments and numerical simulations. Velocity
profiles will be compared for a subset of trajectories considered in this study against the numerical
simulations of Hammer (2017), which were performed with the NACA 0009 airfoil using the same
approach employed in his extensive study of Hammer (2016) for the NACA 0012 airfoil. The
experiments of Bohl & Koochesfahani (2009), a subset of those in Bohl (2002), are also used to for
validation purposes. Lastly, the mean streamwise force is compared across a wider range of studies
than mentioned above.

3.1.1

Velocity Profiles

The control volume method utilizes the mean and fluctuating velocity profiles in the wake of the
airfoil to calculated the mean streamwise force acting on the airfoil. A starting point of comparison
54

k
4
6
8
10



∆ hui/U∞
Max RMS
0.1 0.02
0.05 0.01
0.03 0.01
0.04 0.01



∆ hvi/U∞
Max RMS
0.03
0.02
0.02
0.01
0.01 <0.005
0.01 <0.005



∆ hui0/U∞
Max RMS
0.05 0.01
0.05 0.01
0.04 0.01
0.06 0.01



∆ hvi0/U∞
Max RMS
0.06 0.01
0.04 0.01
0.04 0.01
0.02 0.01

Table 3.1: Deviations of rigid airfoil results and numerical simulations of NACA 0009. The
maximum and RMS deviation between the velocity profiles extracted at x/c = 1 for the rigid airfoil
and the numerical simulations of the NACA 0009 airfoil of Hammer (2017).

before calculating the thrust/drag, is the velocity profiles themselves. Figure 3.1 shows the velocity
profiles extracted one chord length (1c) downstream of the trailing edge of the rigid airfoil and
the simulations of Hammer (2017) which also used the NACA 0009 airfoil. The mean streamwise
velocity profile, hui, and its RMS fluctuation, hui0, the mean transverse velocity, hvi, and its RMS
fluctuation, hvi0, are all compared at reduced frequencies of k=4, 6, 8, and 10. Overall there is
excellent qualitative and quantitative agreement between the current study and the simulations; the
details of which, are discussed below.
In comparison to the CFD, the measured mean streamwise velocity, shown in Figure 3.1a, has
a maximum deviation of 0.10U∞ at k = 4, while all other reduced frequencies have a maximum
deviation of ≤0.05U∞ . The maximum deviation is less than 6% of U∞ for all of the other mean
velocity and velocity fluctuation profiles. The average deviation, quantified by the RMS of the
deviation of the entire profile is ≤0.02U∞ for not only hui, but also hvi, hui0, and hvi0. These
deviations are tabulated in Table 3.1 and indicate excellent quantitative agreement between the rigid
airfoil and the numerical simulations at the conditions available for direct comparisons.
Overall, there is very good quantitative agreement between the simulations and experiments,
the source of some of the deviations will be discussed next. The higher spatial resolution of
the simulations is evident at the higher reduced frequencies in the streamwise fluctuation profile
(shown in Figure 3.1b) near y/c = 0. The experiment’s lower spatial resolution and inherent
spatial smoothing makes the small localized peak in the profile just barely visible. The slight

55

(a)

(b)

(c)

(d)

Figure 3.1: Velocity profiles of the rigid NACA 0009 airfoil compared to the numerical simulations
of Hammer (2017). The simulations have a solid line and the experiments have filled circles as
markers with no connecting line with the color indicating the reduced frequency (k).

56

attenuation in the peaks in the RMS profiles of the experiment compared to the simulations is
attributed to both the aforementioned spatial averaging of the experiment and phase averaging
whereas the simulations have no phase averaging: the RMS profile is computed directly from the
96 instantaneous realizations available which were equally spaced in the phase of the motion.
Quantitatively, the experiments and the simulations deviate the most at a reduced frequency of
k=4. Comparing the experiment to the simulation (whose results are nearly perfectly symmetric
about the centerline) at k = 4 the experiment is seen to produce slightly asymmetric profiles. A
number of experimental uncertainties can contribute to asymmetries, including the mean angle
of attack and the asymmetry in the airfoil geometry itself. The flow itself can develop large
asymmetries, which has been studied by Dong et al. (2006), Godoy-Diana et al. (2009), Marais
et al. (2012), Calderon et al. (2014), and Wei & Zheng (2014). Most of these studies show that large
asymmetries in the flow typically occur at high reduced frequency while the experiments of Bohl &
Koochesfahani (2009) at α0 = 2◦ did observe a slight tilt in the flow field developing at k = 11.5.
From the literature, most of the asymmetries in these flows are typically observed at high reduced
frequency, so the small asymmetries observed in the current experiments compared to simulations
are most likely due to the aforementioned angle of attack uncertainty and airfoil geometry.
While the maximum values of hvi are very small for all cases (less than 2.4% of freestream),
there is a marked difference between the experiments and the simulations. The primary difference
is the simulations have no value greater than 0.6% of freestream compared to the 2.4% seen in
the experiments. However, the v-velocity profile is not used in the control volume method, so it
will not directly impact the the estimates for drag/thrust. These differences can be linked to the
simulations being 2D with a domain of over 100c in all directions compared to the AR=2.8 airfoil
with sidewalls located at y/c = Âą2.5 in the experiment. The true 3D nature of the experiments
can lead to differences in v-velocity even with the u-velocity agreeing since there will also be a
w-velocity.
Comparisons with the experiments of Bohl & Koochesfahani (2009) are slightly more difficult as
the current study did not replicate the same exact reduced frequencies of the earlier measurements

57

in addition to the difference in airfoil profile (NACA 0012 in Bohl & Koochesfahani (2009)
compared to the current NACA 0009). The velocity profiles of the current rigid airfoil at k=0
(static airfoil), 5, 6, and 12 will be compared against k=0 (static airfoil), 5.2, 5.7, and 11.5 from
Bohl & Koochesfahani (2009). While the reduced frequencies are close, the flow is highly sensitive
to k and good quantitative agreement would not be expected. Figure 3.2 shows the comparison
of hui, hvi, hui0, and hui0 for these cases. A discussion of each k considered for comparison is
discussed below.
The first observation of the differences between the current study and Bohl & Koochesfahani
(2009) for the case of the static airfoil is the smaller width of the wake deficit region than the
thicker airfoil; which can be observed in the profile of hui shown in Figure 3.2a. This an expected
outcome due to the thinner airfoil section. The fluctuation in the mean velocity (Figure 3.2b)
shows a qualitatively similar profile but is again thinner for the current experiment. The baseline
freestream level of fluctuations is lower in the current experiments by approximately 20%-30%
compared to the previous study. These differences are attributed to improvements in camera sensor
noise between the studies and while the differences are only 2-3% of U∞ , the fluctuation profiles
are used in estimating the streamwise force and hence fluctuation velocity errors would produce
errors in the thrust/drag. However, this problem is unique to the static airfoil case as all other cases
are phase averaged, which minimizes the influence of random noise on the measurements. As will
be seen, for all of the moving airfoil cases, the baseline fluctuation levels between both studies are
consistent with one another.
The next set of cases to compare are for a reduced frequency of k=5 for the current study and
k=5.2 from Bohl & Koochesfahani (2009). Considering the mean streamwise velocity profile,
both studies have a wake profile with a decreased maximum deficit compared the static case with
a consistent value across both studies. Increasing the reduced frequency to k = 5.7, Bohl &
Koochesfahani (2009) has captured the neutral wake, where the mean velocity profile is nearly U∞
across the entire transverse domain. The next case for the current study is k = 6, which has a peak
mean streamwise velocity excess of 1.12U∞ , whose excess is the same as the velocity deficit seen at

58

(a)

(b)

(c)

)
NACA 0009
)
NACA 0012

Figure 3.2: Velocity profiles of the rigid airfoil compared to the experiments of Bohl & Koochesfahani (2009). The previous experiments have a solid line and the current experiments have filled
circles as markers with no connecting line with the color indicating the reduced frequency (k).
The previous experiments (B&K 2009) used the thicker NACA 0012 profile compared to the
NACA 0009 airfoil of the current study.

59

k = 5 (0.88U∞ ), leading one to estimate, based on linear interpolation, the neutral wake condition
of the current experiments to be approximately k = 5.5, slightly lower than what is observed by
Bohl & Koochesfahani (2009) for the NACA 0012 airfoil. While Bohl & Koochesfahani (2009) did
show that the neutral wake, which corresponds to the change in the vortex array pattern switching
orientation from a von KĂĄrmĂĄn wake to a reverse von KĂĄrmĂĄn wake, is not sufficient for predicting
drag-to-thrust crossover frequency, the decrease in k at which the neutral wake occurs in the current
experiment indicates that there will be a difference in the drag-to-thrust crossover frequency. This
is likely connected to the NACA 0009 airfoil section having a lower static drag coefficient than the
NACA 0012 airfoil, which, will be considered in the next section.
At the highest k considered for comparison, k=11.5 of Bohl & Koochesfahani (2009) and
k=12 for the current rigid airfoil continues to show deviations between the two cases. With the
reduced frequencies differing by k=0.5, it is difficult to draw explicit conclusions regarding the
differences between the two experiments. The same general trend, in both hui and hvi0, is seen in
each of the experiments as k is increased. The triple peak signature in the experiments of Bohl &
Koochesfahani (2009) for hui0 for k=11.5 is not seen in the current experiments at k=12.

3.1.2

Streamwise forces

The streamwise force of the of the current experiments is computed using the control volume
method described in Section 2.5.4. The drag/thrust of the rigid airfoil will be compared to the
numerical simulations of the NACA 0009 airfoil by Hammer (2017) using both the control volume
method and the values he obtained from the direct integration of the shear stress and pressure on the
airfoil’s surface. The discussion will shift to comparisons with the wider range of cases available
in Bohl & Koochesfahani (2009) and Hammer (2016) in addition to several additional studies all
of which used the NACA 0012 profile.
The mean streamwise force coefficient for the rigid airfoil at α0 = 2◦ is plotted against k in
Figure 3.3 with the numerical simulations of the NACA 0009 airfoil from Hammer (2017) and the
numerical simulations of the NACA 0012 from Hammer (2017). The static drag coefficient of the

60

Figure 3.3: Mean Ct of the rigid airfoil and numerical simulations at α0 = 2◦ . The mean
thrust coefficient (Ct ) plotted against reduced frequency (k) for the rigid airfoil and compared to
the numerical simulations of a NACA 0009 airfoil by Hammer (2017) and the simulations of a
NACA 0012 airfoil by Hammer (2016). The control volume method is performed on the numerical
simulations at k=4, 6, 8, and 10 for the NACA 0009, indicated by “C.V.” in the legend, compared
to “direct” which is the direct integration of the pressure and shear stress.

rigid airfoil is computed to be Cd = 0.025 Âą 0.006 compared to Cd = 0.031 of the simulations. At
k ≤ 6, the current experiments are consistent with the numerical simulations of the NACA 0009
airfoil, while at higher k there is more deviation. At the reduced frequency of k=8, there is a large
difference between the control volume method and the actual Ct than what is seen at other reduced
frequencies. This shows that the control volume method may be responsible for a slight bias or shift
in the estimated crossover frequency, but more further study is required. Looking at the curve as a
whole, the crossover from drag to thrust of the experimental results appears to be shifted to a higher
k than the numerical simulations beyond what is seen to be the influence of the control volume
method. The numerical simulations (of the NACA 0009 airfoil) yield an estimate of the crossover

61

to be k=7.25 based on the actual Ct and k=8.25 based on the control volume method. The current
experiments leads to an estimate of the crossover frequency at k=9. The crossover frequency of the
simulations of the NACA 0012 airfoil is k=8.1. The NACA 0009 airfoil was expected (based on
the simulations) to crossover at a lower k compared to what is seen in the experiment.
To put in perspective, it is noted that the discrepancy between the experimental and computation
regarding the crossover frequency is within the scatter in the greater body of literature, which is
shown in Figure 3.4. In this figure the thrust coefficient, Ct is plotted against k for the numerical
simulations of Ramamurti & Sandberg (2001), and Young & Lai (2004), joining those of Hammer
(2016, 2017), in addition to the experiments of Bohl (2002), Bohl & Koochesfahani (2009), and
Mackowski & Williamson (2015). The current results of the rigid airfoil fall well within the scatter
in the data from all of these studies, and tend to agree better with the experiments of Bohl (2002)
and Bohl & Koochesfahani (2009) compared to the experiments of Mackowski & Williamson
(2015). Plotting Ct against Strouhal number, St, which is shown in Figure 3.5, shows decent
collapse in the data across different oscillation amplitudes, with the current rigid airfoil results
falling within the scatter of the results from all of these studies for all of the amplitudes in the
current experiments. The current rigid airfoil experiments estimate the drag-to-thrust crossover
frequency to be St = 0.13 − 0.15.
Additional observations can be made utilizing the additional data in the literature to form a
hypothesis regarding the differences seen in the drag-to-thrust crossover frequency between the
rigid airfoil of the current study, the previous experiments of Bohl & Koochesfahani (2009) and the
simulations of Hammer (2016, 2017). The aspect ratio of the numerical simulations is effectively
infinite, while Bohl & Koochesfahani (2009) has an AR = 4, and the current experiments has an
AR = 2.8, and the experiments of Mackowski & Williamson (2015) have an AR = 2.4. Referring
to Figure 3.4, the crossover frequency for of these studies are k=8, k=7.71, k=9, and k=10.1,
respectively for each of the studies in decreasing aspect ratio. This indicates that the aspect ratio
may be playing a significant role in shifting the drag-to-thrust crossover frequency; but this cannot
1 The crossover frequency listed in Bohl & Koochesfahani (2009) is 8.7Âą0.2, but has since been

updated to 7.7Âą0.1(but this has not yet been published).
62

Figure 3.4: Mean thrust coefficient (Ct ) plotted against reduced frequency (k) for various different
studies. The current results (Rigid) and Hammer 2017 utilize a NACA 0009 airfoil and the rest of
the studies utilize a NACA 0012 profile. Color represents a single study, marker shape represents
pitching amplitude Îą0 , and markers outlined in black are numerical simulations while non-outlined
markers are experimental studies. Legend entries have been shorted for space: B&K - Bohl &
Koochesfahani (2009), R&S - Ramamurti & Sandberg (2001), Y&L - Young & Lai (2004), M&W
- Mackowski & Williamson (2015).

be concluded with the data available as the current rigid airfoil is a NACA 0009 compared to the
NACA 0012 of the rest of the studies.
In summary, the current rigid airfoil results compare favorably with previous experiments and
computations. The NACA 0009 airfoil does differ slightly in the wake velocity profiles compared
to the NACA 0012 airfoil, but it is seen as a slight shift in k attributed to the differences in the
drag of the static airfoil. The simulations of Hammer (2017) indicate the NACA 0009 airfoil has a
lower crossover frequency compared to the NACA 0012 (k=7.2 versus k=8.2). The control volume
method was shown, again using the simulations of Hammer (2017), to shift the crossover frequency
for the NACA 0009 airfoil from k=7.2 to k = 8.2 for α0 = 2◦ . The rigid NACA 0009 airfoil of the
current study has a crossover frequency of k=8.9, slightly higher than the simulations, but within

63

Figure 3.5: Mean thrust coefficient (Ct ) plotted against Strouhal number (St) for various different
studies. The current results (Rigid) and Hammer 2017 utilize a NACA 0009 airfoil and the rest of
the studies utilize a NACA 0012 profile. Color represents a single study, marker shape represents
pitching amplitude Îą0 , and markers outlined in black are numerical simulations while non-outlined
markers are experimental studies. Legend entries have been shortened for space: B&K - Bohl &
Koochesfahani (2009), R&S - Ramamurti & Sandberg (2001), Y&L - Young & Lai (2004), M&W
- Mackowski & Williamson (2015).

the uncertainty of the current experiment. At the higher oscillation amplitudes of α0 = 4◦ and
α0 = 6◦ the current results for mean Ct agree with the computations of Hammer (2016) and Bohl
(2002) to within the experimental uncertainty. With the current baseline measurements validated,
the effect of flexibility on the the mean streamwise force will be explored next.

3.2

Effect of Flexibility

The effect of flexibility on the mean velocity profiles and the mean streamwise force is examined
here. The first subsection will examine the velocity profiles extracted at x/c = 1; which are the
same velocity profiles used in computing the streamwise force using the control volume method.
The velocity profiles are compared for the same prescribed motions between the rigid and flexible

64

airfoils, which highlights the direct effect of flexibility on the velocity fields. The following
subsection examines the effect of chordwise flexibility on the streamwise force in terms of the
prescribed motion, and the actual trailing edge amplitude. The effect of chordwise flexibility on the
drag/thrust of a pitching airfoil will be shown to collapse when plotted against the Strouhal number
based on the actual trailing edge amplitude (Stte ).

3.2.1

Velocity Profiles

Figure 3.6 shows the velocity profiles of hui, and hvi while Figure 3.7 plots the velocity profiles of
hui0 and hvi0 for the oscillation amplitude of α0 = 2◦ and reduced frequencies of k=2, 5, 6, and 10.
The discussion below outlines the effect of flexibility on the mean and fluctuating velocity profiles
as a function of the reduced frequency.
At the lowest reduced frequency shown, k=2, the rigid and flexible airfoil have nearly the same
mean and fluctuating velocity profiles. At this reduced frequency, the trailing edge displacement of
the flexible airfoil is only 5.7Âą2.2% (0.36 mm) larger than the rigid airfoil; this small change in h p
does not modify the flow significantly. However, at a k=5, h p of the flexible airfoil is 16.3Âą3.2%
higher than the rigid airfoil and the change in the flowfield is significant. The mean streamwise
velocity for the rigid airfoil has a wake profile while the flexible airfoil has switched to a jet-like
profile, closer to the jet-like profile of k=6 for the rigid airfoil. The fluctuating velocity profiles
convey a similar message with both hui0, and hvi0, of the flexible airfoil at k=5 exhibiting velocity
profiles closer to the rigid airfoil case at k=6.
The highest reduced frequency plotted for comparison, k=10, shows a much more significant
difference between the airfoils than at the lower reduced frequencies previously discussed. The hui
profile for both rigid and flexible airfoils have a strong jet-like profile, but the flexible’s centerline
velocity magnitude is 1.95U∞ compared to 1.6U∞ of the rigid airfoil. The width (in y/c) of
the velocity excess is much larger for the flexible airfoil as well. The increase in width of the
velocity profiles is also observed in the fluctuating velocity profiles. The notable increase in width
is consistent with h p of the flexible airfoil being 48Âą5.7% higher than the rigid airfoil; making the

65

(a)

(b)

Figure 3.6: Selected mean velocity profiles at α0 = 2◦ . Mean velocity profiles extracted at x/c=1
of both the rigid and flexible airfoils at an oscillating amplitude of α0 = 2◦ .

66

(a)

(b)

Figure 3.7: Selected fluctuating streamwise velocity profiles at α0 = 2◦ . Fluctuating streamwise
velocity profiles extracted at x/c=1 of both the rigid and flexible airfoils at an oscillating amplitude
of α0 = 2◦ .

67

effective α0 nearly 3◦ compared to the prescribed amplitude of α0 = 2◦ .
Similar effects as those described for α0 = 2◦ are seen at the higher pitch amplitudes of α0 = 4◦
and α0 = 6◦ . For brevity, only the mean streamwise velocities of both of the higher amplitudes are
compared in Figure 3.8 for the rigid and flexible airfoils. While it appears, by comparing α0 = 4◦
and α0 = 6◦ to α0 = 2◦ , shown in Figure 3.6a, that the effect of flexibility is reduced at the higher
amplitudes, this is not the case. Comparing the same reduced frequency from different oscillation
amplitudes yields qualitatively similar changes in the velocity profiles between the rigid and flexible
airfoils. This is consistent with the relative change in h p (see Figure 2.18) being nearly the same
for different Îą0 as long as k<6.

3.2.2

Streamwise force

The mean streamwise force is computed from the control volume method, described previously in
Section 2.5.4, using the velocity profiles extracted 1c downstream of the airfoil trailing edge for
both the rigid and flexible airfoils. The mean thrust, Ct , is plotted against reduced frequency in
Figure 3.9. Not all cases studied are included in the figure to focus the discussion on the behavior
near the drag-to-thrust crossover; all cases can be located in Appendix D.
The static airfoil measurements were repeated twice at each field of view (FOV), one before all
of the motion cases, and one at the end. The measurements on the flexible airfoil were repeated
twice with a separate set of FOVs to yield a second, independent set of measurements for the
mean streamwise force. Therefore, there are four measurements of the static drag coefficient for
the flexible airfoil. The drag coefficient of the stationary rigid airfoil is Cd =0.026Âą0.006 and
Cd =0.039Âą0.008 for the stationary flexible airfoil. While the flexible airfoil does have, on average,
more drag than the rigid (reasonable based on the extra drooping trailing edge of the flexible airfoil)
the differences are within the experimental uncertainty.
The flexible airfoil is shown to produce more thrust (or less drag) under nearly all reduced
frequencies and amplitudes considered in the current study with a few exceptions at k ≤ 2. This
shifts the drag-to-thrust crossover frequency to a lower frequency for all amplitudes. The highest

68

(a) α0 = 4◦

(b) α0 = 6◦

Figure 3.8: Selected mean streamwise velocity profiles at α0 = 4◦ and α0 = 6◦ . Mean streamwise
velocity profiles extracted at x/c=1 for selected reduced frequencies of both the rigid and flexible
airfoils at an oscillating amplitude of α0 = 4◦ (a) and α0 = 6◦ (b).

69

Figure 3.9: Mean thrust plotted against reduced frequency for rigid and flexible airfoils. The mean
thrust computed from the control volume method using velocity profiles extracted at 1c downstream
of the trailing edge. Error bars indicate a 95% confidence interval.

amplitude’s (α0 = 6◦ ) crossover frequency is shifted the least by only ∆k ≈ 0.1, while α0 = 4◦
is shifted by ∆k ≈ 0.5 and α0 = 2◦ is shifted the most by ∆k = 2.5. While the increase in thrust
for a flexible airfoil compared to its rigid counterpart is well documented in the literature (Dewey
et al., 2013; Gursul et al., 2014), this study shows the influence of chordwise flexibility in shifting
the drag-to-thrust crossover to a lower frequency for the flexible relative to the rigid airfoil.
Figure 3.10 shows the mean thrust coefficient plotted against Strouhal number (St). The rigid
airfoil results collapse well (within the experimental uncertainty) for all Îą0 as shown previously
in Section 3.1.2 but there is considerable scatter for the flexible airfoil results. Figure 3.11
shows Ct plotted against the Strouhal number based on the actual trailing edge displacement,
i.e. Stte = f hP /U∞ , and much better collapse for both rigid and flexible airfoils across all α0 is
observed.
Triantafyllou et al. (1991); Triantafyllou (1993) first proposed the use of Strouhal number as the
parameter to describe the wake dynamics of flapping airfoils. Triantafyllou (1993), more explicitly

70

Figure 3.10: Mean thrust plotted against Strouhal number for rigid and flexible airfoils. The mean
thrust computed from the control volume method using velocity profiles extracted at 1c downstream
of the trailing edge. Error bars indicate a 95% confidence interval.

stated that St should be taken as the maximum excursion of the trailing edge (i.e. h p in the current
measurements) but the direct application of this definition to chordwise flexibility has not been seen
in the literature. Thiria & Godoy-Diana (2010) still used St based only on the prescribed motion, or
St based on other fixed length scales such as c in Heathcote & Gursul (2007), or airfoil thickness in
Godoy-Diana et al. (2008, 2009) and Marais et al. (2012).Dewey et al. (2013) defines the Strouhal
number in this way (based on the actual h p ) but present their data based with k, pointing out that
Stte is an output of the system and therefore, not known a priori. The current study reaffirms
that it is necessary to use the actual trailing edge amplitude to properly characterize these types of
flows. Recently, Monnier et al. (2015) showed Stte is better at describing the characteristics of the
flowfield (in terms of various vortex array properties) compared to St for pitching flexible airfoils of
the same teardrop based design of Heathcote & Gursul (2007) at a lower Reynolds than the current
study of Rec = 2 × 103 .
As previously discussed, the mean thrust coefficient, Ct , when plotted against Stte is shown

71

Figure 3.11: Mean thrust coefficient plotted against Strouhal number based on the actual trailing
edge displacement for rigid and flexible airfoils. The mean thrust computed from the control
volume method using velocity profiles extracted at 1c downstream of the trailing edge. Error bars
indicate a 95% confidence interval.

to collapse the data in Figure 3.11. The scatter in the crossover Strouhal number (Stte ) is well
within the scatter shown previously in Figure 3.5 for several different pitching airfoil studies.
Considering botht the rigid and flexible airfoils for amplitudes studied, the crossover frequency
is Stte =0.13Âą0.02. This collapse based on Stte shows that, to first order, the effect of chordwise
flexibility can be considered as an amplitude effect, i.e. the peak to peak excursion of the trailing
edge is sufficient to capture the effect of chordwise flexibility on the mean thrust. The implications
of this is the mean thrust is not effected by the shape change (camber) of the airfoil in this
study. However, the current study’s chordwise flexibility is not “very” flexible compared to many
other studies, and only the first bending mode is observed, so the implication of mean thrust not
being effected by the shape change of the airfoil may not apply to all other levels of flexibility
exhibiting more complicated structural dynamics. Studies investigating the shape change of thin
(typically flat plate or membrane wings) and their effect on forces and efficiency can be found in

72

the analytical models of Katz & Weihs (1978, 1979) and Moore (2014, 2015) and the experimental
based discussions in Ramananarivo et al. (2011) and Dewey et al. (2013).

3.3

Summary

In this chapter the effect of chordwise flexibility on the streamwise force was investigated.
The measured thrust based on the control volume method was compared for the baseline rigid
NACA 0009 airfoil against several different numerical and experimental studies. The effect of
the NACA 0009 airfoil was shown, via a comparison of 2D numerical simulations, to lower the
crossover frequency relative to the typically studied NACA 0012 airfoil, but the current data were
seen to have a crossover frequency higher than the previous experiments of Bohl & Koochesfahani
(2009) at α0 = 2◦ . These differences are attributed to aspect ratio effects, but requires further
study. The current experiments at higher amplitudes (α0 = 4◦ and α0 = 6◦ ) were observed to agree
well with the available data in the literature in terms of mean Ct . Chordwise flexibility was seen
to increase the thrust of the oscillating airfoil under nearly all conditions. While that observation
is not a new finding, the shift in the crossover frequency due to chordwise flexibility has not been
reported previously in the literature. This study finds that the crossover frequency, at a prescribed
oscillation amplitude, decreases due to chordwise flexibility.

73

CHAPTER 4
WAKE VORTEX PROPERTIES

In this chapter the phase-averaged results will be used to explore the effect of chordwise flexibility
on the vortical patterns in the wake of oscillating airfoils. The chapter focuses primarily on motion
trajectories which produce a single pair of vortices, one of each sign, in the wake per oscillation
cycle of the airfoil. These wakes are typically referred to as either the classical von KĂĄrmĂĄn vortex
street or reverse von KĂĄrmĂĄn vortex street, depending on their arrangement. The chapter begins with
a comparison of the rigid airfoil results to both computational and experimental data to establish
a baseline for this study. The effect of using the thinner airfoil profile, a NACA 0009, rather than
the NACA 0012 used in the comprehensive studies of Bohl (2002), Bohl & Koochesfahani (2009),
and Hammer (2016) is investigated. The chapter continues with a detailed comparison between
the rigid and flexible airfoils in terms of the vortex strength (peak vorticity, ω pk , and circulation,
Γ0 ), size (radius of gyration, rc ), and arrangement (wavelength, a, and transverse spacing, b). The
definitions of these terms and how they are computed for the current experiments was described in
Section 2.5.5. The appropriate scaling of these wake vortex parameters are discussed in the context
of flexibility. The data and scaling parameters are compared, and contrasted, to those seen in the
literature when appropriate.

4.1
4.1.1

Validation of Rigid Airfoil Results
Comparison with 2D Simulations of the NACA 0009 airfoil

The current study utilizes a NACA 0009 airfoil profile while in the recent studies of Bohl (2002),
Bohl & Koochesfahani (2009), and Hammer (2016), who also investigated the wake vortex properties in the wake of an sinusoidally pitching airfoil, used a NACA 0012 airfoil. In the previous chapter
the effect of the thinner airfoil profile on the mean flowfield was described, with the primary effect
seen to be the lower static drag coefficient, which causes, for the oscillating airfoil, the crossover

74

frequency from drag to thrust to occur at a lower reduced frequency than the NACA 0012. The
simulations of Hammer (2017), which were conducted as described in Hammer (2016), but on the
NACA 0009 airfoil are again used in this section to compare the vorticity field both qualitatively
and quantitatively in terms of the wake vortex parameters.
The phase-averaged vorticity of the experiment is compared (at the same phase of the motion)
to the instantaneous vorticity of the simulation in Figure 4.1 at reduced frequencies of k=4, 6, 8,
and 10 at an amplitude of α0 = 2◦ . Qualitatively, the experiments and computations produce very
similar flow structures at each reduced frequency.
At k=4, the shear layer in the experiment is seen to be rolling up into two positive and two
negative vortices by x/c=0.5 while the one in the computation computation doesn’t roll up as
quickly. The braid vorticity is observed in the experiments to be thicker, or more diffuse. The
experiment is phase-averaged over 1/64th of the cycle in contrast to the computation which is an
instantaneous result. The phase-averaging process will result in a small amount of spatial averaging
related to how far the flow structures convect during one phase bin. Cycle-to-cycle variations in
the spatial location of the flow structures also contribute to a certain amount of spatial averaging.
At k=4, the flow structures (vortices) convect 0.0125c in during one phase bin (1/64th of the
cycle), just over one experimental grid point. The flow structures convect faster at higher k, which,
counter-intuitively, means the effective spatial averaging due to the phase-averaging process will be
negligible at higher reduced frequencies and this explains why some of the braid vorticity appears
to be more diffuse in the experiment than in the simulation.
At k = 8 and k = 10 the vortex array pattern is seen to tilt downwards in the experiment but not
in the computation. A slight tilt in the wake is commonly observed in experimental studies; it can
be observed in the classic study by Koochesfahani (1989), and Bohl & Koochesfahani (2009) also
observe a slight tilt in the wake pattern at k=11.5. In general however, the wake of a sinusoidally
pitching airfoil is expected to remain mostly symmetrical at these reduced frequencies as seen in
the studies of Koochesfahani (1989), Bohl & Koochesfahani (2009), Godoy-Diana et al. (2008),
Godoy-Diana et al. (2009), and others. The slight asymmetry in the baseline results is contributed to

75

Computation

Experiment

Figure 4.1: Comparison of experiment and simulations of phase-averaged vorticity fields. The
phase-averaged vorticity at the same phase of oscillation: computations (left) and experiment
(right). The reduced frequency increases from top to bottom, with the reduced frequency listed
as the title of each plot. The experiment and simulation have identical colormaps for the same
reduced frequency, but the colormap limits increase as k increases. The phase shown corresponds
to α = 0◦ while the airfoil is pitching down (trailing edge moving upwards).

76

either by a slight non-zero mean angle of attack (the current experiment determined αm = 0◦ ±0.1◦ )
and/or the asymmetry in airfoil profile discussed in Section 2.3.3.
Another observation between the experiments and the computations is the reduction in the rates
at which the peak vorticity decreases as the vortices convect downstream. In the study by Bohl
& Koochesfahani (2009), the streamwise variation in the peak spanwise vorticity is attributed to
the combined effects of vortex viscous diffusion and vortex stretching due to spanwise velocity
gradients, with the latter being the result of the axial velocity. Since the computations are 2D, they
do not replicate this feature of the real flows. Further analysis of this effect is beyond the scope of
the current study.
The wake vortex properties of the computations were analyzed with the same post-processing
tools used in the current study to ensure a 1:1 comparison in the computed parameters. The ω pk ,
Γ0 , rc , and b/a were computed at x/c=0.5 for both the rigid airfoil and the simulations of the
NACA 0009 airfoil of Hammer (2017). The results are shown in Figure 4.2 for the average value of
the positive and negative vortices (except in the case of b, since it is the transverse distance between
the positive and negative sign vortex). A discussion of the differences between the positive and
negative vortex properties is relegated to Appendix F, since in general, the differences are small.
Overall, quantitative comparison shows that the experiments and computations agree favorably as
discussed in detail below.
At a reduced frequency of k=6, 8, and 10 there is less than 3% difference between the experiments
and the computations for Γ0 and a/c (not shown in Figure), with less than 9% difference for rc . The
transverse spacing normalized by the wavelength (b/a), shown in Figure 4.2(d), agrees with the
computations within the experimental uncertainty. The peak vorticity in the experiment however,
is ∟10-20% less than the computations. The accuracy in measuring vorticity is not expected to
contribute to this discrepancy; the remapping procedure and subsequent finite difference method
used for calculating the vorticity of the experimental data is expected, based on the analysis of
Cohn & Koochesfahani (2000), to be less than 2% of the peak vorticity and therefore, is also not
significant enough to explain the differences. A brief study in the effect of spatial resolution on

77

(a)

(b)

(c)

(d)

Figure 4.2: Vortex Properties of Rigid Airfoil and Simulations of NACA 0009. The peak vorticity
(a), circulation (b), radius of gyration (c), and the transverse spacing normalized by the wavelength
(d) are the average value of the (+) and (-) sign vortex. Both the experiments (Rigid) and the
numerical simulations (Hammer 2017) use the NACA 0009 airfoil profile. For the experimental
dataset invisible error bars (e.g. k < 12 for Γ0 /cU∞ ) indicate the error bar is less than or equal to
the symbol size.

78

the extracted vortex parameters was conducted with the simulation data. The resolution of the
simulation was reduced from 0.0025c to 0.01c to match the current experiments and the vorticity
was calculated with the same method as used in the experiments (an explicit fourth-order-central
finite difference method). Additional details of this analysis can be found in Appendix G. When
the resolution is reduced to that of the experiment, the peak vorticity decreases by 1.5% to 8%.
The remainder of the difference in peak vorticity is attributed to the previously discussed vortex
stretching mechanism that the computations do not capture. Koochesfahani (1989) observed that
the magnitude of the axial flow increases with increasing k, which is consistent with the growing
discrepancy between the experiments and computations as k increases.

4.1.2

Effect of NACA 0009 on vortex properties

In preparation of validating the current experiments with the experiments of Bohl & Koochesfahani
(2009), the effect of the thinner airfoil on the vortex properties is investigated with the numerical
simulations. The computations of Hammer (2016) and Hammer (2017) are used for this purpose
as the he used the NACA 0012 and NACA 0009 airfoils in each of the respective data sets. The
vortex properties of ω pk , Γ0 , rc , and b/a are shown in Figure 4.3 for the two computations and the
rigid airfoil of the current study for reference.
The following discussion is based on the comparison between the NACA 0009 and the
NACA 0012 airfoil of the numerical simulations, and is limited to α0 = 2◦ . The peak vorticity at k=6, 8, and 10 is higher by 4-8% for the thinner airfoil than the thicker airfoil. The circulation
is 3.5% less for k=6, but 1.9% and 2.6% greater for k=8 and k=10, respectively, for the thinner
airfoil. The core size is similar to circulation in that the thinner airfoil has a smaller vortex at k=6
(by 1.6%) but is larger at k=8 and k=10 (by 2.7% and 2.1% respectively). The b/a spacing tells a
more important (but still small) difference between the two airfoils: the switch in the vortex pattern
from a von KĂĄrmĂĄn to a reverse von KĂĄrmĂĄn wake structure occurs at a lower k for the thinner
airfoil. This is consistent with thinner airfoil having a lower drag-to-thrust crossover frequency
which was seen (for the numerical simulations) in Chapter 3.

79

(a)

(b)

(c)

(d)

Figure 4.3: Vortex Properties of the rigid airfoil compared to numerical simulations of the NACA
0009 and NACA 0012 airfoils. The peak vorticity (a), circulation (b), radius of gyration (c), and
the transverse spacing normalized by the wavelength (d) are the average value of the (+) and (-)
sign vortex. Both the experiments (Rigid) and the numerical simulations of Hammer 2017 use the
NACA 0009 airfoil profile while the numerical simulations of Hammer 2016 utilize the NACA
0012 profile. For the experimental dataset invisible error bars (e.g. k < 12 for Γ0 /cU∞ ) indicate
the error bar is less than or equal to the symbol size.

80

There is conflicting observations of the current experiment’s rigid NACA 0009 airfoil regarding the reduced frequency at which the vortex pattern switches and the drag-to-thrust crossover
frequency for α0 = 2◦ . In Chapter 3 the drag-to-thrust crossover frequency for the rigid airfoil was shown to be higher than the numerical simulations of both a NACA 0009 airfoil and a
NACA 0012 airfoil, and the experiments of Bohl & Koochesfahani (2009) of a NACA 0012 airfoil.
The switch in vortex pattern of the current experiment’s NACA 0009 airfoil however, is consistent
with the simulations of the NACA 0009, with the switch in the vortex pattern occurring at lower
k compared to the thicker NACA 0012 airfoil. Bohl & Koochesfahani (2009) pointed out that the
switch in the vortex pattern is not sufficient in determining the drag-to-thrust crossover frequency.
This study confirms this observation and notes the difference between drag-to-thrust crossover and
vortex pattern switch is also not strictly associated with the static drag coefficient of the airfoil.
The drag coefficient of the rigid airfoil was measured to be less than the drag coefficient of the
experiment of Bohl & Koochesfahani (2009) on the NACA 0012 airfoil (which is consistent with
basic knowledge, see Abbot & von Doenhoff (1959), of the effect of airfoil thickness on drag). If
the difference between in the frequency at which the vortex pattern switches and the drag-to-thrust
crossover frequency was strictly a function of the static airfoil’s drag, this drag-to-thrust crossover
frequency should have been lower than what Bohl & Koochesfahani (2009) measured, but, this is
was not observed.
The larger difference between the frequencies at which the vortex pattern switches and drag-tothrust occurs must be related to the few remaining differences between the current experiment, the
experiments of Bohl & Koochesfahani (2009) and the numerical simulations. It is hypothesized
that the likely reason is the aspect ratio, with the current study being AR=2.8 compared to AR=4
of Bohl & Koochesfahani (2009) and the infinite AR of the simulations of Hammer (2016, 2017).
Unfortunately, the limitations imposed on the current study’s primary goal of studying flexibility
prohibited the aspect ratio to be matched between the current experiments and those of Bohl &
Koochesfahani (2009). Further study of this effect is recommended.

81

4.1.3

Comparison with Experiments and 2D Simulations of the NACA 0012

In Section 4.1.1 the instantaneous vorticity field of numerical simulations on the NACA 0009 airfoil
were shown to agree well with those of the rigid airfoil. In the previous section the effect of the
thickness of the oscillating airfoil, the NACA 0009 compared to a NACA 0012, was seen to shift
the b/a curve to a slightly lower frequency, but otherwise, there was little impact on the rest of the
quantified vortex properties. Moving forward, only the computational results of the NACA 0012
from Hammer (2016) will be used rather than including the results from the NACA 0009 to reduce
clutter in the plots since there is little difference between the vortex properties of the two airfoils
extracted from the numerical simulations. In this section, the experimental study reported in
Bohl (2002) and Bohl & Koochesfahani (2009) will be used for comparison in addition to the
aforementioned computational study. The vortex properties are computed in the same fashion and
at the same downstream location of x/c = 0.5 as performed in Bohl & Koochesfahani (2009).
Overall, quantitative comparisons between the current experiment and the previous experiment of
Bohl & Koochesfahani (2009) and Bohl (2002) are in good agreement and will be discussed below.
Figure 4.4 shows the values of ω pk , Γ0 , rc , and b/a for three different pitching amplitudes,
α0 = 2◦, 4◦ , and 6◦ . The peak vorticity (ω pk ), Figure 4.4(a), is consistent with, and agrees
well overall with the other experiment and the computation for all oscillation amplitudes. At the
lowest amplitude of α0 = 2◦ , the peak vorticity agrees more with the other experiment than the
computation, which supports the previous explanation of the differences between the computation
and the current experiments being primarily due to the 2D computations not being able to capture
the vortex stretching phenomena and the finite aspect ratio of the experiments. There is good
agreement, albeit with some scatter, for both α0 = 4◦ and α0 = 6◦ in terms of the peak vorticity.
Circulation, core size, and transverse spacing normalized by wavelength, Figure 4.4(b), (c), and
(d), respectively, all agree very well across all of the studies. The overall quantitative agreement
between the baseline measurements of the rigid airfoil and the literature gives a strong indication
that the current study is faithfully reproducing, and measuring, the desired flow conditions. The
rest of this chapter is focused on the influence of flexibility on these vortex properties.

82

(a)

(b)

(c)

(d)

Figure 4.4: Comparisons of vortex properties of the rigid airfoil and other experiments and
simulations of a NACA 0012 airfoil. The peak vorticity (a), circulation (b), radius of gyration (c),
and the transverse spacing normalized by wavelength (d). are the average value of the (+) and
(-) sign vortex. The current experiments (Rigid) use a NACA 0009 airfoil profile while the other
studies, Bohl (2002), Bohl & Koochesfahani (2009), and Hammer (2016), all utilize a NACA 0012
airfoil profile. For the experimental dataset, invisible error bars indicate the error bar is less than
or equal to the symbol size.
83

4.2

Effect of flexibility on vorticity field

The phase-averaged vorticity field is shown for α0 = 2◦ in Figure 4.5 for several different
reduced frequencies (2 ≤ k ≤ 10). The differences between the rigid and flexible flowfields are
subtle, but will be described next. The peak vorticity of the flexible airfoil appears to be higher
than the corresponding vortex of the rigid case (for all k). The transverse vortex spacing (b) is
notably higher at k = 8 and 10 for the flexible airfoil compared to the rigid airfoil. The vortex array
at k=10 also appears to tilt slightly downwards in the rigid case (as mentioned previously) and
slightly upwards in the flexible case. This difference is attributed to the asymmetries of the airfoil
profiles (discussed in Section 2.3.3) and when the airfoils were installed in the test section the slight
curvatures of their trailing edges’ were oriented in opposite directions. The current study does not
have enough information to confirm this is the only issue that may be contributing to the different
directions of the the wake deflection. However, the slight tilt of the vortex array is not significant
near x/c=0.5 where the vortex properties will be extracted, so the difference in the slight upwards
vs. slight downwards tilt should be negligible regarding the vortex properties under consideration.
Recalling the relative change in trailing edge amplitude is greatest at α0 = 2◦ , the changes in the
flowfield are expected to be even more subtle at α0 = 4◦ and α0 = 6◦ , and therefore, comparisons at
other amplitudes is relegated to Appendix E. The next section will make quantitative comparisons
of the vortex properties to illicit the effect of flexibility on the vorticity field.

84

Rigid

Flexible

Figure 4.5: Phase-averaged vorticity field of rigid (left) and flexible (right) airfoils. The reduced
frequency increases from top to bottom, with the reduced frequency listed as the title of each plot in
addition to the Strouhal number (St) and the Strouhal number based on the measured trailing edge
displacement (Stte ). The rigid and flexible plots have identical colormaps for the same reduced
frequency, but the colormap limits increase as k increases. The phase shown corresponds to α = 0◦
while the airfoil is pitching down (trailing edge moving upwards).

85

4.3

Wake vortex properties of rigid and flexible airfoils

In this section the effect of chordwise flexibility on the wake vortex properties will be examined.
The effect on peak vorticity (ω pk ), circulation (Γ0 ) core size (rc , rγ ), convective speed (Uc ), wavelength (a), and transverse spacing (b) will be considered. The appropriate non-dimensionalization
for each term will be evaluated and the appropriate scaling with motion trajectory parameters will
investigated. Specifically, the current study will demonstrate that the influence of flexibility can be
captured as an effective amplitude effect by using the actual trailing edge amplitude, (h p ).

4.3.1

Peak Vorticity

The peak vorticity nondimensionalized by c/U∞ is plotted against k for the rigid and flexible airfoils
in Figure 4.6a. There are three distinct groupings visible for each Îą0 , with the flexible airfoil data
at α0 = 2◦ starting to deviate from the rigid airfoil more significantly as k increases. When plotted
against St there remains three separate groupings for each amplitude (Figure 4.6b), and plotting
against Stte makes the data collapse look even worse (Figure 4.6c). The effect of flexibility on the
peak vorticity is not entirely clear; for all k at α0 = 2◦ the flexible airfoil has a higher peak vorticity,
but at the higher pitching amplitudes, this is not always true. When plotted against Stte , it is clear
that for all cases, the peak vorticity of the flexible airfoil is less than its rigid counterpart. Further
investigation in scaling ω pk is warranted to yield a consistent message and to at least, collapse the
data with respect to Îą0 .
In the literature there has been little success in scaling the peak vorticity. Akkala et al. (2015)
scaled vorticity by the frequency (ω/ f ) and observed better scaling than with the convective time
(ωc/U∞ ) for the peak vorticity values. However, they reached this conclusion based on a visual
assessment of the vorticity field, not by quantifying the actual peak vorticity for different cases.
This scaling (ωc/U∞ ) was tested with the current measurements, but no improvement was found.
The alternative scalings of ω pk h p /U∞ , ω pk h p /(c f ) were also evaluated but did not yield an overall
collapse in the data for different oscillation amplitudes. It should be noted however, that the scaling

86

(a)

(b)

(c)

Figure 4.6: Non-dimensional peak vorticity for different oscillation amplitudes of the rigid and
flexible. Non-dimensional peak vorticity plotted against k (a), St (b), and Stte (c).

87

of ω pk h p /U∞ plotted against Stte did a decent job collapsing each oscillation amplitude between
the rigid and flexible airfoils, as it is shown in Figure 4.7, but doesn’t collapse the data for different
oscillation amplitudes.

Figure 4.7: Scaling of the peak vorticity by h p /U∞ . The peak vorticity of both the rigid and
flexible airfoils are scaled as ω pk h p /U∞ and plotted against Stte .

The peak vorticity non-dimensionalized as ω pk c/U∞ is plotted against kSt, and kStte in Figure
4.8. Plotting against kSt and kStte show excellent collapse of both rigid and flexible airfoils and
for all oscillation amplitudes. While plotting against kSt appears to collapse the data better for the
entire range, plotting against kStte is better at collapsing the data for kStte < 1. This was determined
by comparing the norm of the residuals from the quadratic fit to the data, which is shown in the
Figures as well.
Wu & Wu (1993) extended the analysis of Lighthill (1963) to describe the vorticity production
at a moving wall. The analysis of Wu & Wu (1993), once sufficiently simplified for a 2D,
incompressible flow, and for pitching motion, the vorticity production becomes a function of the
wall acceleration and the pressure gradient. The term kStte is proportional to the maximum
trailing edge acceleration of the pitching airfoils. From a dimensional analysis of flapping wing

88

aerodynamics, Kang et al. (2011) also arrive at using kSt as their non-dimensional term for
acceleration. Further analysis of this scaling term is warranted, but beyond the scope of the current
investigation.
(b)

(a)

Figure 4.8: Peak Vorticity scaling with trailing edge acceleration. The non-dimensional peak
vorticity for different oscillation amplitudes of the rigid and flexible airfoils are plotted against kSt
(a) and kStte (b). The dashed lines shown are quadratic fits (in a least squares sense) to all of the
cases.

89

4.3.2

Circulation

In Sections 4.1.1 and 4.1.3 the rigid airfoil measurements of the current study was shown to agree
well across different studies for Γ0 /cU∞ at each oscillation amplitudes of α0 = 2◦ , α0 = 4◦ ,
and α0 = 6◦ . There was not however, good collapse across α0 when plotted against the reduced
frequency of k. Triantafyllou et al. (1991) proposed that the wake dynamics of oscillating airfoils is
dominated by the Strouhal number, St. Ramamurti & Sandberg (2001), as discussed in the previous
chapter, showed that St collapsed CT for α0 = 2◦ and α0 = 4◦ . The dependency on Strouhal number
for wake vortex properties are discussed in the yet-to-be published follow up analysis from Bohl
& Koochesfahani (2009), and the published works of Monnier et al. (2015) and Hammer (2016).
These works show St to be effective in collapsing Γ0 for different oscillation amplitudes. Monnier
et al. (2013, 2015) showed Stte was better than k and St for collapsing circulation for chordwise
flexible airfoils, albeit with a significant amount of scatter in the data.
Consistent with the previous discussion, there is no collapse across α0 and between the rigid/flexible airfoils for Γ0 /cU∞ when plotted against k, shown in Figure 4.9a. Better collapse is seen in
Figure 4.9b where Γ0 /cU∞ is plotted against St. There remains however, a systematic deviation in
the flexible cases as St increases (when flexibility becomes more important). Figure 4.9c shows a
complete collapse in Γ0 /cU∞ when plotted against Stte which contrasts the result of Monnier et al.
(2015) as they observed considerable scatter. Their study looked at a wide range of flexibilities
including many cases where there was “too much flexibility”, which may have contributed to the
scatter. The scatter can also be explained, in party, by their large uncertainty in determining the
trailing edge displacement, and therefore, a large uncertainty in Stte .
Monnier et al. (2015) found Γ0 /a to be a parameter that characterized the mean velocity excess
from their extended analysis of the vortex array model of Naguib et al. (2011). In an effort to
show better collapse in Γ0 , they considered scaling the circulation as Γ0 /aU∞ and observed slightly
better collapse in their data (but there was still significant scatter). This scaling was tested with the
current measurements, shown in Figure 4.9d, but no such collapse in the data was found.

90

(a)

(b)

(c)

(d)

Figure 4.9: Nondimensional circulation plotted against k and St, Stte . Circulation in the nondimensional form of Γ0 /cU∞ (a-c), and circulation nondimensionalized as Γ0 /aU∞ plotted against
Stte (d).

91

Buchholz et al. (2011) non-dimensionalized circulation as
∗ =
ΓPG

Γ
f h2p

,

(4.1)

and combined it with the similarity parameter for the pressure coefficient they developed for small
aspect ratio pitching panels to arrive at a scaling of circulation to be
∗
ΓGS

=

∗
ΓPG




hp
1+β
.
S

(4.2)

They found this scaling is effective at collapsing their low aspect ratio cases using β = 2, while
their higher aspect ratio (AR = 2.4) collapses, within their experimental uncertainty, for β = 7.
This scaling is investigated, with a β = 7, with the current data but does not collapse the data as
shown in Figure 4.10 for different pitching amplitudes or between rigid/flexible airfoils. In the
experiments of Buchholz et al. (2011), they saw not only collapse in the data for different Îą0 but
∗ was nearly a constant value. Their efforts were focused on scaling circulation for pitching
the ΓGS

3D panels, a large departure from the current experiment that strives for two dimensionality. They
also compute circulation differently seeking to capture all of the shed vorticity from one half of
a cycle compared to the current study which focuses on isolated vortices. The scaling parameter
does however show a reduced variation in circulation as a function as Stte compared to Γ0 /cU∞ .

4.3.3

Vortex core size

The size of the vortices is determined by calculating the radius of gyration of the vortex when the
ω pk is located at x/c=0.5. The core sizes of the rigid and flexible airfoils are plotted against k and
St in Figure 4.11 and against Stte in Figure 4.12a. For pitching amplitude of α0 = 2◦ and α0 = 4◦
the flexible airfoil is seen to increase the size of the vortex core size for k>6. At lower k there isn’t
any significant difference between rigid and flexible airfoils (except for the one outlier at k=5 where
the core size of the flexible airfoil is 10% smaller than the rigid airfoil). The core size collapses
only slightly better when plotted against Stte for each Îą0 but there is still quite a bit of scatter and
no overall collapse for different pitching amplitudes is observed.

92

Figure 4.10: Circulation normalized following the scaling of Buchholz et al. (2011). Circulation
scaling using a β = 7 plotted against Stte

The measurement of rc is very sensitive to the domain of integration. In the current study, the
domain was fixed at 0.083c for comparison with the previous experiments of Bohl (2002), Bohl &
Koochesfahani (2009), and the computations of Hammer (2016) (which was discussed previously
in Section 4.1.3). The fixed domain size may be biasing the measurement for higher amplitude
cases as rc is shown to increase with increasing Îą0 (see Figure 4.11a). A less sensitive metric of the
core size is one based on a direct fit to a 2D Gaussian function of the vorticity field, designated at rÎł .
The two different core size metrics are plotted against Stte in Figure 4.12a-b for rc and rÎł . The core
size estimated by rÎł consistently computes a smaller core size compared to rc . While neither rc
p
nor rÎł collapse the data when plotted against Stte , when the core size is scaled by rÎł / h p /c, better
collapse of the data across Îą0 is observed when plotted against Stte for both the rigid and flexible
p
airfoils. Figure 4.12c shows the original rc (radius of gyration) with the 1/ h p /c scaling plotted
against Stte but it does not collapse the data well, again supporting the notion that the radius of
gyration with a fixed integration domain is a very sensitive parameter and fails to capture all of the
details of the vortex core size in experiments with different oscillation amplitudes. The influence

93

of flexibility is again captured by considering the actual trailing edge amplitude as an important
scaling parameter to to collapse the vortex core size across all pitching amplitudes (Îą0 ) considered.

(a)

(b)

Figure 4.11: Vortex core size plotted against k and St. The vortex core size, rc , plotted against k
(a), and St (b) for both rigid and flexible airfoils.

94

(a)

(b)

(c)

(d)

Figure 4.12: Core size of rigid and flexible airfoils plotted against Stte . The vortex core size
calculated by the radius of gyration, rc , (a,c), and calculated from a 2D
p Gaussian fit to vorticity
field, rÎł , (b,d), all plotted against Stte . The core size is scaled by 1/(c h p /c) in (c) and (d). The
legend shown in (d) is common to all plots (a-d).

95

4.3.4

Convective Speed and Wavelength

The convective speed, Uc , of the vortices is calculated from tracking the location of ω pk over
x/c = 0.5 ± 0.25 and finding a linear fit the peak’s location vs. phase. The wavelength is
computed from a = Uc / f . Estimating the wavelength in this way is warranted for two reasons: (1)
it provides a localized estimate of the wavelength, and (2) removes the ambiguity that arises when
wanting to extract the vortex properties close to the airfoil. If the vortex formation location (x f ) is
farther downstream than xm − a/2 where xm is the extraction location, it would be impossible to
measure a centered at xm . Estimating the wavelength by a = Uc / f circumvents this issue and has
been used previously in Monnier et al. (2015).
The convective speed is plotted against k, St, and Stte in Figure 4.13 and acceptable collapse in
the data is observed for Uc plotted against Stte , but the collapse is not as good as what was observed
for Γ0 and ω pk . Monnier et al. (2015) observed a similar amount of scatter in their data for
Uc /U∞ plotted against Stte , which can be seen in Figure 4.14. However there is a more significant
discrepancy in Uc in the current study compared to Monnier et al. (2015), with the current study
agreeing more closely with Bohl & Koochesfahani (2009), and Hammer (2016) which are also
plotted in Figure 4.14. The discrepancy between the studies may be due to Reynolds number and
airfoil shape differences. The airfoil used in Monnier et al. (2015) differs significantly from a the
NACA 0009 airfoil used in the current study; its geometry was similar to Heathcote & Gursul
(2007) with a rigid “head” and a flexible “tail”. The Reynolds number in the study of Monnier et al.
(2015) also was much lower than the current study, Rec =2.3×103 for compared to Rec =1.2×104
for the current study.
The wavelength normalized by the chord, a/c, is plotted against k, St, and Stte in Figure
4.15a-c. The wavelength shows better collapse with k than St and scatters even more with Stte
(since Stte ≥ St for all cases in the current study). However, referring to Figure 4.15b, a/c is seen
to group well with Îą0 with the flexible airfoil deviating systematically away from the rigid airfoil
as St increases. This motivates seeking a more appropriate scaling of a, and h p is certainly an
appropriate length scale. Figure 4.15d shows the a/h p collapses the data when plotted against Stte

96

(a)

(b)

(c)

Figure 4.13: Vortex convective speed of the rigid and flexible airfoils. The convective speed, Uc ,
normalized by the freestream, U∞ , and plotted against k (a), St (b), and Stte (c).

97

Figure 4.14: Vortex convective speed of the rigid and flexible airfoils compared to selected
literature. Vortex convective speed of the rigid and flexible airfoils compared to selected literature.
The convective speed normalized by the freestream, Uc /U∞ , and plotted against Stte . The data
from Bohl & Koochesfahani (2009), Monnier et al. (2015), and Hammer (2016) are included. The
Rec of Monnier et al. (2015) is 2.3 × 103 compared to nominally 1.2 × 104 for the rest of the studies.

for all Îą0 for both rigid and flexible airfoils. With very good collapse seen for a/h p , it is reasonable
to go back to convective speed and re-scale it. Scaling the convective speed as Uc /(U∞ Stte ) would
show considerably better collapse in data compared to what was observed previously in Figure
4.13 however, this scaling simply reduces Uc /U∞ Stte to a/h p since Uc = a f , resulting in no new
information.
This scaling of a/h p is tested with the results from Bohl (2002), Bohl & Koochesfahani (2009),
and Hammer (2016). Figure 4.16 shows a/h p plotted against Stte for these studies and the current
rigid and flexible airfoil measurements. The scaling of wavelength by h p collapses all the data
from all studies. Overall, this scaling is not a surprising outcome since b/h p will be shown in the
next section to quickly approach a constant value of 1 and, in their yet-to-be published follow-up
study, Bohl & Koochesfahani (2009) show b/a, the vortex array aspect ratio, tends to collapse to
a single value as well. However, this study demonstrates that yet again, the effect of chordwise
flexibility is captured as an amplitude effect that is quantified by measuring the actual trailing edge

98

(a)

(b)

(c)

(d)

Figure 4.15: Wavelength scaling. The wavelength is normalized by c and plotted against k (a), St
(b), and Stte (c). The wavelength normalized by h p and plotted against Stte (d).

99

Figure 4.16: Wavelength normalized by the trailing edge displacement.. Nondimensional wavelength, a/h p , plotted against Stte including the data from Bohl (2002), Bohl & Koochesfahani
(2009), and Hammer (2016).

displacement for scaling both Uc and a (since a = Uc f ).

4.3.5

Transverse Spacing

Two different scalings for the transverse spacing (b) of the vortices will be considered. The first,
b/h p is plotted against k, St, and Stte in Figure 4.17. Considering α0 = 2◦ , b/h p in Figure 4.17a is
shown to increase with increasing k up until b/h p ≈ 1, where it appears to asymptote to that value.
As mentioned previously, this effect was observed in the yet-to-be-published paper by Bohl and
Koochesfahani where they investigate the effect of amplitude using α0 = 2◦ and α0 = 4◦ data from
Bohl (2002). This asymptotic behavior is shown more clearly when the data is plotted against St
and Stte , as it is in Figure 4.17b-c. The higher amplitude cases fall onto one curve with the α0 = 2◦
data better for Stte as the abscissa compared to St. The limit of b/h p = 1 is reasonable from
the simplistic viewpoint of the vortex formation being connected with the extrema of transverse
position of the trailing edge. However, at low Stte , the vortex spacing is nearly −2h p , indicating
that at these cases, the transverse spacing is dependent on a more complicated mechanism related

100

to the when in the oscillation cycle the vortex actually begins to roll up.
In terms of flexibility, the flexible airfoil’s b/h p at k = 4, 5 and 6 is higher than its rigid
counterpart when oscillating with α0 = 2◦ . This is consistent with the flexible airfoil having a
lower thrust crossover than the rigid airfoil, which was shown in Chapter 3. When b/h p is plotted
against Stte this shift is no longer observed as both rigid and flexible airfoils appear to fall on the
same (steep) curve for Stte < 0.1 before reaching the asymptotic value.
(b)

(a)

(c)

Figure 4.17: Transverse spacing normalized by trailing edge displacement of the rigid and flexible
airfoils. Nondimensional transverse spacing, b/h p is plotted against k (a), St (b), and Stte (c).
Figure 4.18 shows b/h p plotted against Stte for both the rigid and flexible airfoils of the current
101

study, the experiments of Bohl (2002) and Bohl & Koochesfahani (2009), and the simulations
of Hammer (2017). The transverse spacing normalized by the trailing edge displacement, b/h p ,
is consistent across all of the studies. Considering all of the studies presented, b/h p appears to
approach 1 around Stte = 0.15, but as Stte increases beyond this, b/h p appears to gradually decrease
to b/h p = 0.9 as Stte → 0.4.

Figure 4.18: Transverse spacing normalized by trailing edge displacement compared with selected
literature. The transverse spacing normalized by trailing edge displacement b/h p plotted against
Stte of the rigid and flexible airfoils are plotted with the additional studies of Bohl (2002), Bohl &
Koochesfahani (2009), and Hammer (2016).

The second scaling that is considered for the transverse spacing is b/a, shown in Figure 4.19,
which is the vortex array aspect ratio. When b is normalized by a, the data still remain collapsed
across Îą0 and between rigid/flexible airfoils, but does not asymptote as quickly as b/h p . When
plotted against Stte , b/h p asymptotes near Stte = 0.13 while b/a asymptotes at a higher range of
0.2 < Stte < 0.25 (the scatter in the data at higher Stte makes selecting this value difficult). The two
are not necessarily related either as the first relates to the trailing edge displacement, and the latter
is the aspect ratio of the free vortex array. It is however, interesting that b/h p asymptotes close
to the thrust to drag crossover Stte while the aspect ratio of the vortex array, b/a doesn’t reach its
asymptotic value until a much higher Stte .
102

(a)

(b)

(c)

Figure 4.19: Transverse spacing normalized by wavelength for the rigid and flexible airfoils. The
transverse spacing normalized by the wavelength, b/a plotted against k (a), St (b), and Stte (c).

103

Figure 4.20: Transverse spacing normalized by wavelength compared with selected literature. The
vortex array aspect ratio, b/a, is plotted against Stte with the additional studies of Bohl (2002),
Bohl & Koochesfahani (2009), Hammer (2016), Hammer (2017).

The asymptotic behavior of the transverse spacing normalized by the wavelength was shown
previously by Bohl & Koochesfahani (2009) to be b/a → 0.2 for α0 = 2◦ and Hammer (2016)
also observed the same asymptotic value of 0.2 for an extended range of α0 = 2◦ , 4◦ , and 6◦ . The
current study also observes this behavior, but over the extended Stte range the asymptotic value
of 0.2 is less clear; this can be seen by considering Figure 4.19c for Stte > 0.2, which shows the
asymptotic value is greater than 0.2. Results from both of these previous studies and Hammer
(2017) are plotted together in Figure 4.20 in addition to the current experiments. Considering all of
these studies, the asymptotic value for b/a does appear to be greater than 0.2 when considering the
larger range of Stte of the current experiments and the simulations of Hammer (2016) than those
of Bohl & Koochesfahani (2009).

104

4.4

Summary

In this chapter, the baseline (rigid airfoil) cases of the current study were validated against
previous experiments and computations, showing excellent overall agreement both qualitatively, in
terms the vorticity fields, and quantitatively, in terms of the extracted vortex parameters of cases
where a single pair of vortices were shed per oscillation. The effect of chordwise flexibility was
shown to be captured as an amplitude effect measured by the actual trailing edge displacement
(h p ). As a length scale, h p was shown to be useful for normalizing b, a, and rÎł , and in the
Strouhal number based on the actual trailing edge displacement, Stte . When Stte is used as the
abscissa, nearly all of the individual vortex properties and the vortex array properties collapse for
all oscillation amplitudes for both rigid and flexible airfoils. Table 4.1 summaries all of the vortex
configuration parameters and their dimensionless forms found to collapse the data for both rigid
and flexible airfoils across all pitching motions considered in this study.
Base Parameter

ordinate

abscissa

Γ0
ω pk
rÎł
b
b
a

Γ0 /(cU∞ )
ω pk c/U
p ∞
(rÎł /c)/ h p /c
b/h p
b/a
a/h p

Stte
kStte
Stte
Stte
Stte
Stte

Table 4.1: Vortex parameter nondimensionalizations and scaling variable. A summary of vortex
parameters and the corresponding parameters used to scale and collapse the data.

105

CHAPTER 5
CONCLUSIONS

In this work, the influence of chordwise flexibility on the thrust and vortical wake structure produced
by a harmonically pitching airfoil at a chord Reynolds number of Rec = 1.24 × 104 was investigated
experimentally. While this effect has received renewed interested recently, this study takes a
different approach from what is commonly seen in the literature (flat plates, membranes, flaps)
and studies a single geometric shape, a standard NACA 0009 airfoil, that is either rigid or flexible.
Molecular tagging velocimetry was used for planar two-component measurements of the velocity
field at the mid-span in the wake of the oscillating airfoil. A rigid and a chordwise flexible airfoil
were oscillated at multiple amplitudes and frequencies. The range of oscillation frequencies and
amplitudes considered in this study were such that the trailing edge displacement increased with
increasing k, and Îą0 was kept low enough such that leading edge vortices were not formed. The wake
vortex properties were quantified in terms of their strength (peak vorticity, ω pk , and circulation,
Γ0 ), core size (radius of gyration, rc , and 2d-gaussian core size, rγ ), and spacing (wavelength a,
and transverse spacing b) for oscillation frequencies and amplitudes combinations that resulted in
a von KĂĄrmĂĄm or reverse von KĂĄrmĂĄn vortex pattern.
The rigid airfoil results were compared against numerical simulations of a NACA 0009 and
both experiments and simulations of a NACA 0012 airfoil to establish a baseline of comparison
to previous studies. Both the individual vortex properties and mean velocity profiles were compared. While there are subtle differences between the thinner NACA 0009 airfoil compared to the
NACA 0012 used more extensively in the literature, the differences were not significant compared
to the effect of chordwise flexibility, the primary effect under consideration in this study.
Chordwise flexibility was shown to increase thrust production at a given oscillation frequency
and amplitude under the conditions considered in this study. The greater body of literature shows that
chordwise flexibility can increase or decrease the thrust production, depending on the oscillation
frequency and the flexibility of the airfoil. It was found that when plotted against the Strouhal
106

number based on the actual trailing edge displacement, Stte , that thrust collapses equally well for
both rigid and flexible airfoils for the range of parameters studied.
Little attention in the literature is given to the wake vortex properties of chordwise flexible
airfoils. Unique to this study, the properties of the vortices and their arrangement was studied
in detail for a chordwise flexible airfoil at a chord Reynolds number of Rec = 1.25 × 104 . The
scalings and non-dimensionalizations that collapsed the rigid airfoil data for different oscillation
amplitudes were shown, when plotted against Stte , to not only collapse the flexible airfoil, but to
collapse the data to the rigid airfoil data as well. The measured trailing edge displacement (h p )
was found to be the proper length scale for non-dimensionalizing, except in the case of circulation
(Γ0 ). When trailing edge displacement was used in scaling the transverse spacing b, wavelength
a, and core size (rÎł ) the data collapse when plotted against Stte for both rigid and flexible airfoils
for all oscillation amplitudes studied (α0 = 2◦ , 4◦ , and 6◦ ). The peak vorticity was also found
to collapse for both rigid and flexible airfoils for all Îą0 considered when plotted against kStte , a
non-dimensional term proportional to the maximum trailing edge acceleration. Overall, chordwise
flexibility was shown to act primarily through alteration of the trailing edge amplitude for the range
of parameters investigated through.

107

APPENDICES

108

APPENDIX A
TWO CAMERA CALIBRATION

The imaging system used for the current work utilizes two 1392×1024 pixel PCO.Pixelfly cameras
which are mounted underneath the water tunnel as shown in Figure A.1. The camera assembly is
connected to the streamwise and vertical traverses such that it also moves with the entire optical
assembly. It is on a separate cross-stream traverse since it must move at a different rate than the
optical system due to the refraction of the UV beam passing through the quartz windows at an
angle. A single 105 mm f/2.0 Nikon Nikkor Lens, a shortened F-mount to C-mount adaptor, and
a standard 1” C-mount cube beam splitter from Edmund Optics completes the setup. The compact
two-camera system relies on a calibration routine based on the Affine transform to compensate for
large (greater than 10 pixel) sensor misalignments, which is described in detail in below.
An affine transformation is defined as any transformation that preserves collinearity and ratios
of distances; in general an affine transformation is a composition of rotations, translations, dilations,
and shears (Shapiro & Stockman (2001)). The use of affine transformations is very common in
image processing and is generalized by the Equation A.1 where x 0 is the transformed coordinates
and x is the original coordinates.
x 0 = T(x) = Ax + B

(A.1)

Interpolation is typically employed after computing the transformed coordinates to create output a
new image with the original spacing. The most common interpolation methods are either nearest
neighbor, bi-linear, and bi-cubic.
A single affine transformation needs to be determined to remap the image taken from one
camera into the pixel coordinate system of the second camera. The transform is found by the
linear least squares method using a large number of known control points that are imaged from
both cameras simultaneously (as seen in Figure A.2), and determining the displacement of each
control point (grid intersection) using our 2c-MTV correlation program. The displacement field,
shown in Figure A.3, indicates a slight translation, but is dominated by angular rotation between the
109

Figure A.1: MTV Imaging System. Photograph of dual camera MTV imaging system. Two
PCO.Pixelfly cameras (blue) are flush mounted to a backing plate and a cube beam splitter (highlighted in red). A single lens is mounted to a shortened c-mount to Nikon F-mount adapter
(highlighted green). The trigger boxes can be seen in the upper right of the image. The cameras
are located underneath the water tunnel looking up towards the test section. The gray pipe in the
bottom of the image is the return pipe of the water tunnel.

two sensors. The traditional multi-camera approach would have used mechanically corrected these
errors using translation and rotational stages such that there was less than 1 px misalignment for
all control points. However, in this application, the computed displacements of the control points
are used to determine the optimal affine transformation using the built-in functions in MATLAB’s
image processing toolbox.
The application of the affine transform to overcome the misalignment of the two camera’s
sensors is assessed by re-computing the displacement field shown in Figure A.3 with the original
image of Camera #1 (Figure A.2 left) and the transformed image of Camera #2 (not shown). If
the true transformation between the imaging planes were Affine (they aren’t due to higher order

110

Figure A.2: Targets of each camera. (left) Camera #1 target image (right) Camera #2 target image.
The red line is straight relative to the paper; the relative shift may be seen by comparing the distance
of an intersection to the red line on both images.

distortions) and there was no error in the displacement calculations of the control points (there
isn’t due to imaging and processing noise), a perfect transformation should achieve a mean of
zero displacement across the entire image. A mean zero displacement to the 1st decimal place is
achieved as shown in Figure A.4 with an overall RMS (of 256 images and 357 intersections) of 0.04
and 0.01 pixels in the horizontal and vertical directions respectively. The displacement field still
has a distinct pattern, indicating there may be room for improvement in either the affine transform
parameters or by using a higher order transform to drive the error further to zero. The perspective
transform was also tested but did not perform significantly better than the simpler affine transform.
However, this residual displacement field is now a known bias error and can be subtracted from
other results.
In order to demonstrated the feasibility to use this technique to correct for the sensor misalignment a known target is displacement and imaged by both cameras simultaneously. The true
displacement of the target can be computed from correlating the target at location 1 to location 2
independently from each camera. When this is performed, the displacement of the target is: Âąx =
0.75 pixels and ∆y = 15.23 pixels. Now, the affine transformation computed previously is used to
transform one set of images and the target image is correlated at location 1 from one camera to the

111

Mean ∆x:
-1.70 px
Std ∆x: +2.565 px
Mean ∆y: +1.036 px
Std ∆y: +3.501 px

Figure A.3: Uncorrected displacement field between the two cameras. Initial displacement field.
The vector field shows the displacement from Camera #1 (c1) at Location 2 (L2) to Camera #2 (c2)
also at Location 2 (L2). The displacement field is computed by correlating the two images together
that were looking at the same target as shown in Figure A.2. Reference vector shown in upper right
of plot.
transformed images of camera 2 at location 2. The result yields an error in the mean displacement
field of 0.07% in x (0.01 pixels), and 0.4% in y (0.61 pixels), with RMS values of 0.047 pixels in
both directions. Now, the residual displacement field due to the not-perfect image transformation,
as shown in Figure A.4, is subtracted. The results are more than satisfactory: -0.060% in x (<0.01
pixels), and -0.024% in y (<0.04 pixels), with 0.019 pixel and 0.031 pixel RMS in spatial variation,
respectively.
These results demonstrate that the sensor misalignment in the current setup can be completely
corrected for by applying the affine transform determined from a calibration target to one camera’s
images prior to correlation and subtracting the residual misalignment to the displacement fields
with a subpixel accuracy of 0.01 Âą 0.02 pixels RMS. The results can be improved with a higher
spatial density calibration target such that residual displacement field is more finely known; which
is done in all further experiments where the number of control points is increased by over a factor
of 2 of the those computed here.

112

Mean ∆x:
Std ∆x:
Mean ∆y:
Std ∆y:

-0.013 px
+0.043 px
+0.005 px
+0.014 px

Figure A.4: Corrected displacement field using affine transform. Corrected displacement field
using affine transform. A reference vector, 0.1 px in length, is displayed in the upper right of the
plot. The vector field shows the displacement from Camera #1 at Location 2 to the transformed
image of Camera #2 also at Location 2. This displacement field is the residual error left after
applying the affine transformation. Reference vector shown in upper right of plot.

113

APPENDIX B
UNCERTAINTY ESTIMATES

B.1

Vorticity

The uncertainty analysis outlined here for the phase averaged spanwise vorticity extends what
was presented in Bohl & Koochesfahani (2009) based on the fourth order accurate finite difference
scheme discussed in Cohn & Koochesfahani (2000). The current analysis includes the spatial
dependence of the uncertainty of the velocities. The spanwise vorticity, ωz , is calculated using a
fourth order accurate central finite difference scheme, given by:

 
∂v
∂u
(ωz )i, j =
−
∂ x i, j
∂ y i, j


where


and


(B.1)


−vi+2, j + 8vi+1, j − 8vi−1, j + vi−2, j
∂v
,
=
∂ x i, j
12h

(B.2)


−ui, j+2 + 8ui, j+1 − 8ui, j−1 + ui, j−2
∂u
.
=
∂ y i, j
12h

(B.3)

The data spacing, h, is a known constant; therefore, the vorticity is a function of four transverse
velocities and four streamwise velocities; each with uncertainty v and u , respectively. The uncertainty of each of the transverse and streamwise velocity components are a function of both the
subpixel accuracy of the velocimetry and the variation in the phase averaging process. Standard uncertainty analysis allows the determination of the uncertainty in spanwise vorticity, ωz , according
to
ωz


2

2
∂ωz
=
 χn
∂ χn


(B.4)

where χn and  χn represent the nth independent variable and its uncertainty, respectively. Carrying
out the uncertainty analysis using the vorticity definition and the finite difference formulation of

114

velocity derivatives given above leads to the final result
2 
2 
2 
2


2
−1
2
−2
1
ωz i, j =
v
v
v
v
+
+
+
12h i+2, j
3h i+1, j
3h i−1, j
12h i−2, j
2 
2 
2 
2

−1
2
−2
1
+
u
u
u
u
+
+
+
12h i, j+2
3h i, j+1
3h i, j−1
12h i, j−2

B.2

(B.5)

Circulation

The circulation of a vortex, Γ0 , is calculated by numerical integration (for a single sign of
vorticity at a time), given by:
Γ0 =

ÕÕ
i

ωi, j dA,

(B.6)

j

where the differential area, dA, is a known constant, and ωi, j is summed only for the region of
interest, typically a radial distance of R from the center of the vortex. The uncertainty in circulation
is given by:
( Γ0

)2

=

ÕÕ 
i



ωi, j dA ωi, j

 2

,

(B.7)

j

where the uncertainty of each measurement of ωi, j is given by Equation B.5. The quadrature
summation is derived from standard uncertainty analysis assuming each measurement of ωi, j is
independent.

B.3

Radius of Gyration

The vortex core sized is calculated by the radius of gyration, rc , given by:
Í Í 2
i j r i, j ωi, j
,
rc = Í Í
i j ωi, j dA

(B.8)

where r is the distance from the center of the vortex and the integration domain is a finite size of
radius R. The uncertainty rc is derived from standard uncertainty analysis:

2
∂rc
2
(rc ) =
 χn ,
∂ χn

(B.9)

where χn and  χn are the nth independent variable and its corresponding uncertainty. In the case
of the radius of gyration the independent variables are ωi, j located within the integration radius R
115

surrounding the center of the vortex. The differentiation of Equation B.8 with respect to each ωi, j
is:

2
∂rc
dA © Γ0ri, j − 1 ª
= √ ­
ÂŽ,
∂ωi, j 2 rc
Γ02
ÂŤ
ÂŹ
therefore, the uncertainty in radius of gyration is:

Õ Õ ©© dA
(rc )2 =
­­ √
2 rc
i
j ÂŤÂŤ

B.4

2

(B.10)

2

© Γ0ri, j − 1 ªª
ÂŞ
­
®® ωi, j ® .
2
Γ0
ÂŤ
ÂŹ
ÂŹÂŹ

(B.11)

Control volume method for estimating the streamwise force

The control volume method for computing the streamwise force is described by Bohl & Koochesfahani (2009), whose improvements in using the fluctuating velocity profiles demonstrated a marked
improvement compared to the traditional analysis of only using the mean velocity profile. In the
context of the current measurements, the integral form of the equation is approximated using the
numerical quadrature method at a specific downstream location of i = I and over the transverse
region of j | y = ÂąH denoted by j = [0, M]. This leads to the expression:
!
!
(
M
hui I, j
2 Õ hui I, j hui I, j
Ct =
−1 +ξ
−1
c
U∞
U∞
U∞
j=1
!)
!2
!2
hui0 I, j
hvi0 I, j
Uo2
1
+
−
+
∆y,
1−
2
U∞
U∞
2
U∞

(B.12)

where


1
U0
Ξ=
1−
.
2
U∞
Following standard uncertainty analysis the uncertainty of Ct is expressed as

2
∂Ct
2
(Ct ) =
 χn ,
∂ χn

(B.13)

(B.14)

where the partial derivatives of Ct with respect to each of the independent variables ( χn ) needs to
be computed. The independent variables are: hui I, j , hui0I, j , hvi0I, j , and U0 . The partial derivatives
of Ct with respect each variable are:
∂Ct
∂huii, j

=

4huii, j − U0
2
U∞

116

!
1
−
∆y
2U∞

(B.15)

!
2 hui0i, j
∂Ct
=
∆y
2
∂hui0i, j
U∞
!
2 hvi0i, j
∂Ct
= −
∆y
2
∂hvi0i, j
U∞
!
hui
∂Ct
1
= −
∆y
+
2
∂U0
2U∞
2U∞

(B.16)

(B.17)

(B.18)

The uncertainty of the mean streamwise velocity, huii, j , can be estimated by the mean standard
error,

1.96 hui0i, j
,
 huii, j =
√
nφ

(B.19)

where nφ is the number of phase bins used in creating the phase averaged time series (64 in the
current study). The factor of 1.96 is used for the error to represent the 95% confidence interval,
which assumes the variable has a normal distribution. Similarly, the error in U0 is estimated by the
mean standard error since it it is estimated from averaging hui over p samples at the extent of the
mean velocity profile as given by
p

M

Õ
1 ©Õ
ÂŞ
hui I , j +
hui I , j ÂŽ .
U0 =
(B.20)
­
2p
j=0
j=M−p+1
ÂŤ
ÂŹ
In the current study, p = 5, and the standard deviation of U0 is the standard deviation of hui over
the region described above, denoted as U0,rms . This leads to the uncertainty in U0 being
U0 =

(1.96)U0,rms
p
2p

(B.21)

for a 95% confidence interval.
The uncertainties in hui0I, j and hvi0I, j required a more detailed examination of the underlying
flow and its measurement. The measured velocity is
u(x, y, t) = u(x, y) + uφ (x, y, φ) + u01 (x, y, t) + u02 (x, y, t) + u03 (x, y, t)

(B.22)

where u is the time average, uφ is the oscillatory velocity fluctuation, u01 is the cycle to cycle variation,
u02 is the random fluctuations from turbulence, and u03 is the fluctuations due to measurement error.
Phase averaging this by
hui (x, y, φ) =

∫ φ
2
φ1

117

(u(x, y)) dφ

(B.23)

yields, at each phase bin φ,
uφ (x, y, φ) = hui (x, y, φ) + hui0φ ,

(B.24)

where hui0φ is the phase-bin-rms. The uncertainty in hui0i, j should be computed using the average
variance from all of the individual phase-bin-rms results,
r
1 Õ  0 2
hui φ,i
Shui0 =
n
i, j

(B.25)

where the actual confidence interval is based on the precision interval of the sample variance. The
following expression is given in Figliola & Beasley (2006):
νSx2
2
χ0.025

≤ σ2 ≤

νSx2

(B.26)

2
χ0.975

where ν is the degrees of freedom, Sx2 is the measured sample variance of measured variable x
(in this case either hui0I, j or hvi0I, j ), and σ 2 is the true variance of the population. The average
difference between the upper and lower bounds and the measured Sx2 is used such that it can be cast
into the the uncertainty analysis being carried out, which utilizes a single (Âą) uncertainty value.
The final expressions for uncertainty in hui0I, j and hvi0I, j are therefore:
ν ∗ Shui0
1Š
I, j
ÂŞ
− hui0I, j ® +

= ­
2
χ2
ÂŤ 0.025
ÂŹ
ν ∗ Shvi0
1Š
I, j
ÂŞ
 hvi0I, j = ­
− hvi0I, j ® +
2
2
χ
ÂŤ 0.025
ÂŹ
where ν = (nφ − 1), which is 63 in the current study.
hui0I, j

1Š 0
­ hui I, j −
2
ÂŤ
1Š 0
­ hvi I, j −
2
ÂŤ

ν ∗ Shui0

I, j ÂŞ

2
χ0.975

ν ∗ Shvi0

ÂŽ
ÂŹ

I, j ÂŞ

2
χ0.975

ÂŽ
ÂŹ

(B.27)

(B.28)

The final uncertainty (95% confidence interval) of Ct is therefore the summation of each of the
partial derivatives multiplied by each of their corresponding uncertainties, as previously indicated
by Equation B.4. It should be noted that this indicates that each term is summed over the extent
of the profile (i.e. summation over j for a single value of i = I), which in the current study is 121
elements, for a total of 484 total terms for the four flow variables used in computing Ct .
Overall, the uncertainty estimates computed here are compared to two sets of mean thrust
coefficient measurements which were acquired from independent data sets. Figure B.1 shows the
118

two sets of measurements plotted against reduced frequency. The two sets of measurements, agree
with each other within the computed 95% confidence intervals.

Figure B.1: Independent measurements of Ct plotted against reduced frequency (k). Comparison
of independent measurements of Ct plotted against reduced frequency.

119

APPENDIX C
EFFECT OF VORTICITY CUTOFF LEVEL

A vorticity cutoff level (ωcut ) is used throughout this study to minimize the background noise in
the vorticity from impacting parameters calculated from the vorticity fields such as circulation and
the vortex core size. The cutoff functions as:
ωz (x, y) = 0, where: |ωz (x, y)| ≤ ωcut .

(C.1)

A convenient starting point for investigating the sensitivity of the vorticity cutoff is considering an
isolated Gaussian vortex whose radial (r) vorticity distribution is given by:
2

ω(r) = ω pk e−(r/rc ) , r ≥ 0,

(C.2)

where ω pk , is the peak vorticity and rc , is the core radius. The core radius in this context is defined
as the radial distance at which the vorticity has reduced to e−1 of its initial value (nominally 37%
of ω pk ). In general, circulation is defined as the line integral around a closed curve of the velocity
and application of Stokes’ Theorem yields the relationship between circulation and vorticity as the
area integral of the vorticity.
Γ=

∎
∂S

V ¡ dl =

∏

ω · dS

(C.3)

S

The circulation of the isolated Gaussian vortex is computed from integrating Equation C.2 as
r → ∞ to yield
Γ0 = πω pk rc2,

(C.4)

and the radial distribution of circulation is
Γ(r) = Γ0



2
1 − e−(r/rc )



,

(C.5)

which may also be expressed as
Γ(r) =

πω pk rc2



2
1 − e−(r/rc )

120



.

(C.6)

These relationships show the interconnectedness of the peak vorticity, core radius, and circulation
of a Gaussian vortex, and are useful in relating more complicated vortical flow fields to idealized
vortices. Currently, these relationships are exploited to examine the sensitivity to a cutoff level of
vorticity and its impact on the calculation of circulation and core size. Any vorticity cutoff is shown
(Figure C.1) to reduce the calculated values of Γ and rc , with a greater impact on rc . The reduction
in circulation is linear with a slope of 1 over the extent examined such that ωcut = 0.04ω pk yields
a 4% error in Γ0 and a 6% error in rc .

Figure C.1: Effect of vorticity cutoff on Γ and rc . Effect of vorticity cutoff on calculated circulation
and radius of gyration of a model Gaussian vortex. The cuttoff vorticity level always decreases both
circulation and radius of gyration. The circulation will be reduced by 1:1 and radius of gyration
1:1.4 based on the vorticity cutoff’s percentage of the peak vorticity over for up to 10% of the peak
vorticity.

The effect of the vorticity cutoff is considered for three different vorticity cutoff levels on the
rigid airfoil data and summarized in Table C.1. The change in Γ0 for three different levels (0, 3, and
5 non-dimensional vorticity) is consistent with the relationship established above for an isolated
Gaussian vortex based on the cutoff level and the peak vorticity. The change in rc is more sensitive
121

than predicted for the highest cutoff level (5), but closer for the lower level of (3). The increased
sensitivity to rc is contributed to the spatial resolution of the measurements (1%) compared to rc
varying between 3% and 4% of the chord. Overall for the cutoff level of 3 used in this study rc
varies less than 12% for α0 = 2◦ and less than 5% for both α0 = 4◦ and α0 = 6◦ compared to using
either no cutoff or a level of 5. Similarly, Γ varies less than 7% for α0 = 2◦ and less than 4% for
both α0 = 4◦ and α0 = 6◦ .

α0 = 2◦

α0 = 4◦

α0 = 6◦

k
5
6
7
8
9
10
12
15
5
6
7
9
10
4
5
6
8

St
0.08
0.10
0.12
0.13
0.15
0.17
0.20
0.25
0.17
0.20
0.25
0.30
0.33
0.20
0.25
0.30
0.40

ω pk
43
68
80
111
122
143
192
253
90
134
167
223
250
96
145
166
243

ωcut = 0
0.037
0.036
0.034
0.033
0.034
0.034
0.033
0.033
0.041
0.039
0.040
0.040
0.040
0.044
0.044
0.044
0.044

rc
ωcut = 3
0.035
0.035
0.031
0.031
0.031
0.032
0.032
0.032
0.040
0.039
0.040
0.040
0.040
0.043
0.044
0.044
0.044

Γ
ωcut = 5
0.033
0.033
0.028
0.028
0.029
0.029
0.030
0.031
0.038
0.037
0.039
0.039
0.040
0.042
0.043
0.044
0.044

ωcut = 0
0.15
0.18
0.21
0.25
0.29
0.33
0.42
0.56
0.34
0.42
0.57
0.73
0.84
0.40
0.55
0.71
1.03

ωcut = 3
0.14
0.18
0.20
0.24
0.28
0.32
0.41
0.56
0.33
0.41
0.57
0.73
0.84
0.40
0.55
0.71
1.02

ωcut = 5
0.13
0.17
0.18
0.23
0.27
0.31
0.40
0.55
0.32
0.40
0.57
0.73
0.84
0.39
0.55
0.71
1.02

Table C.1: Effect of vorticity cutoff level on rc and Γ. Circulation and radius of gyration for
different vorticity cutoff levels. The peak vorticity for each case is also shown for comparison to
the cutoff level. The average of the positive and negative vortex is displayed and all parameters
have been non-dimensionalized (i.e. ω pk c/U∞ , Γ/(cU∞ ), rc /c).

122

APPENDIX D
STREAMWISE FORCES

Included in this Appendix are several figures that show the mean streamwise thrust coefficient
plotted for k, St, and Stte for either the full range of the data or for a selected range of the data.
These plots supplement the range of the parameters that were included in Chapter 3.

Figure D.1: Mean thrust coefficient plotted against k, full range. Mean thrust coefficient plotted
against k for the full range of the data.

123

Figure D.2: Mean thrust coefficient plotted against k, selected range. Mean thrust coefficient
plotted against k for a selected range of data.

124

Figure D.3: Mean thrust plotted against St, full range. Mean thrust coefficient plotted against St
for the full range of the data.

125

Figure D.4: Mean thrust coefficient plotted against St, selected range. Mean thrust coefficient
plotted against St for a selected range of data.

126

Figure D.5: Mean thrust coefficient plotted against Stte , full range. Mean thrust coefficient plotted
against Stte for the full range of the data.

127

Figure D.6: Mean thrust coefficient plotted against Stte , selected range. Mean thrust coefficient
plotted against Stte for a selected range of data.

128

APPENDIX E
VORTICITY FIELDS FOR RIGID AND FLEXIBLE AIRFOILS

The vorticity fields for the rigid and flexible airfoils are compared at selected reduced frequencies
for each of the different oscillation amplitudes considered. Figure E.1 compares the phase averaged
vorticity field of the rigid and flexible airfoils for four different reduced frequencies at α0 = 2◦ .
Figure E.2 compares the phase averaged vorticity field of the rigid and flexible airfoils for four
different reduced frequencies at α0 = 4◦ . Figure E.3 compares the phase averaged vorticity field
of the rigid and flexible airfoils for four different reduced frequencies at α0 = 6◦ .

129

α0 = 2◦
Rigid

Flexible

Figure E.1: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 2◦ . The phase
averaged vorticity at the same phase of the rigid (left) and flexible (right) airfoils. The reduced
frequency increases from top to bottom, with the reduced frequency listed as the title of each plot in
addition to the Strouhal number (St) and the Strouhal number based on the measured trailing edge
displacement (Stte ). The rigid and flexible plots have identical colormaps for the same reduced
frequency, but the colormap limits increase as k increases.

130

α0 = 4◦
Rigid

Flexible

Figure E.2: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 4◦ . The phase
averaged vorticity at the same phase of the rigid (left) and flexible (right) airfoils. The reduced
frequency increases from top to bottom, with the reduced frequency listed as the title of each plot in
addition to the Strouhal number (St) and the Strouhal number based on the measured trailing edge
displacement (Stte ). The rigid and flexible plots have identical colormaps for the same reduced
frequency, but the colormap limits increase as k increases.

131

α0 = 6◦
Rigid

Flexible

Figure E.3: Phase averaged vorticity field of rigid and flexible airfoil for α0 = 6◦ . The phase
averaged vorticity at the same phase of the rigid (left) and flexible (right) airfoils. The reduced
frequency increases from top to bottom, with the reduced frequency listed as the title of each plot in
addition to the Strouhal number (St) and the Strouhal number based on the measured trailing edge
displacement (Stte ). The rigid and flexible plots have identical colormaps for the same reduced
frequency, but the colormap limits increase as k increases.

132

APPENDIX F
POSITIVE AND NEGATIVE VORTEX PROPERTIES

The average value of the vortex properties of Γ0 , ω pk , a, Uc , and rc were reported in this study
rather than the values for each individual positive and negative vortex. For most vortex properties
the differences between the positive and negative vortex are less than the experimental uncertainty.
When the differences do exceed the experimental uncertainty, the differences are small relative to
the change of the variable seen in the parameter space. Overall, the use of the average compared
to the individual positive and negative vortices do not impact the interpretation of the data of this
study.
The peak vorticity extracted from the both the positive and negative vortices are plotted against
k in Figure F.1a for the rigid airfoil and Figure F.1b for the flexible airfoil. Overall, 10 out of
16 cases for the rigid airfoil and 11 out of 15 cases, 67% of all cases studied had no significant
difference between the positive and negative vortices. The peak vorticity did not vary by more than
10% for 87% of all of the cases.
Similarly, Figure F.2 shows the circulation extracted for both positive and negative vortices
plotted against k in for both the rigid and flexible airfoils. There is less difference between the pair
of vortices in terms of their circulation than what is observed in was observed for the peak vorticity.
Over 80% of all cases have less than a 4% difference in circulation between the pair of opposite sign
vortices. The maximum deviation was 18%, which appears relatively high due to how circulation
is scaled.

133

(a)

(b)

Figure F.1:
 
 Peak
  vorticity of both positive and negative sign vortices. The magnitude of the peak
vorticity,  ω pk , for both the positive and negative vortices are plotted against k for the rigid (a)
and flexible (b) airfoils.

134

(a)

(b)

Figure F.2: Circulation of both positive and negative sign vortices. The magnitude Circulation,
|Γ0 |, for both the positive and negative vortices are plotted against k for the rigid (a) and flexible
(b) airfoils.

135

APPENDIX G
SPATIAL RESOLUTION AFFECT ON VORTEX PROPERTIES

The numerical simulations of Hammer (2017) on the NACA 0009 airfoil are used to examine the
impact of spatial resolution on the extracted vortex parameters of peak vorticity, ω pk , circulation,
Γ0 , core size rc , transverse spacing b, wavelength, a, and the vortex array aspect ratio, b/a. The
resolution of the numerical simulations is 0.0025c in both the x and y directions in the wake of the
airfoil. The spatial resolution of the experiment is 0.01c in both x and y directions. The vortex
properties are computed, using the same methods outlined in Section 2.5.5 for the full resolution
numerical simulation results, and also for a reduced resolution of the numerical simulations. The
reduced resolution of the numerical simulations was achieved by keeping only every 4th data grid
point, and computing the vorticity on this new grid, the same data spacing as the experiments.
The percent change in the vortex properties between the full resolution and lower resolution
results are summarized in Table G.1. Peak vorticity will be attenuated due to the lower spatial
resolution, as well as the circulation. The change in b, a, and rc is not deterministic due to the
lower spatial resolution. Overall, ω pk and b are affected the most by spatial resolution, leading to
errors of 1%-8% for the reduced frequencies considered. The error in circulation is less than 3%
for k=4 and less than 1% for the rest of the cases. The error in wavelength, a, is less than 1% for all
cases. The error in b/a is driven by the error in b, and is between 1% and 8%. The effect on both
the (+) and (-) sign vortices is included in the table but don’t significantly differ from one another.

136

k

∆b/c
(%)

∆a
(%)

∆b/a
(%)

(+) ∆ω pk
(%)

(-) ∆ω
(%)

(+) ∆Γ0
(%)

(-) ∆Γ
(%)

(+) ∆rc
(%)

(-) ∆rc
(%)

4
6
8
10

4.72
7.70
3.38
-1.08

0.86
0.08
0.35
0.08

3.83
7.61
3.01
-1.16

-1.53
-4.36
-7.77
-3.85

-1.65
-4.47
-7.88
-3.51

-2.78
-0.39
-0.74
-0.27

-2.93
-0.34
-0.41
-0.38

-4.06
0.09
-1.80
-0.51

-3.93
0.16
-0.96
-0.89

Table G.1: Change in the vortex properties due to reduced spatial resolution. Change in the
vortex properties due to reduced spatial resolution. Using the numerical simulations of Hammer
(2017), the % change in the vortex parameters between the full resolution (0.00025x/c) and reduced
resolution equal to the experiments (0.001x/c).

137

BIBLIOGRAPHY

138

BIBLIOGRAPHY

Abbot, Ira H. & Albert E. von Doenhoff. 1959. Theory of Wing Sections: Including a Summary of
Airfoil Data. Mineola, NY, USA: Dover Publications Inc. 2nd edn.
Akkala, James M., Azar Eslam Panah & James H.J. Buchholz. 2015. Vortex dynamics and
performance of flexible and rigid plunging airfoils. Journal of Fluids and Structures 54. 103–
121.
Alben, Silas. 2008. Optimal flexibility of a flapping appendage in an inviscid fluid. Journal of
Fluid Mechanics 614. 355.
Betz, A. 1912. Ein beitrag zur erklaerung des segelfluges. Z. Flugtech. Motorluftsch 3. 269–270.
Bohl, Douglas G. 2002. Experimental Study of the 2-D and 3-D Structure of a Concentrated Line
Vortex Array: Michigan State University dissertation.
Bohl, Douglas G. & Manoochehr M. Koochesfahani. 2009. MTV measurements of the vortical
field in the wake of an airfoil oscillating at high reduced frequency. Journal of Fluid Mechanics
620. 63–88.
Buchholz, James H. J., Melissa A. Green & Alexander J. Smits. 2011. Scaling the circulation shed
by a pitching panel. J. Fluid Mech 688(2011). 591–601.
Calderon, D. E., D. J. Cleaver, I. Gursul & Z. Wang. 2014. On the absence of asymmetric wakes for
periodically plunging finite wings On the absence of asymmetric wakes for periodically. Physics
of Fluids 26(071907).
Chopra, M. G. 1976. Large amplitude lunate-tail theory of fish locomotion. Journal of Fluid
Mechanics 74(01). 161.
Cleaver, D. J., D. E. Calderon, Z. Wang & I. Gursul. 2016. Lift enhancement through flexibility of
plunging wings at low Reynolds numbers. Journal of Fluids and Structures 64. 27–45.
Cleaver, D. J., I. Gursul, D. E. Calderon & Z. Wang. 2014. Thrust enhancement due to flexible
trailing-edge of plunging foils. Journal of Fluids and Structures 51.
Cohn, R. K. & M. M. Koochesfahani. 2000. The accuracy of remapping irregularly spaced velocity
data onto a regular grid and the computation of vorticity. Experiments in Fluids 29(7). S061–
S069.
Combes, S. A. 2003. Flexural stiffness in insect wings I. Scaling and the influence of wing venation.
Journal of Experimental Biology 206(17). 2979–2987.
Dewey, P. A., B. M. Boschitsch, K. W. Moored, H. a. Stone & a. J. Smits. 2013. Scaling laws for
the thrust production of flexible pitching panels. Journal of Fluid Mechanics 732. 29–46.

139

Dong, H., R. Mittal & F. M. Najjar. 2006. Wake topology and hydrodynamic performance of
low-aspect-ratio flapping foils. Journal of Fluid Mechanics 566. 309.
Egan, Brendan C., Cody J. Brownell & Mark M. Murray. 2016. Experimental assessment of
performance characteristics for pitching flexible propulsors. Journal of Fluids and Structures
67. 22–33.
Figliola, Richard S. & Donald E. Beasley. 2006. Theory and Design for Mechanical Measurements.
Hoboken, New Jersey: John Wiley & Sons Inc. 4th edn.
Freymuth, Peter. 1988. Propulsive vortical signature of plunging and pitching airfoils. AIAA
Journal 26(7). 881–883.
Garrick, I. E. 1937. Propulsion of a Flapping and Oscillating Airfoil. Tech. rep. NACA 567 Langley
Field, VA, United States.
Gendrich, C. P. & M. M. Koochesfahani. 1996. A spatial correlation technique for estimating
velocity fields using molecular tagging velocimetry (MTV). Experiments in Fluids 22(1). 67–
77.
Gendrich, C. P., M. M. Koochesfahani & D. G. Nocera. 1997. Molecular tagging velocimetry and
other novel applications of a new phosphorescent supramolecule. Experiments in Fluids 23(5).
361–372.
Godoy-Diana, Ramiro, Jean-Luc Aider & JosĂŠ Eduardo Wesfreid. 2008. Transitions in the wake of
a flapping foil. Physical Review E 77(1). 016308.
Godoy-Diana, Ramiro, Catherine Marais, Jean-Luc Aider & JosĂŠ Eduardo Wesfreid. 2009. A
model for the symmetry breaking of the reverse BénardâĂŞvon Kármán vortex street produced
by a flapping foil. Journal of Fluid Mechanics 622. 23.
Gursul, I., D.J. Cleaver & Z. Wang. 2014. Control of low Reynolds number flows by means of
fluid-structure interactions. Progress in Aerospace Sciences 64. 17–55.
Hammer, Patrick R. 2016. COMPUTATIONAL STUDY OF THE EFFECT OF REYNOLDS NUMBER AND MOTION TRAJECTORY ASYMMETRY ON THE AERODYNAMICS OF A PITCHING
AIRFOIL AT LOW REYNOLDS NUMBER: Michigan State University dissertation.
Hammer, Patrick R. 2017. personal communication.
Heathcote, Sam & Ismet Gursul. 2007. Flexible Flapping Airfoil Propulsion at Low Reynolds
Numbers. AIAA Journal 45(5). 1066–1079.
Hussain, A. K. M. F. & W. C. Reynolds. 1970. The mechanics of an organized wave in turbulent
shear flow. Journal of Fluid Mechanics 41(02). 241.
James, E. C. 1973. A LINEARIZED THEORY FOR THE UNSTEADY MOTIONS. Tech. Rep.
August Naval Ship R. & D. Center Report 4098 Bethda, Maryland.

140

Jones, K. D., C. M. Dohring & M. F. Platzer. 1998. Experimental and computational investigation
of the Knoller-Betz effect. AIAA Journal 36(7). 1240–1246.
Kang, C.-K., H. Aono, C. E. S. Cesnik & W. Shyy. 2011. Effects of flexibility on the aerodynamic
performance of flapping wings. Journal of Fluid Mechanics 689. 32–74.
von KĂĄrmĂĄn, Th. & J. M. Burgers. 1935. General Aerodynamic Theory-Perfect Fluids. In
William Frederick Durand (ed.), Aerodynamic theory, Springer Berlin.
von KĂĄrmĂĄn, Th. & W. R. Sears. 1938. Airfoil Theory for Non-Uniform Motion. Journal of the
Aeronautical Sciences (Institute of the Aeronautical Sciences) 5(10). 379–390.
Katz, J. & D. Weihs. 1978. Hydrodynamic propulsion by large amplitude oscillation of an airfoil
with chordwise flexibility. Journal of Fluid Mechanics 88(03). 485.
Katz, J. & D. Weihs. 1979. Large amplitude unsteady motion of a flexible slender propulsor.
Journal of Fluid Mechanics 90(04). 713.
Katzmayr, R. 1922. Effect of Periodic Changes of Angle of Attack on Behavior of Airfoils. Tech.
rep. NACA TM-147.
Knoller, R. 1909. Die gesetze des luftwiderstandes. Flug Motortech 3(21). 1–7.
Koochesfahani, Manoochehr M. 1989. Vortical patterns in the wake of an oscillating airfoil. AIAA
Journal 27(9). 1200–1205.
Koochesfahani, Manoochehr M. & Daniel G. Nocera. 2007. Molecular Tagging Velocimetry.
In C. Tropea, AL Yarin & JF Foss (eds.), Handbook of experimental fluid mechanics, chap.
Molecular Tagging Velocimetry. Springer-Verlag Berlin Heidelberg.
Lai, J. C. S. & M. F. Platzer. 1999. Jet Characteristics of a Plunging Airfoil. AIAA Journal 37(12).
1529–1537.
Lighthill, M J. 1969. HYDROMECHANICS OF AQUATIC ANIMAL PROPULSION. Annual
Review of Fluid Mechanics 1. 413–446.
Lighthill, MJ. 1963. Introduction. Boundary Layer Theory. In Rosenhead L (ed.), Laminar
boundary layers, chap. Introducti. New York, NY: Oxford University Press.
Lighthill, Sir James. 1975. Mathematical Biofluiddynamics. Philadelphia: Society for Industrial
and Applied Mathematics.
Lissaman, P B S. 1983. Low-Reynolds-Number Airfoils. Annual Review of Fluid Mechanics 15(1).
223–239.
Mackowski, A W & C H K Williamson. 2015. Direct measurement of thrust and efficiency of an
airfoil undergoing pure pitching. Journal of Fluid Mechanics 524–543.
Marais, C., B. Thiria, J. E. Wesfreid & R. Godoy-Diana. 2012. Stabilizing effect of flexibility in
the wake of a flapping foil. Journal of Fluid Mechanics 710. 659–669.

141

Mitchell, M. R., R. E. Link, A. W. Mix & A. J. Giacomin. 2011. Standardized Polymer Durometry.
Journal of Testing and Evaluation 39(4). 103205.
Monnier, B., A. M. Naguib & M. M. Koochesfahani. 2013. Investigation of the Wake Vortex Pattern
of Rigid and Flexible Airfoils Undergoing Harmonic Pitch Oscillation. In 51st aiaa aerospace
sciences meeting January, 1–14. AIAA 2013-0841.
Monnier, B., A. M. Naguib & M. M. Koochesfahani. 2015. Influence of structural flexibility on
the wake vortex pattern of airfoils undergoing harmonic pitch oscillation. Experiments in Fluids
56(4).
Moore, M. Nicholas J. 2014. Analytical results on the role of flexibility in flapping propulsion.
Journal of Fluid Mechanics 757. 599–612.
Moore, M. Nicholas J. 2015. Torsional spring is the optimal flexibility arrangement for thrust
production of a flapping wing. Physics of Fluids 27(9).
Murray, M.M. 2000. Hydroelasticity Modeling of Flexible Propulsors Durham: Duke University
dissertation.
Naguib, A. M., J. Vitek & M. M. Koochesfahani. 2011. Finite-Core Vortex Array Model of the
Wake of a Periodically Pitching Airfoil. AIAA Journal 49(7). 1542–1550.
Olson, David Arthur. 2011. FACILITY AND FLOW DEPENDENCE ISSUES INFLUENCING THE
EXPERIMENTAL CHARACTERIZATION OF A LAMINAR SEPARATION BUBBLE AT LOW
REYNOLDS NUMBER: Michigan State University dissertation.
Platzer, Max, Kerrin Neace & Chung-Kiang Pang. 1993. Aerodynamic analysis of flapping wing
propulsion. In 31st aerospace sciences meeting Aerospace Sciences Meetings, American Institute
of Aeronautics and Astronautics.
Prempraneerach, P, F S Hover & Michael S Triantafyllou. 2003. The effect of chordwise flexibility
on the thrust and efficiency of a flapping foil. In International symposium on unmanned untethered
submersible technology, .
Ramamurti, Ravi & William Sandberg. 2001. Simulation of Flow About Flapping Airfoils Using
Finite Element Incompressible Flow Solver. AIAA Journal 39(2). 253–260.
Ramananarivo, Sophie, Ramiro Godoy-Diana & Benjamin Thiria. 2011. Rather than resonance,
flapping wing flyers may play on aerodynamics to improve performance. Proceedings of the
National Academy of Sciences of the United States of America 108(15). 5964–9.
Shapiro, Linda G. & George Stockman. 2001. Computer Vision. Upper Saddle River, NJ, USA:
Prentice Hall PTR 1st edn.
Shyy, W., H. Aono, S. K. Chimakurthi, P. Trizila, C. K. Kang, C. E S Cesnik & H. Liu. 2010. Recent
progress in flapping wing aerodynamics and aeroelasticity. Progress in Aerospace Sciences 46(7).
284–327.

142

Spagnolie, Saverio E., Lionel Moret, Michael J. Shelley & Jun Zhang. 2010. Surprising behaviors
in flapping locomotion with passive pitching. Physics of Fluids 22(4). 1–20.
Spedding, G R & P B S Lissaman. 1998. Technical aspects of microscale flight systems. Journal
of Avian Biology 29(4). 458–468.
Stratasys. 2017. PolyJet Materials Data Sheet.
Theodorsen, Theodore. 1934. General Theory of Aerodynamic Instability and the Mechanism of
Flutter. Tech. rep. NACA TR 496.
Thiria, Benjamin & Ramiro Godoy-Diana. 2010. How wing compliance drives the efficiency of
self-propelled flapping flyers. Physical Review E 82(1). 015303.
Triantafyllou, GS. 1993. Optimal thrust development in oscillating foils with application to fish
propulsion. Journal of Fluids and Structures 7(2). 205–224.
Triantafyllou, M. S, G S Triantafyllou & R Gopalkrishnan. 1991. Wake mechanics for thrust
generation in oscillating foils. Physics of Fluids A: Fluid Dynamics 3(12). 2835.
Wei, Z. & Z.C. Zheng. 2014. Mechanisms of wake deflection angle change behind a heaving airfoil.
Journal of Fluids and Structures 48. 1–13.
Wu, Jie-Zhi & Jain-Ming Wu. 1993. Interactions between a solid surface and a viscous compressible
flow field. Journal of Fluid Mechanics 254. 183–211.
Wu, T. Yao-Tsu. 1961. Swimming of a waving plate. Journal of Fluid Mechanics 10(03). 321.
Young, John & Joseph C. S. Lai. 2004. Oscillation Frequency and Amplitude Effects on the Wake
of a Plunging Airfoil. AIAA Journal 42(10). 2042–2052.
Zhu, Qiang. 2007. Numerical Simulation of a Flapping Foil with Chordwise or Spanwise Flexibility.
AIAA Journal 45(10). 2448–2457.

143