THEORETICAL AND COMPUTATIONAL ANALYSIS OF UNSTEADY FLOW
OVER AN ASYMMETRIC AIRFOIL
By
Shiwei Qin

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirement
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2012

ABSTRACT

THEORETICAL AND COMPUTATIONAL ANALYSIS OF UNSTEADY FLOW
OVER AN ASYMMETRIC AIRFOIL
By
Shiwei Qin

Various methods are employed for studying steady and unsteady flows over an asymmetric
SD7003 airfoil. First, linear stability methods are used for the investigation of airfoil flow
stability. Second, flows over the stationary airfoil with steady and unsteady freestream conditions
are computed by the large eddy simulation (LES) method.
The laminar separation bubble (LSB) region of the airfoil flow at Reynolds number of
60,000 and angle of attack (AoA) of 4°is analyzed by the local and global stability methods.
Both methods correctly predict that the flow is unstable and transitions to turbulence at these
conditions. The maximum growth rates predicted by the two methods are in agreement.
The LES results for several steady freestream flows are shown to be consistent with the
available experimental and numerical data, and then three types of unsteady freestream flows are
simulated by the LES method. In the first type, the AoA is fixed at 4° while the freestream
,
velocity magnitude varies harmonically with a mean Reynolds number of 60,000, reduced
frequency from π/8 to 2π, and amplitudes of 0.183 and 0.366 of the mean velocity. In the second
type of unsteady flows, the freestream velocity magnitude is fixed at Reynolds number of 60,000,
but the AoA varies harmonically with a mean value of 4° reduced frequency from π/8 to 2π, and
,
amplitudes of 4°and 8° In the third type of unsteady flow, a wind gust model is employed.
.
For the flows with oscillating freestream velocity magnitude, the effect of freestream
unsteadiness on the aerodynamic forces, vorticity field, velocity fluctuations, and boundary layer

separation and reattachment are studied. The vorticity plots show that the size of the recirculation
region changes slightly during a freestream cycle. The mean lift and drag coefficients are nearly
the same as those obtained for the steady mean freestream flow. However, there is a phase shift
between the aerodynamic forces and the freestream velocity which is mainly affected by the
reduced frequency. Higher reduced frequencies and freestream amplitudes cause “wider”
oscillations in the lift and drag forces, separation and reattachment locations. For the case with
the highest frequency and the largest amplitude, the spanwise velocity fluctuations are much
lower than those of the streamwise and normal velocities, while the fluctuation in the streamwise
velocity is slightly more than that of the normal velocity.
For the flows with oscillating freestream AoA, higher reduced frequencies and AoA
amplitudes cause wider oscillations in the lift and drag forces. The amplitude of the lift force is,
however, about one order of magnitude larger than that of the drag force. Compared to the steady
mean freestream flow, there is little change in the mean lift, while the mean drag is reduced due
to Katzmayr effect. The size of the recirculation region changes significantly during a freestream
cycle at low reduced frequencies and high AoA amplitudes. Higher reduced frequency decreases
the amplitude of separation point oscillations, but higher AoA amplitude increases it. Significant
oscillation in the surface friction occurs near the trailing edge at high reduced frequencies and
AoA amplitudes, resulting from strong vortex shedding at the trailing edge. Compared to the
steady mean freestream flow, the mean separation point moves downstream and the mean
reattachment point moves upstream when the freestream velocity magnitude or AoA oscillates.
The wind gust simulation shows a rapid increase in the aerodynamic loads on the airfoil at
the beginning of the gust. It is also shown that small fluctuations in the gust can cause significant
fluctuations in the aerodynamic forces.

Copyright by
SHIWEI QIN
2012

ACKNOWLEDGEMENTS

I would like to thank my advisors Profs. Farhad Jaberi and Mei Zhuang for their endless
support, advice and patience during my PhD studies. I greatly appreciate their faith in allowing
me to work on this project, which has been very challenging and rewarding. I would also like to
thank other members of my PhD committee, Profs. Alejandro Diaz, Zhengfang Zhou, and
Manoochehr Koochesfahani. I am sincerely grateful for their very helpful suggestions.
Additionally, I like to thank very much the Air Force Research Laboratory at WrightPatterson Air Force Base and Drs. Jack Benek and Miguel Visbal for supporting this project and
for their valuable comments and suggestions. I also want to thank Marshall Galbraith for his
assistance.
Finally, I would like to thank my friends and especially my family for their support during
my PhD studies at Michigan State University. Special thanks to my wife, Yanming.

v

TABLE OF CONTENTS

LIST OF TABLES ....................................................................................................................... viii
LIST OF FIGURES ....................................................................................................................... ix
LIST OF SYMBOLS ................................................................................................................. xviii
1.

Introduction .............................................................................................................................. 1

2.

Linear Stability Analysis ......................................................................................................... 6
2.1 Global Stability Formulation .......................................................................................... 6
2.1.1 Base flow .............................................................................................................. 7
2.1.2 Perturbation equations .......................................................................................... 7
2.2 Local Stability Formulation .......................................................................................... 10
2.3 Base Flow Results......................................................................................................... 11
2.4 Local Stability Results .................................................................................................. 12
2.5 Global Stability Results ................................................................................................ 13

3.

Simulation Method ................................................................................................................ 20
3.1 Background ................................................................................................................... 20
3.2 Filtered Equations ......................................................................................................... 21
3.3 Subgrid-Scale Model .................................................................................................... 22
3.4 Grid and Boundary Conditions ..................................................................................... 24
3.5 Separation, Transition, Reattachment, and Force Coefficients .................................... 27
3.6 Comparison of 2D and 3D Simulations ........................................................................ 29

4.

Large Eddy Simulations of Steady Freestream Flows over SD7003 Airfoil ......................... 41

5.

Large Eddy Simulation of Unsteady Freestream Flows over SD7003 Airfoil ...................... 50
5.1 Flows with Oscillating Freestream Velocity Magnitude .............................................. 50
5.1.1 Aerodynamic forces............................................................................................ 51
5.1.2 Vorticity field ..................................................................................................... 60
5.1.3 Velocity fluctuations .......................................................................................... 71
5.1.4 Boundary layer separation and reattachment ..................................................... 76
5.2 Flows with Oscillating Angle of Attack (AoA) ............................................................ 80
5.2.1 Aerodynamic forces............................................................................................ 82
5.2.2 Vorticity field ..................................................................................................... 92
5.2.3 Velocity fluctuations ........................................................................................ 100
5.2.4 Boundary layer separation and reattachment ................................................... 104
5.3 A Discussion Regarding Effect of Freestream Unsteadiness on Mean Separation and
Reattachment Points ................................................................................................... 108
5.4 Flow with Wind Gust.................................................................................................. 112

6.

Summary and Conclusions .................................................................................................. 118

APPENDICES ............................................................................................................................ 122
vi

APPENDIX A ...................................................................................................................... 123
APPENDIX B ...................................................................................................................... 125
APPENDIX C ...................................................................................................................... 129
BIBLIOGRAPHY ....................................................................................................................... 133

vii

LIST OF TABLES

Table 1: Mean quantities obtained from 2D and 3D solutions for the flow with steady
freestream and Re=60,000 and AoA=4°...........................................................................30
.
Table 2: Mean quantities obtained from 2D and 3D solutions for the flow with steady
freestream and Re=60,000 and AoA=8°.............................................................................32
.
Table 3: Computed and measured time-averaged data for the flow with steady freestream,
Re=60,000 and AoA=4° LES refers to our simulations with WALE SGS model. ............42
.
Table 4: Computed time-averaged data for the flow with steady freestream, Re=60,000
and AoA=8° ........................................................................................................................45
.
Table 5: Time-averaged data for the flow with steady freestream, Re=60,000, and
AoA=4° 8° 12° 16° ..........................................................................................................47
, ,
,
.
Table 6: Separation and reattachment points obtained by LES and TMUAL experiments...........48
Table 7: Mean quantities for cases with oscillating velocity magnitude. ......................................51
Table 8: Time-averaged data for cases with oscillating AoA. .......................................................81
Table 9: Measured separation and reattachment points for AoA=4°and Re=60,000 from
different experimental facilities. ........................................................................................108
Table 10: Constants of the wind gust model. ...............................................................................113

viii

LIST OF FIGURES

Figure 1: The SD7003 airfoil. ..........................................................................................................5
Figure 2: Base flow around a SD7003 airfoil, Re=60,000, α=4° [7]. Contours represent
time-averaged x-velocity. The rectangular solid line shows the domain used in the
global stability analysis (For interpretation of the references to color in this and all other
figures, the reader is referred to the electronic version of this dissertation)............................12
Figure 3: Time-averaged x-velocity (U) profile at x/C=0.5. .........................................................13
Figure 4: Local stability analysis results for time-averaged x-velocity (U) profile at x/C=0.5;
(left) growth rate vs. wavenumber, (right) angular frequency vs. wavenumber. ....................13
Figure 5: The four most unstable modes at wavenumber k=20, obtained with different base
flow grids. ................................................................................................................................14
Figure 6: The most unstable modes vs. wavenumber; (upper) maximum growth rate, (lower)
angular frequency. ...................................................................................................................16
Figure 7: Turbulent-kinetic-energy frequency spectra at x/C=0.5 (about center of LSB) [7]. ......16
Figure 8: Reconstructed total x-velocity using the conjugate pair of most unstable modes for
k=1; (a) t=0, (b) t=T/4, (c) t=T/2, (d) t=3T/4...........................................................................18
Figure 9: Reconstructed velocity disturbances using the conjugate pair of most unstable
modes for k=1 at t=T/2. ...........................................................................................................19
Figure 10: Instantaneous normalized z-velocity w U m on a z-section, α=4°, Re=60,000. ..........19
Figure 11: Comparison of time- and spanwise-averaged pressure coefficient over the
SD7003 airfoil, obtained by our LES with different SGS models and compared with
those obtained by Galbraith et al. [7]. The flow Reynolds number based on freestream
velocity and airfoil chord is 60,000, and AoA=4° KET stands for Kinetic Energy
.
Transportation model. ..............................................................................................................24
Figure 12: Computational grid: (a) grid on XY plane, (b) grid in the vicinity of airfoil, (c)
grid near the trailing edge, (d) 3D airfoil. ................................................................................25
Figure 13: Schematic of flow over airfoil, showing the freestream and the coordinate system.
.................................................................................................................................................28
Figure 14: Time-averaged surface pressure coefficient obtained by 2D and 3D simulations
of the flow with steady freestream and Re=60,000 and AoA=4°...........................................30
.

ix

Figure 15: Time-averaged skin friction coefficient on the upper surface of the airfoil
obtained by 2D and 3D simulations of the flow with steady freestream and Re=60,000
and AoA=4°............................................................................................................................31
.
Figure 16: Instantaneous normalized z-vorticity ( z C U m ) contours obtained by 2D and
3D simulations of the flow with steady freestream and Re=60,000 and AoA=4° .................31
.
Figure 17: Time-averaged surface pressure coefficient obtained by 2D and 3D simulations
of the flow with steady freestream and Re=60,000 and AoA=8°...........................................32
.
Figure 18: Time-averaged skin friction coefficient on the upper surface of the airfoil
obtained by 2D and 3D simulations of the flow with steady freestream and Re=60,000
and AoA=8°............................................................................................................................33
.
Figure 19: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=10,000. .............................................................................................36
Figure 20: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations
of the flow with oscillating AoA and Re=10,000 (colored by u/Um for better visibility). ......36
Figure 21: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=30,000. .............................................................................................37
Figure 22: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations
of the flow with oscillating AoA and Re=30,000. ...................................................................37
Figure 23: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=60,000. .............................................................................................38
Figure 24: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations
of the flow with oscillating AoA and Re=60,000. ...................................................................38
Figure 25: Pressure coefficient obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=60,000. .............................................................................................40
Figure 26: Instantaneous contours of x-vorticity ( x C U m ) obtained by 3D simulations of
the flow with steady freestream and Re=60,000, AoA=16°...................................................40
.
Figure 27: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for the flow with
steady freestream, Re=60,000, AoA=4° 8° 12°and 16° ......................................................42
, ,
.
Figure 28: Time-averaged surface pressure coefficient for the flow with steady freestream,
Re=60,000 and AoA=4° .........................................................................................................44
.
Figure 29: Time-averaged skin friction coefficient on the upper airfoil surface for the flow
with steady freestream, Re=60,000 and AoA=4°...................................................................44
.

x

Figure 30: Time-averaged surface pressure coefficient for the flow with steady freestream,
Re=60,000 and AoA=8° .........................................................................................................45
.
Figure 31: Time-averaged skin friction coefficient on the airfoil upper surface for the flow
with steady freestream, Re=60,000 and AoA=8°...................................................................46
.
Figure 32: Time-averaged surface pressure coefficient for the flow with steady freestream at
different AoAs and Re=60,000. ...............................................................................................47
Figure 33: Time-averaged lift and drag coefficients obtained by LES for flows with steady
freestream condition. ...............................................................................................................49
Figure 34: Time variations CL for flows with oscillating freestream speed (cases U1-5). ............53
Figure 35: Time variations of CD for flows with oscillating freestream speed (cases U1-5). .......53
Figure 36: Time variations of CL for flows with oscillating freestream speed (cases U6-10).
The phase shift for one of the cases is shown as an example. .................................................54
Figure 37: Time variations of CD for flows with oscillating freestream speed (cases U6-10). .....54
Figure 38: Time variations of CL* for flows with oscillating freestream speed (cases U6-10). ....55
Figure 39: Time variations of CD* for flows with oscillating freestream speed (cases U6-10). ....55
Figure 40: Amplitude of CL and CD for flows with oscillating freestream velocity magnitude.
.................................................................................................................................................56
Figure 41: Phase shift between force coefficients and freestream velocity for flows with
oscillating freestream velocity magnitude. ..............................................................................57
Figure 42: Pressure coefficient on the airfoil surface at the time of maximum acceleration of
freestream velocity for case U6 with k  2 and   0.366 and case U7 with k  
and   0.366 ..........................................................................................................................59
Figure 43: Pressure coefficient around the airfoil at the time of maximum acceleration of
freestream velocity for case U6 with k  2 and   0.366 . ................................................59
Figure 44: The n-s coordinate system around the airfoil surface. ..................................................62
Figure 45: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one residence
time for the flow with steady freestream, Re=60,000, AoA=4° (T is the residence time
.
C U m , same for Figures 46-47). ............................................................................................63
Figure 46: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one
residence time for the flow with steady freestream, Re=60,000, AoA=4°.............................63
.

xi

Figure 47: Integrals of spanwise-averaged vorticity magnitude |Ω| for the flow with steady
freestream, Re=60,000, AoA=4° (a) Integrated along n axis; (b) integrated along n and
.
then s axes. ...............................................................................................................................64
Figure 48: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for
case U1 with k  2 and   0.183 . (t0 is the beginning of the cycle when U equals
Um and is increasing. T is the period. Same for Figures 49-60) .............................................64
Figure 49: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one
cycle for case U1 with k  2 and   0.183 . .......................................................................65
Figure 50: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U1 with k  2
and   0.183 . (a) Integrated along n axis; (b) integrated along n and then s axes. ...............65
Figure 51: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for
case U5 with k   8 and   0.183 . .....................................................................................66
Figure 52: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one

cycle for case U5 with k   8 and   0.183 . ......................................................................66

Figure 53: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U5 with k   8
and   0.183 . (a) Integrated along n axis; (b) integrated along n and then s axes. ...............67
Figure 54: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for
case U6 with k  2 and   0.366 . ......................................................................................67
Figure 55: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one
cycle for case U6 with k  2 and   0.366 . .......................................................................68
Figure 56: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U6 with k  2
and   0.366 . (a) Integrated along n axis; (b) integrated along n and then s axes. ...............68
Figure 57: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for
case U10 with k   8 and   0.366 ....................................................................................69
Figure 58: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one

cycle for case U10 with k   8 and   0.366 . ...................................................................69

Figure 59: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U10 with
k   8 and   0.366 . (a) Integrated along n axis; (b) integrated along n and then s
axes. .........................................................................................................................................70
Figure 60: Integral of spanwise-averaged vorticity magnitude |Ω| along n axis for cases U1,
U5, U6, and U10. .....................................................................................................................71

xii

Figure 61: Integrated values of velocity RMS along n axis for the flow with steady
freestream, Re=60,000 and AoA=4° ......................................................................................73
.
Figure 62: Integrated values of velocity RMS along n axis for case U6 with k  2 and
  0.366 . ...............................................................................................................................74
Figure 63: Integrated values of velocity RMS along n and s axes for case U6 with k  2
and   0.366 . .........................................................................................................................75
Figure 64: Integrated values of turbulence kinetic energy (TKE) along n and s axes for case
U6 with k  2 and   0.366 . ..............................................................................................75
Figure 65: Comparison of integrated values of velocity RMS along n axis between case U6
with k  2 and   0.366 and the flow with steady freestream with Re=60,000 and
AoA=4° (a) u RMS; (b) v RMS; (c) w RMS..........................................................................76
.
Figure 66: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for the flow with steady freestream, Re=60,000 and AoA=4°
....................................78
Figure 67: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case U5 with k   8 and   0.183 . ...................................................................79
Figure 68: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case U1 with k  2 and   0.183 . ....................................................................79
Figure 69: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case U6 with k  2 and   0.366 . ....................................................................80
Figure 70: Time variations of CL for flows with oscillating AoA, cases A1-5. ............................84
Figure 71: Time variations of CD for flows with oscillating AoA, cases A1-5. ............................85
Figure 72: Time variations of CL for flows with oscillating AoA, cases A6-10. ..........................85
Figure 73: Time variations of CD for flows with oscillating AoA, cases A6-10. ..........................86
Figure 74: Time variations of CL* for flows with oscillating AoA, cases A6-10. .........................86
Figure 75: Time variations of CD* for flows with oscillating AoA, cases A6-10. .........................87
Figure 76: Amplitude of oscillations in CL and CD for flows with oscillating freestream AoA
and the prediction by Theodorsen’s theory for a pitching airfoil in the steady mean
freestream. ...............................................................................................................................87
Figure 77: Contours of pressure coefficient around the airfoil at the time of maximum
acceleration of freestream AoA for case A1 with k  2 and  f  4 . ................................88

xiii

Figure 78: Pressure coefficient on the airfoil surface at the time of maximum acceleration of
freestream AoA for case A1 with k  2 and  f  4 ..........................................................89
Figure 79: Phase shift of CL and CD for flows with oscillating AoA and the prediction by
Theodorsen’s theory for a pitching airfoil in the steady mean freestream. The two
Theodorsen curves overlay each other. ...................................................................................91
Figure 80: Computed time-averaged lift coefficient for flows with oscillating freestream
AoA and the prediction by Theodorsen’s theory for a pitching airfoil in the steady mean
freestream. The two Theodorsen curves overlay each other. ..................................................91
Figure 81: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for flows with steady
freestream, Re=60,000 and AoA=4° 8°and 12° ...................................................................94
,
.
Figure 82: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for flows
with steady freestream, Re=60,000 and AoA=4° 8°and 12°................................................94
,
.
Figure 83: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A1 with

k  2 and  f  4 . At t0, α(t) = 4°; at t0+T/4, α(t) = 8°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = 0°. ....................................................................................................................95
Figure 84: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A1
with k  2 and  f  4 . At t0, α(t) = 4°; at t0+T/4, α(t) = 8°; at t0+T/2, α(t) = 4° at
;
t0+3T/4, α(t) = 0°. ....................................................................................................................95
Figure 85: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A1 with k  2
and  f  4 . (a) Integrated along n axis; (b) integrated along n and then s axes. ..................96
Figure 86: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A5 with

k   8 and  f  4 . At t0, α(t) = 4°; at t0+T/4, α(t) = 8°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = 0°. ....................................................................................................................96
Figure 87: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A5
with k   8 and  f  4 . At t0, α(t) = 4°; at t0+T/4, α(t) = 8°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = 0°. ....................................................................................................................97
Figure 88: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A5 with k   8
and  f  4 . (a) Integrated along n axis; (b) integrated along n and then s axes. ..................97

xiv

Figure 89: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A6 with

k  2 and  f  8 . At t0, α(t) = 4°; at t0+T/4, α(t) = 12°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = -4° ..................................................................................................................98
.
Figure 90: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A6
with k  2 and  f  8 . At t0, α(t) = 4°; at t0+T/4, α(t) = 12°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = -4° ..................................................................................................................98
.
Figure 91: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A6 with k  2
and  f  8 . (a) Integrated along n axis; (b) integrated along n and then s axes. ..................99
Figure 92: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A10 with
k   8 and  f  8 . At t0, α(t) = 4°; at t0+T/4, α(t) = 12°; at t0+T/2, α(t) = 4°; at
t0+3T/4, α(t) = -4° ..................................................................................................................99
.
Figure 93: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case
A10 with k   8 and  f  8 . At t0, α(t) = 4°; at t0+T/4, α(t) = 12°; at t0+T/2, α(t) =
4°; at t0+3T/4, α(t) = -4°.......................................................................................................100
.
Figure 94: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A10 with

k   8 and  f  8 . (a) Integrated along n axis; (b) integrated along n and then s axes. ..100
Figure 95: Integral of velocity RMS along airfoil surface normal direction, n for case A6
with k  2 and  f  8 . .....................................................................................................102
Figure 96: Integral of velocity RMS along n and s axes for case A6 with k  2 and

 f  8 . .................................................................................................................................103
Figure 97: Integral of turbulence kinetic energy (TKE) along n and s axes for case A6 with

k  2 and  f  8 . .............................................................................................................103
Figure 98: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A5 with k   8 and  f  4 ....................................................................105
Figure 99: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A3 with k   2 and  f  4 . ..................................................................106
Figure 100: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A1 with k  2 and  f  4 . ....................................................................106
xv

Figure 101: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A6 with k  2 and  f  8 . ....................................................................107
Figure 102: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A10 with k   8 and  f  8 . .................................................................107
Figure 103: Comparison of the separation and reattachment points for freestream turbulence
intensity of 1.0% and 1.5% at Re=20,000 [29, 52]. ..............................................................110
Figure 104: Measured wind speed at (upper) a sea mast in Vindeby offshore wind farm [56]
and (lower) a 10 meter tower [57]. ........................................................................................112
Figure 105: Time variations of the wind gust model parameters, velocity and lift and drag
coefficients induced by the wind gust. CL and CD are computed by normalizing the lift
and drag forces with U0 and U(t), both. .................................................................................116
Figure 106: Instantaneous iso-surfaces of z-vorticity ( z C U0  10 ) simulated by LES
with wind gust at different times. The numbers on figures match the times in Figure 105. .117
Figure 107: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U2 with
k   and   0.183 . ............................................................................................................125
Figure 108: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U3 with
k   2 and   0.183 . ........................................................................................................126
Figure 109: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U4 with
k   4 and   0.183 . ........................................................................................................126
Figure 110: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U7 with
k   and   0.366 . ............................................................................................................127
Figure 111: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U8 with

k   2 and   0.366 . ........................................................................................................127

Figure 112: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U9 with
k   4 and   0.366 . ........................................................................................................128
Figure 113: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A2 with

k   and  f  4 . ...............................................................................................................129

xvi

Figure 114: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A3 with

k   2 and  f  4 . ...........................................................................................................130
Figure 115: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A4 with
k   4 and  f  4 . ...........................................................................................................130
Figure 116: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A7 with
k   and  f  8 . ...............................................................................................................131
Figure 117: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A8 with
k   2 and  f  8 .............................................................................................................131
Figure 118: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A9 with

k   4 and  f  8 .............................................................................................................132

xvii

LIST OF SYMBOLS

Symbols
A, B

Assembly mass matrices

C

Acoustic speed in fluid, or Chord length

CD

Drag coefficient

CF

Skin friction coefficient

C ij

Cross stress tensor

CL

Lift coefficient

CP

pressure coefficient

Cw

WALE model constant

D

Drag force per unit span of the airfoil

D1  D7

Constants in wind gust model

G

Filter function

G1  G3

Functions in wind gust model

I

Unit matrix

L

Lift force per unit span of the airfoil

Ma

Mach number

Re

Reynolds number

Rij

Reynolds stress tensor

S

Acoustic speed in the fluid

S ij

Filtered rate-of-strain tensor

xviii

d
Sij

Tensor in WALE LES model

T

Period

TR

Residence time

U

Freestream velocity in x direction

V

Freestream velocity in y direction, or Volume of a computational cell

W

Freestream velocity in z direction

V0

Reference velocity

XS

Separation location

XR

Reattachment location

XT

Transition location

f

Fluid body force (e.g. gravitation)

g ij

Filtered velocity gradient

i

Imaginary unit

k

Wavenumber or reduced frequency

n

Boundary surface normal

p

Pressure

s

axis tangent to the airfoil surface

t

Time

u, v, w

Velocity components in x, y, and z directions

v

Velocity vector

x, y, z

Cartesian coordinates



Filter width
xix



Infinity

Greek Symbols


Vector of eigenfunctions

|Ω|

Vorticity Magnitude



Angle of attack

 ij

Kronecker symbol



Von Karman constant



Kinematic viscosity

t

Turbulent viscosity



Density

0

Reference density

 ij

Subgrid-scale stress tensor



Angular frequency or vorticity

Subscripts
f

Fluctuation quantities

i

Imaginary part of complex number

i, j , k

Coordinate index

m

Mean quantities

r

Real part of complex number

t

Turbulent quantities

*

Two-dimensional operator

xx

Superscripts



Averaged or filtered quantities

^

Amplitude of disturbance quantities

'

Disturbance quantities or relative coordinate vector

'', ''''

Second- and fourth-order derivatives

Operators
2

Laplace Operator

xxi

1.

Introduction

There has been a growing interest in small Unmanned Air Vehicles, including Micro Air
Vehicles (MAVs), in recent years due to their growing capability to perform a wide range of
missions. These vehicles typically operate at low to moderate Reynolds numbers in the order of
4

5

10 to 10 due to low speed and small vehicle size. At this range of Reynolds numbers, the flow
near the airfoil leading edge is still laminar, and flow separation occurs at some point due to
substantial adverse pressure gradients. At not very high Reynolds number and AoA, the flow can
reattach to the airfoil. The flow streamlines show a laminar separation bubble (LSB) [1-3]
bounded by the separation and reattachment points. Based on its size, a LSB can be categorized
as either a small bubble or a large bubble. A small bubble covers a small portion of the airfoil
surface and does not have a significant effect on the velocity and pressure distributions. A large
bubble covers a considerable portion of the surface and significantly changes the pressure
distribution and the velocity field. A separation bubble, especially a large one, decreases the lift
and increases the drag, and consequently reduces the efficiency of the airfoil. The onset and
successive breakdown of the LSB at low Reynolds number flows is known to be detrimental to
the performance, endurance, and stability of MAVs.
To better understand the stability characteristics of the flow over a 2D wing at low
Reynolds number, local and global stability methods are used in this study to analyze the flow
instability over the SD7003 airfoil with a considerable LSB. Stability analysis is considered to be
useful for better understanding of the instability of attached and separated flows [4]. Hammond
and Redekopp [5] used a modeled separation bubble with the assumption of quasi-parallel flow
and performed linear local stability analysis with a one-dimensional eigenvalue and OrrSommerfeld solver to determine the conditions for absolute instability of the flow over an airfoil

1

with LSB. Theofilis et al. [6], on the other hand, used non-parallel two-dimensional steady flows
as base flows in their “BiGlobal” linear stability analyses. Since stability analysis is sensitive to
the base flow, it is important that the base flow is carefully quantified. Thus, experimental
measurements or “reliable” numerical simulations are necessary for doing an accurate stability
analysis. In the current study, the base flow is computed by LES [7]. Results of local and global
stability analyses are compared, and stability features of the flow over the SD7003 airfoil are
discussed.
Many studies have been done on the characterization and control of LSB over symmetric
and asymmetric airfoils. A commonly studied asymmetric airfoil for low-moderate Reynolds
number flows is the SD7003 airfoil. Experimental study and measurements of flows around the
SD7003 airfoil have been reported by a number of researchers using different facilities and
methods. Radespiel et al [8] conducted experiments in a water channel and a low-noise wind
tunnel at Technical University of Braunschweig (TU-BS) and obtained high resolution velocity
and Reynolds stress data. Detailed particle image velocimetry (PIV) measurements were
conducted by Ol et al. [9] at the Air Force Research Laboratory (AFRL) water channel.
Measurements in a water channel using Molecular Tagging Velocimetry (MTV) were conducted
by Katz et al. [10-13] at Turbulent Mixing and Unsteady Aerodynamics Laboratory (TMUAL) at
Michigan State University. Numerical tools were also applied to simulate flows over the SD7003
airfoil. Reynolds-averaged Navier-Stokes (RANS) [14-17] and LES [7, 18-21] models were both
employed to capture the flow dynamics and to predict the aerodynamic forces.
Unsteady aerodynamics of maneuvering airfoils has also been studied by experimental and
numerical methods [22-24]. Both rigid and flexible airfoils were considered. The airfoil motions
studied include pitching, plunging, or combination of the two. Bohl and Koochesfahani [22, 25]

2

reported measurements of the vortical field in the wake of a pitching airfoil. LES computation of
flows past a plunging airfoil was reported by Visbal [23]. A comparison between measurement
and computation of an airfoil in pitching and plunging motions was given by McGowan et al.
[24]. Inviscid theory for this type of flows has a 75 years history and has provided valuable
insight to many questions regarding unsteady flows over moving airfoils. For example,
Theodorsen [26] presented the lift and pitching moment for a two-dimensional, flat-plate airfoil
under harmonic pitching and plunging motions. Von Ká n and Sears [27] developed a general
rmá
theory for unsteady aerodynamics of an airfoil in non-uniform motion. An excellent review of
works in this area is given in reference [28].
Discrepancy of the measured mean separation and reattachment points in flows over the
SD7003 airfoil, obtained from different facilities has been reported [9, 11, 29]. The main reason
for these discrepancies is suggested to be the different levels of freestream turbulence in the
experimental facilities. Experiments by Olson [12] confirmed that the increased freestream
turbulence delays the separation and triggers earlier reattachment.
For MAVs typically flying at low speeds, freestream unsteadiness may cause significant
changes in aerodynamic forces and thus possible flight instability. Some analytical works based
on inviscid theory were developed by Isaacs [30] and Greenberg [31] for unsteady flows over
airfoils. By extending the method of Theodorsen [26], they calculated the lift and moment of an
airfoil at a fixed AoA in a freestream with harmonic velocity magnitude oscillation. Although
different assumptions are applied to the wake, the predicted results by Isaacs and Greenberg are
shown to be very similar. There are not many experimental studies on fixed airfoils operating in
an unsteady freestream flow. This might be partly due to difficulty of controlling the flow in
experiments. The frequency of freestream flow oscillation is limited by the fact that the

3

freestream uniformity across the test section becomes unacceptable even at relatively low
frequencies [32]. Therefore, experiments of this type usually have very low reduced frequencies.
To verify the Katzmayr effect [33], Toussaint et al. [34] used an array of pitching blades in front
of the tested airfoil to create a freestream with oscillating flow direction. In this experiment, the
maximum reduced frequency is less than 0.5, and freestream properties (e.g. level of turbulence)
behind the oscillating blades were not given. Williams et al. [35] measured flows over a semicircular wing in a wind tunnel that generates sinusoidal oscillations in the freestream velocity
magnitude with reduced frequencies less than 0.7. They found that the measured lift force lags
the freestream velocity by a noticeable amount even at low reduced frequencies, and the phase
shift does not show a strong dependence on the mean AoA. The amplitude of oscillations in the
lift coefficient was shown to be affected by both the AoA amplitude and the frequency of
oscillations in freestream AoA. However, the measured amplitude and phase shift in lift
coefficient differed significantly from predictions by Greenberg’s theory. Using k-ω turbulence
model and a 2D mesh, Gharali and Johnson [36] conducted unsteady RANS simulations of flows
over fixed S809 airfoil with oscillating freestream velocity direction. In these simulations, the
reduced frequency ranges from 0.026 to 18. The results for a case with a very low reduced
frequency of 0.026 were compared to experimental results of a pitching airfoil with some level of
agreement. Their study also showed that, at high reduced frequencies, the airfoil experiences
oscillations in lift and drag coefficients with very large amplitude.
In this work, we study the details of the unsteady flow over the SD7003 airfoil by unsteady
numerical simulations with the LES method. The simulated oscillation frequencies extend from a
relatively low value to high values, well beyond the highest found in experiments conducted with
freestream flow oscillations. This provides a better understanding of the unsteady flow physics

4

over two-dimensional airfoils. Three types of unsteady freestream flows are considered. In the
first type, the freestream velocity direction and AoA are fixed, but the freestream velocity
magnitude oscillates harmonically in time. Frequency and amplitude of the oscillation are varied
(within a range possible in our simulation) to investigate the flow response to different levels of
flow unsteadiness. In the second type of unsteady flows considered in this study, the freestream
velocity magnitude is fixed but the freestream velocity direction or AoA is changed harmonically
in time. Again, the flow direction is changed over a range of frequencies and amplitudes. Effects
of these two types of freestream flow oscillations on the aerodynamic forces, vorticity field,
velocity fluctuations, boundary layer separation and reattachment are studied. The computed lift
coefficients are compared to those predicted by the inviscid theory. Finally in the third unsteady
flow type considered in this study, a wind gust model is used for studying the flow over the
airfoil exposed to more realistic “wind type” flow conditions.
The SD7003 airfoil is used in all analyses presented in this dissertation. Figure 1 shows a
schematic 2D view of this airfoil, which has a maximum thickness and camber of 8.5% and
1.48% of the chord length, respectively. This airfoil is chosen because of the availability of
experimental and computational data.

Figure 1: The SD7003 airfoil.

5

2.
2.1

Linear Stability Analysis

Global Stability Formulation
The global stability analysis conducted here is based on the numerical solutions of a three-

dimensional partial-derivative eigenvalue problem, which describes small-amplitude threedimensional disturbances developing on a two-dimensional steady or time-averaged base flow
around the airfoil. The base flow is assumed to be the same in the spanwise direction (denoted by
z) such that /z=0 for the base flow. Following Ding and Kawahara [37], slight compressibility
is assumed so that the singularity in the spectrum of the eigenproblem is avoided. The singularity
arises from the continuity constraint of incompressible flow as explained by Malkus [38]. The
nondimensional scales for length, velocity, time, density, and kinematic pressure are taken to be
C , U 0 , C U 0 , 0 , and SU 0 , where S is the sound speed in the fluid. The dimensionless form

of the continuity equation is
D
   V  0 ,
Dt

where

(1)

D 
   V   is the substantial derivative. With the assumption of small
Dt t

compressibility, the kinematic pressure p (which is the pressure divided by density) is a
function of density as [37],
Dp
1 D

Dt  M a Dt

(2)

where M a  U 0 S is the Mach number of the flow. By substituting Equation (2) into Equation
(1), a modified continuity equation is obtained

6

Dp
1

  V 0 .
Dt M a

(3)

Under small compressibility conditions, considering the kinematic viscosity  of the fluid to be
constant and f to be a constant body force (e.g. gravitational force), the nondimensional
momentum equation can be written as

Dv
1
1

p 
Dt
Ma
Re

 2 1

 v  3     v    f ,



(4)

where Re  U 0 C  is the Reynolds number.
2.1.1

Base flow
For global stability analysis, the airfoil is assumed to be infinitely long in the spanwise

direction, so that the base flow becomes two-dimensional and steady. Compressibility is
neglected in computing the base flow. Therefore, the base flow equations become
V  0

(5)

1 2
 V *  V  *P  R * V ,

(6)

e

where V and P are the velocity and kinematic pressure in the base flow, respectively. *
represents the two-dimensional gradient operator.
2.1.2

Perturbation equations
A flow variable may be decomposed into a steady base flow and an unsteady perturbation.

The perturbed flow variables may be described in the following form

v  x, y, z, t   V( x, y)  v '  x, y, z, t 

(7)

p  x, y, z, t   P ( x, y)  p '  x, y, z, t  .

(8)

7

The perturbations u ', v ', w ', p ' are assumed to be inhomogeneous in x and y, have a harmonic
dependence on the spanwise coordinate z, and grow exponentially in time as,
ˆ
p '  ip  x, y  eikz t ,

(9)

ˆ
u '  iu  x, y  eikz t ,

(10)

ˆ
v '  iv  x, y  eikz t ,

(11)

ˆ
w '  w  x, y  eikz t .

(12)

In Equations (9) to (12), i is the imaginary number, k  2 C Lz is the nondimensional spanwise
wavenumber where Lz is the spanwise wavelength, and   r  ii denotes the complex
growth rate. The reason for using the imaginary amplitude in the normal modes is to avoid
complex arithmetic in the calculation of eigenproblem [37].
Substituting Equations (7) and (8) into the continuity Equation (3) and the momentum
Equation (4), subtracting the base flow Equations (5) and (6), neglecting the high order terms of
perturbations and replacing the perturbed variables with Equations (9) to (12) result in the
following eigenproblem for perturbations with the growth rate as the eigenvalue and perturbation
ˆ ˆ ˆ ˆ
amplitudes u, v, w, p as the eigenfunctions.





(13)





(14)

ˆ
ˆ ˆ
u   V *  u   v * U  M a 1

ˆ
p
1 
 2
ˆ
ˆ
ˆ 
 Re1  *  k 2 u 
*  v  kw
x
3 x



ˆ
ˆ ˆ
v   V *  v   v * V  M a 1

ˆ
 2

p
1 
ˆ
ˆ
ˆ
 Re1  *  k 2 v 
*  v  kw
y
3 y












2
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
 w   V *  w  M a 1kp  Re1  *  k 2 w  k  *  v  kw 

1
3

ˆ
ˆ ˆ
ˆ
ˆ
 p   V *  p   v *  P  M a 1  *  v  kw  0

8

(15)
(16)

For the stability analysis of the separated flow, the following boundary conditions are
implemented: all disturbance velocity components are zero on the solid wall and in the
freestream. On the solid wall, the compatibility conditions (Equations (17) and (18) below) are
derived from the equations for the pressure perturbation. In Equation (17), I* is the twodimensional identity matrix. All disturbances at the inflow boundary are zero, which, from a
physical point of view, serves to ensure that no perturbations other than those due to potentially
self-excited global eigenmodes enter the separated flow region. Linear extrapolation from the
interior of the domain is applied for obtaining the disturbances on the outflow boundary.

 1
1 
1

ˆ
ˆ
ˆ
ˆ
pI*   v  I*  *  v  kw    n  0

Re 
3

 Ma

ˆ
w  n  0

(17)
(18)

Finite element method is used to discretize Equations (13)-(16). With this method, the
imposition of the boundary conditions results in a real non-symmetric generalized eigenvalue
problem in the form of

AΦ  BΦ ,

(19)

where A and B are the assembly mass matrices calculated from the base flow and

ˆ ˆ ˆ ˆ
  u, v, w, p is the assembly vector of the eigenfunctions. Details on how A and B are
computed are presented in APPENDIX A. According to linear stability theory, the stability
properties of the flow depend on the eigenvalue  . When  is real, the disturbances either grow
or decay monotonically and the critical Reynolds number is one for which   0 . When  is
complex, the neutral condition is r  0 , and the onset of instability is oscillatory with
dimensionless angular frequency i . This normal mode form also includes the time-dependent
two-dimensional instability of the steady flow for which k  0 .

9

With eigenmodes obtained from stability computation, the total three-dimensional flow
field can be reconstructed on the basis of the two-dimensional base flow with an arbitrary
number of disturbance eigenmodes. When these eigenmodes are chosen to be conjugate pairs of
most unstable modes of an eigenspectrum, the three-dimensional flows can be expressed by the
following superposition of eigenfunctions (the pressure and y-velocity v have a similar form as
that of u )
N

 

n

 

ˆ
ˆn
u  U  2  etr cos  kz  uin cos tin  ur sin tin 




n1

N

 

n

 

ˆn
ˆ
w  2  etr cos  kz   wr cos tin  win sin tin 




n1

(20)

(21)

where subscripts r , i , n , N denotes the real and imaginary parts of the eigenfunctions, the nth
pair of eigenmodes and the total number of conjugate pairs of eigenmodes, respectively.
2.2

Local Stability Formulation
As the first validation step of global stability analysis, local stability analysis is carried out

for the same base flow, and its results are compared to those of global stability analysis. The base
flow closely satisfies ∂/∂x<<∂/∂y at all locations, therefore it is reasonable to assume that in a
short streamwise distance, the flow over airfoil is a parallel flow, which justifies the application
of local stability analysis to local x-velocity profiles U at different streamwise locations.
For a parallel viscous flow U  y  , with perturbations assumed to be in form of Equation (22)
below, the temporal local stability variables are governed by the Orr-Sommerfeld Equation (23).
For simplification, the viscosity is neglected and the Orr-Sommerfeld equation is reduced to the
Rayleigh equation (24) [39] with boundary conditions described in Equation (25).

10

ˆ
ˆ
ˆ
 u ', v ', w ', p '  u  y  , v  y  , w  y  , p  y  exp ikx  t 
ˆ











i

2
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
v '''' 2k 2v '' k 4v  0
 U  i  v '' k v  U '' v 
k
k


(22)
(23)


U '' 
ˆ
ˆ
v ''  k 2 
v  0
U  i k 


(24)

ˆ
ˆ
v  0   0 and v  y     ek y

(25)

A smoothing-spline method is developed in this study to enable exact fitting of any base
flow profile obtained from experiments or computations with smooth first- and second-order
derivatives. Simple analytical expressions of the velocity profiles were used by others [40-42] to
perform local stability analysis. However, it is not practical to find simple analytical profiles for
all experimental or simulation data. Inaccurate representation of the base flow could cause
significant error in the stability analysis. The smoothing-spline method gives a nearly perfect
fitting and is used in the current study. This fitting method is implemented in the eigensolver
developed in this study. A shooting method [43] is used in the eigensolver to solve the Rayleigh
equation.
2.3

Base Flow Results
The base flow used in the present work is the same as that obtained by Galbraith et al. [7]

with high-order implicit large eddy simulation (LES) method. The airfoil flow has a Reynolds
number of 60,000 and an AoA of 4° The computational domain (shown in Figure 2 as an
.
enclosed rectangle) for the global stability analysis covers the domain in the range of X/C ∈
(0.1941, 1.0156) and Y/C ∈ (0.0011, 0.2914), which ensures that the whole LSB is enclosed for

11

the global stability analysis. For local stability analysis, x-velocity profiles at different
streamwise locations in the same base flow are used.
0.4

0.2
Y
0
-0.2
-0.2

0.4
0.6
0.8
1
X
Figure 2: Base flow around a SD7003 airfoil, Re=60,000, α=4° [7]. Contours represent timeaveraged x-velocity. The rectangular solid line shows the domain used in the global stability
analysis (For interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this dissertation).

2.4

0

0.2

Local Stability Results
To obtain the most unstable modes, local stability analysis is applied to a series of x-

velocity profiles at different streamwise locations with an increment of 0.05C. The most unstable
mode is found to occur with the velocity profile at x/C=0.5 shown in Figure 3. The Rayleigh
Equation (24) is solved with boundary condition described in Equation (25) for the eigenmodes.
Details of the numerical solution procedure can be found in reference [42]. The most unstable
mode of is calculated to be ω=17.94+41.48i (Figure 4), corresponding to a most unstable
nondimensional frequency (normalized by C/U0) of 6.6.

12

0.25

y/C

0.2
0.15
0.1
0.05
0

0

0.5
U/U
Figure 3: Time-averaged x-velocity (U) profile at x/C=0.5.

80

r=17.94

18

1

60

i

r

16
14

40

12
40

i= 41.48

20
60
80
100
Streamwise wave number k

120

40

60
80
100
Streamwise wave number k

120

Figure 4: Local stability analysis results for time-averaged x-velocity (U) profile at x/C=0.5; (left)
growth rate vs. wavenumber, (right) angular frequency vs. wavenumber.

2.5

Global Stability Results
A computational model/code is developed based on the formulation in Section 2.1. To

validate our code, a series of computations on the global stability of the lid-driven cavity flow
are performed with the same equations, different base flows and boundary conditions. The
results are found to be in agreement with those obtained by Ding and Kawahara [37] and
Theofilis [44].

13

The effect of grid resolution on the solution is tested by considering six different meshes
with the following resolutions: 29× 41× 47× 56× 70× and 76× (Figure 5). The
21,
30,
35,
41,
51,
80
largest difference between the results obtained by meshes 70× and 76× is about 3%, which
51
80
indicates that the effect of mesh resolution is relatively small for the 70× mesh. This non51
uniform mesh is used for global stability analysis.
70
60
50

i

40

mesh 29 21
mesh 41 30
mesh 47 35
mesh 56 41
mesh 70 51
mesh 76 80

30
20
10
0
0

5

10

r

15

20

Figure 5: The four most unstable modes at wavenumber k=20, obtained with different base flow
grids.

For each wavenumber k, the most unstable mode is sought. A wide range of wavenumber
values (0.2, 90) is considered. It should be noted that the upper limit of k is relatively high due to
the fact that the airfoil chord is used as the length scale in our analysis ( k  2 C Lz ). If a smaller
length scale is used, the wavenumbers will be proportionally reduced. This wavenumber range
corresponds to a range of wavelength normalized by the airfoil chord Lz C (0.070, 31.4). The
relation between the most unstable eigenvalue, the one with the largest real part, and the
wavenumber is shown in Figure 6. It is shown that the growth rate decreases as wavenumber

14

increases. In Equations (13) to (15), the terms -k and  are on opposite sides. Therefore, for k
2

larger than 90, the growth rate of the most unstable mode is expected to decrease further. Since
the real part of the growth rate of the most unstable mode is positive for all k values, the flow is
unstable. Moreover, the angular frequency increases as the wavenumber increases, indicating
that the shorter wavelength perturbations are less unstable. With wavenumber k being close to
zero, the most unstable mode reaches a plateau of ω=18+59.4i at which the angular frequency
59.4 corresponds to a nondimensional frequency of 9.5. The maximum growth rates from the
global and local analyses are roughly in agreement, and both methods predict an unstable flow.
However, the most unstable frequency from the global analysis is higher than that predicted by
the local analysis. This difference is caused by the combined effect of using local velocity
profiles and neglecting viscosity in the local stability analysis. Galbraith et al. [7] monitored in
their simulation the time variation of velocity at a point directly above the airfoil surface at
x/C=0.5, and plotted the turbulent kinetic energy spectrum at that point (Figure 7). They pointed
out that the nondimensional frequency of 9 may be the most unstable frequency of the flow even
though it is not the dominant frequency observed in the velocity energy spectrum. Their later
stability analysis [20] also yields a most unstable frequency close to 9.

15

r

20

15

10
0
64

20

40

60

80

100

40
60
80
Spanwise wave number k

100

i

62
60
58
0

20

Figure 6: The most unstable modes vs. wavenumber; (upper) maximum growth rate, (lower)
angular frequency.

0.3
0.25
0.2
E

0.15
0.1
0.05
0
5

10
15
20
ω
Figure 7: Turbulent-kinetic-energy frequency spectra at x/C=0.5 (about center of LSB) [7].

16

According to Equations (20) and (21), the three-dimensional velocity field can be
reconstructed by adding the summation of a number of eigenmodes to the base flow. Since linear
instability is assumed, reconstructed velocities do not necessarily resemble the real flow. But the
emphasis here is on the distribution of velocity disturbances and their relation. Figure 8 shows
the reconstructed x-velocity given by Equation (20) for k=1, obtained with the conjugate pair of
the most unstable eigenmodes. The eigenvalue for the most unstable mode is ω=18+59.4i with a
nondimensional period of T=0.1058. Figures 8 shows the velocity at times 0, T/4, T/2, and 3T/4.
Since the growth rate of the most unstable mode is greater than zero, the stability analysis
predicts that the disturbances will amplify in time and become dominant over the base flow. The
increasing range of the color bar in Figure 8 confirms this. Figure 9 shows the reconstructed
disturbances (without the base flow velocity) of the three velocity components at T/2 obtained
with the same eigenmodes. Evidently, the three velocity disturbances share some similar features.
First, the disturbances are zero at inlet, at airfoil surface and at far-field boundaries, as enforced
by the boundary conditions. It should be noted that disturbances are also zero at the outlet, which
is not restricted by the boundary condition. Second, all disturbances are confined to the region
0.2<x/C<0.65 and y/C<0.2, which corresponds to the LSB region where reverse flow occurs in
the base flow. The maximum magnitude and maximum gradient of disturbances occur around
x/C=0.4. It is found that the z-velocity disturbance is much smaller in magnitude than the
disturbances of x- and y-velocities. Similar differences in disturbances were also found in the
stability analysis of lid-driven cavity flow by Ding and Kawahara [37]. For wavenumbers other
than 1, the qualitative structure of the velocity disturbances, calculated based on the pair of the
most unstable eigenmodes, is similar to that shown in Figure 9. Therefore, the above remarks are
valid at different wavenumber values.

17

The spanwise eigenmodes indicate the existence of oscillatory flow structure perpendicular
to the base flow plane. This can be verified by the instantaneous z-velocity contours shown in
Figure 10 at a z-section plane which are obtained from the 3D LES data of a fully developed
flow.

u
t=0

1

0.1
Y/C

0.12

1.5

t = T/4

2
1

0.08

0

0.06

0.5

0.08

0.1
Y/C

0.12

u

-1

0
0.06
-0.5
0.04
0.2

0.4
X/C

0.04
0.2

0.6

0.4
X/C

0.6

u
2
1

Y/C

0.1

0.12

-1

0
0.08
-2

-2

0.06

0.06
-4

-3
0.04
0.2

0.4
X/C

4

0.1

0
0.08

t = 3T/4

2
Y/C

0.12

t = T/2

u

0.04
0.2

0.6

0.4
X/C

0.6

Figure 8: Reconstructed total x-velocity using the conjugate pair of most unstable modes for k=1;
(a) t=0, (b) t=T/4, (c) t=T/2, (d) t=3T/4.

18

u

0.2
t = T/2

0
0.1

2

0.15

Y/C

Y/C

t = T/2

1

0.15

v

0.2

2

-1

1
0

0.1

-1

-2
0.05
0.2

0.05

-3
0.4

0.6

-2

0.8

0.2

X/C

0.4

0.6

0.8

X/C
w

0.2

0.05

t = T/2

Y/C

0.15
0
0.1
-0.05
0.05
-0.1
0.2

0.4

0.6

0.8

X/C

Figure 9: Reconstructed velocity disturbances using the conjugate pair of most unstable modes
for k=1 at t=T/2.

w Um

Figure 10: Instantaneous normalized z-velocity w U m on a z-section, α=4°, Re=60,000.

19

3.
3.1

Simulation Method

Background
In addition to stability analysis discussed in the previous chapter, detailed simulations of

steady and unsteady flows over SD7003 airfoil are conducted with the large eddy simulation
(LES) method. Reynolds-averaged Navier-Stokes (RANS) models are much computationally
less demanding than LES but they can only predict the ensemble- or time-averaged behavior of
the flow. Furthermore, RANS solutions greatly depend on empirical model constants and often
fail to capture some of the key features of the flow such as flow separation and reattachment. In
LES, large eddies of the flow are explicitly computed by solving the “filtered” Navier-Stocks
equations, while smaller eddies are modeled using a subgrid-scale (SGS) model. The rationale
behind LES is that by modeling less, the error caused by the modeling of turbulence can be
reduced. It is also believed to be easier and more accurate to model small scales, since they are
more isotropic and less affected by the boundary conditions. Although LES usually requires
substantially more computational resources than RANS, it is computationally less demanding
than direct numerical simulation (DNS) method, and hence can be applied to relatively high
Reynolds number flows.
The governing equations for LES are obtained by filtering the time-dependent NavierStokes equations. The filtering process effectively removes eddies with sizes smaller than the
filter width or grid spacing used in the computation. In some of LES calculations such as those
conducted by the FLUENT code [45], the filtering is implicit, which means that the filtered
equations do not explicitly show the information on the filter shape or width. The information is
implicitly contained in the SGS model and the discretization of the equations.

20

3.2

Filtered Equations
In any LES method, the resolved or filtered value of variable  is computed as

  x      x ' G  x, x 'dx ' ,
D

(26)

where D is the flow domain, and G is the filter function which satisfies the normalization
function  G  x, x 'dx '  1 . There are a variety of filter functions such as Gaussian filter, Fourier
D

cut-off filter, or box filter [46] available for LES. In the adopted code, the finite-volume
discretization itself implicitly provides the box filtering operation as

  x 

1
V

V   x 'dx ' ,

(27)

where V is the volume of a computational cell. The implied box filter function is over the
physical space, and it has the following form
x 'v
1 V ,
G  x, x '   
.
 0, otherwise

(28)

The filtered continuity and momentum equations for an incompressible and Newtonian fluid
take the following form [47]:
ui
0
xi





ui  ui u j
1 p




t
x j
 xi x j

(29)

 ui   ij


 x j  x j



(30)

where

 ij  uiu j  ui u j
is the unclosed subgrid-scale (SGS) stress tensor.

21

(31)

Often, the SGS turbulence models used for  ij are based on Boussinesq hypothesis [48]
and are expressed by the following equation
1
3

 ij  ij kk  2 t S ij ,

(32)

where  t is the subgrid-scale eddy viscosity and S ij is the resolved-scale strain rate tensor,

1  u i u j
S ij  

2  x j xi







(33)

With a model for  t , Equation (30) is closed and may be solved for the resolved variables.
3.3

Subgrid-Scale Model
The subgrid-scale model used in this study is the Wall-Adapting Local Eddy-viscosity

(WALE) model of Nicoud et al. [49]. Compared to dynamic Smagorinsky and Kinetic Energy
Transport SGS models, WALE model is relatively simple yet sufficiently accurate for the flows
considered in this study. The model is based on the square of the velocity gradient tensor. The
eddy viscosity is obtained from the following equation

 Sijd Sijd 
2
 t   Cw  
,
52
d d 54
 S ij S ij    Sij Sij 
32

(34)

d
where Sij is defined as
d
Sij 





2
2
1 2
1
g ij  g ji   ij g kk ,
2
3

g ij 

22

ui
x j

(35)

(36)

In Equation (34), Δ is the filter width computed as the cubic root of the computational cell
volume, and Cw is the WALE model constant. A Cw value of 0.325 has been found to yield
satisfactory results for a wide range of flows [45].
The WALE SGS model is a mixed model based on local rate of strain and rotational
tensors, thus, suitable for prediction of the SGS kinetic energy dissipation. The eddy-viscosity
goes naturally to zero in the model in the vicinity of the wall, and the model is able to return the
correct wall asymptotic ( y 3 ) behavior for wall bounded flows; therefore no dynamic procedure
is required. Additionally, the WALE model generates zero eddy viscosity in laminar flows, while
reproducing laminar to turbulent transition through the growth of linear instability modes.
Guleren [50] tested several SGS models and concluded that WALE model performs generally
better than other tested SGS models.
A comparison between WALE and other SGS models is made in Figure 11, where the mean
surface pressure coefficient over SD7003 airfoil as computed with different SGS models are
compared with the LES results generated by Galbraith et al. [7], whose results are in good
agreement with the experiment. Figure 11 shows that the WALE model yields a very good
agreement. Based on these findings, all simulations considered below are conducted with the
WALE SGS model.

23

1

-Cp

0.5

0
Galbraith
Smagorinsky
Dynamic Smagorinsky
KET
Dynamic KET
WALE

-0.5

-1
0

0.2

0.4

0.6

0.8

1

X/C

Figure 11: Comparison of time- and spanwise-averaged pressure coefficient over the SD7003
airfoil, obtained by our LES with different SGS models and compared with those obtained by
Galbraith et al. [7]. The flow Reynolds number based on freestream velocity and airfoil chord is
60,000, and AoA=4° KET stands for Kinetic Energy Transportation model.
.

3.4

Grid and Boundary Conditions
Initially, a structured 2D O-grid (Figure 12a) of size 315×
151 is generated for the SD7003

airfoil with a rounded trailing edge (Figure 12c). To generate accurate solution, the mesh is made
to be highly orthogonal and is stretched rapidly from the surface of the airfoil outward (about
70% of the grid points are within a chord distance away from the airfoil center) to ensure high
resolution of the flow near the airfoil (Figure 12c). The far field boundary is made to be almost a
circle with the center at the chord mid-point and a radius of about 27 chords. The 2D mesh is
copied in the spanwise (z) direction to form a 3D grid. The 3D grid has 101 nodes in the
spanwise direction, and thus the total grid size is 315×
151×
101. The grids are evenly spaced in

24

the spanwise direction with a spanwise extent of 0.2C. The grid is similar to that used by
Galbraith [20].

20

0.4

Freestream
10

0.2

Y 0

Y
Out flow

-10

0
-0.2
-0.4

-20
-20

-10

0

X

10

20

30 (a)

0

0.5

1

X

(b)

0.004
0.002
Y 0
-0.002

0.986

0.99

X

0.994

(c)

0.2C

(d)

Figure 12: Computational grid: (a) grid on XY plane, (b) grid in the vicinity of airfoil, (c) grid
near the trailing edge, (d) 3D airfoil.

Solutions for the steady flow over SD7003 airfoil were obtained using the 3D 315×
151×
101
grid and three other coarser and finer grids [20]. It was found that the 315×
151×
101 grid
captures the flow features well and generates results that closely match those obtained based on
substantially finer grids, even though the substantially finer grids provide slightly better details

25

of the flow in the transition region. The 315×
151×
101 grid provides an adequate flow resolution
for the quantities reported here and is used in this study.
Boundary conditions are specified in the following manner. On most of the far-field
boundary, freestream velocity conditions with no perturbations are prescribed with steady or
unsteady flow conditions for velocity and pressure. For wake region of the outer boundary, we
use the outflow boundary condition in which variables on the boundary are extrapolated from the
interior of the computational domain. The airfoil surface is modeled with a no-slip adiabatic wall
boundary condition. A spanwise periodic boundary condition is imposed on the side faces,
which, equivalently, requires the two side faces to be frictionless and have the same flow
properties (e.g. velocity and pressure).
In all simulations considered in this study, sufficient time is given to the flow to become
fully developed before any data are recorded or analyzed. A flow with steady freestream is
considered fully developed when lift and drag forces become stabilized and their moving
averages over two residence times show no change in time. For a flow with periodic freestream
to be considered fully developed, the monitored lift and drag forces are required to become
periodic and their peak values become steady. In order to reduce the computational time
associated with initial development of the 3D flow, a fully developed 2D flow is generated first,
and the 2D data is then transferred to the 3D mesh and used as the initial condition. The non-4

dimensional time step ΔtU∞/C is kept constant at 1.37×
10

which is sufficiently small for

resolving the temporal variations in the flow.
The momentum equations are discretized using a second-order bounded central difference
scheme, and the continuity equation is discretized using a second-order central scheme. With
very fine mesh near the airfoil surface, implicit time integration is chosen over explicit methods

26

in order to avoid very small time steps. The transient terms are discretized using second-order
implicit scheme.
3.5

Separation, Transition, Reattachment, and Force Coefficients
The mean separation and reattachment points are where the mean skin friction coefficient

changes signs. To define the transition point we follow the approach described in reference [8].
Laminar to turbulent transition is assumed to occur where the normalized Reynolds shear stress
reaches 0.1% of the free steam velocity and grows significantly after that [8].
For flows with steady freestream, the x-direction is defined as the steady freestream
direction, and the y-direction is perpendicular to the x-direction, as shown in Figure 13. The lift
force is defined as the total force on the airfoil in the y-direction, and the drag force is defined as
the total force on the airfoil in the x-direction. The lift and drag coefficients are defined as
2
2
CL  L ( U C 2) and CD  D ( U C 2) where L is the total lift force per unit span length,

D is the total drag force per unit span length, and U  is the freestream velocity. The pressure
2
coefficient CP is calculated as CP  ( p  p ) ( U  2) , where p is the freestream pressure.

The skin friction coefficient CF is calculated as CF  

u
2
( U  2) .
n

For flows with harmonically oscillating freestream velocity magnitude or direction, the xdirection is defined as the mean freestream direction, and the y-direction is perpendicular to the
x-direction. The lift and drag forces are also obtained by decomposing the total aerodynamic
force along the x-direction. The same definition of lift and drag for flows with oscillating
freestream velocity direction was used by other researchers [34, 36]. For flows with a steady
freestream and an oscillating airfoil (not considered in this study), it is not uncommon to
decompose the drag and lift along the freestream direction even though the effective AoA is
27

constantly changing due to the airfoil oscillation [23, 54]. Thus, CL, CD, CP, and CF are defined
the same as for flows with steady freestream, except that U  is now changed to U m , the mean
freestream velocity magnitude. Unlike flows with steady freestream where pressure gradient is
zero, flows with unsteady freestream have time-dependent pressure gradient related to the
freestream velocity acceleration in either velocity magnitude or direction. Therefore for flows
with unsteady freestream, the reference pressure p in the CP definition is chosen to be the
pressure at the point (0, 25C, 0.1C), on the y-axis and 25 chords away from the airfoil leading
edge.
All CL, CD, CP, and CF calculated in our study are based on either both time- and spanwiseaveraged data or just spanwise-averaged data. Figure 13 shows a schematic view of the flow
directions for both steady and unsteady flows.

y

x


(t)
t U or U

f
f

Constan

Figure 13: Schematic of flow over airfoil, showing the freestream and the coordinate system.

28

3.6

Comparison of 2D and 3D Simulations
The effect of flow three-dimensionality flow on the aerodynamic forces, the separation and

reattachment of the boundary layer, and the flow structure is discussed in this section before
considering the 3D steady and unsteady LES results in detail. For the flow with steady
freestream, Re=60,000 and AoA=4° 2D and 3D solutions are compared in Table 1 and Figures
,
14-16. Table 1 shows that the aerodynamic forces obtained from 2D solution are similar to that
of 3D solution, though 2D simulation predicts a slightly longer LSB. Surface pressure coefficient
obtained from the 2D solution also agrees reasonably well with the 3D solution. However, the
separation occurs earlier in 2D simulation and the reattachment point is predicted by the 2D
model to further downstream than that predicted by the 3D model. Also, the skin friction
coefficient from the 2D solution does not increase to the same level as that of the 3D solution
downstream of the reattachment point. The maximum difference in friction coefficient between
2D and 3D solutions is 5.8% of the total variation of CF and it occurs at X/C=0.67. Figure 16
shows that the difference observed in the spanwise-averaged instantaneous 2D and 3D z-vorticity
contours is significant. In the 2D solution, the separated boundary layer develops into a series of
large-scale connected vortices which do not breakdown to smaller structures, consistent with the
observations made by Galbraith [20] and Hodge et al. [51]. The 2D and 3D results for larger
angle of attack of AoA=8° (but same Reynolds number) in Table 2 and Figures 17-18 show
trends similar to those for flows with AoA=4° even though the differences are more noticeable
,
at higher angle of attack.

29

Separation Transition Reattachment Lift coeff. Drag Coeff.
Xs C
Xt C
Xr C
CL
CD
LES, 2D 0.23
0.44
0.67
0.55
0.021
LES, 3D 0.25
0.50
0.64
0.59
0.020
Date Set

Table 1: Mean quantities obtained from 2D and 3D solutions for the flow with steady freestream
and Re=60,000 and AoA=4°
.

2D
3D

1

-CP

0.5

0

-0.5

-1
0

0.2

0.4

0.6

0.8

1

X/C

Figure 14: Time-averaged surface pressure coefficient obtained by 2D and 3D simulations of the
flow with steady freestream and Re=60,000 and AoA=4°
.

30

0.03
2D
3D

0.025
0.02

CF

0.015
0.01
0.005
0
-0.005
-0.01
0

0.2

0.4

0.6

0.8

1

X/C

Figure 15: Time-averaged skin friction coefficient on the upper surface of the airfoil obtained by
2D and 3D simulations of the flow with steady freestream and Re=60,000 and AoA=4°
.

Figure 16: Instantaneous normalized z-vorticity ( z C U m ) contours obtained by 2D and 3D
simulations of the flow with steady freestream and Re=60,000 and AoA=4°
.

31

Separation Transition Reattachment
Xs C
Xt C
Xr C
LES, 2D 0.037
0.15
0.29
LES, 3D 0.042
0.16
0.25
Date Set

Lift coeff.
CL
0.91
0.93

Drag Coeff.
CD
0.042
0.040

Table 2: Mean quantities obtained from 2D and 3D solutions for the flow with steady freestream
and Re=60,000 and AoA=8°
.

3
2D
3D

2.5
2

-CP

1.5
1
0.5
0
-0.5
-1
0

0.2

0.4

0.6

0.8

1

X/C

Figure 17: Time-averaged surface pressure coefficient obtained by 2D and 3D simulations of the
flow with steady freestream and Re=60,000 and AoA=8°
.

32

0.03
2D
3D
0.02

CF

0.01

0

-0.01

-0.02
0

0.2

0.4

0.6

0.8

1

X/C

Figure 18: Time-averaged skin friction coefficient on the upper surface of the airfoil obtained by
2D and 3D simulations of the flow with steady freestream and Re=60,000 and AoA=8°
.

The differences between 2D and 3D solutions are further examined for unsteady freestream
flows over the airfoil. Comparison is made among three cases with oscillating freestream AoA.
In each case, the velocity direction or AoA varies sinusoidally in the x-y plane as

  t   4  8  sin  2U mt C  where U m is the constant freestream velocity magnitude. The zvelocity component is zero at freestream. The Reynolds number is calculated as Re  U mC  ,
and it is 10,000, 30,000 and 60,000 for the three cases simulated, respectively. The 2D and 3D
values of CL, CD and z-vorticity iso-surfaces obtained by LES are compared in Figures 19-24.
As explained in Section 3.5, CL and CD are calculated by normalizing the aerodynamic forces
with U m . In the z-vorticity iso-surface figures 20, 22, and 24, t0 is the beginning of a cycle when
the AoA is at its mean value and increasing. When showing the vorticity iso-surfaces for the 2D
solutions, the 2D domain is extended in the spanwise direction to the same spanwise size as the

33

3D domain for easier comparison. For all three Reynolds numbers, 2D and 3D solutions provide
similar CL and CD profiles except in the peaks and valleys. In the 2D profiles, the peaks and
valleys are somewhat “chopped off”, and the amount “chopped off” increases slightly at higher
Reynolds numbers. This is directly related to the difference in static pressure between 2D and 3D
solutions. For the airfoil flows studied here, static pressure contributes vast majority of the lift
(or drag).
For the unsteady case with oscillating AoA and Re=60,000, the pressure coefficient CP on
the airfoil surface at different times of a cycle is plotted in Figure 25. It can be seen that at times
t0 and t0+T/2 when there is a significant acceleration in freestream AoA, the CP curves of the 2D
solution are narrower than those of the 3D solution, which causes the “chop-off” in the 2D CL
and CD curves. There is no significant difference in the overall shape of CP curves obtained from
2D and 3D solutions at times t0+T/4 and t0+3T/4. This indicates that lack of spanwise instability
in the 2D simulations prevents the pressure in the separation region from reaching the same low
and high values as in the 3D simulations at times of fast changing freestream AoA. This
difference between 2D and 3D solutions also happens at relatively low Reynolds number of
10,000, which indicates that the spanwise instability has a significant effect on pressure even in
laminar flows.
The z-vorticity iso-surfaces obtained from 3D simulations clearly show that the flow
changes from laminar to turbulent as Reynolds number increases, while 2D solutions fail to
capture this trend. At lower Re=10,000, the 2D and 3D solutions show very similar flow
structures that are laminar and similar in spanwise direction. At higher Re=30,000, the 3D
simulation shows spanwise fluctuation, and breakdown of vortices near the trailing edge while a
large portion of the flow over the airfoil remains laminar. As Reynolds number further increases

34

to 60,000, the 3D simulation shows that the large-scale vortices break down earlier. At both
Reynolds numbers of 30,000 and 60,000, the 2D simulations predict the flows to be laminar. As
Reynolds number increases, the flow structures predicted by 2D and 3D simulations become less
similar; at Re=30,000, the similarity is still noticeable, but it is much less so at Re=60,000.

35

4

0.25
0.2
0.15

2D
3D

3

CD

CL

2
1
0

0.1
0.05
0

-0.05
-0.1

-1
-2
0

2D
3D

1

2

-0.15
0

3

1

t/T

2

3

t/T

Figure 19: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=10,000.

u/Um
3D

2D

2.00
1.75
1.38
1.00
0.62
0.25
-0.12
-0.50

Figure 20: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations of
the flow with oscillating AoA and Re=10,000 (colored by u/Um for better visibility).

36

4

0.25
0.2
0.15

2D
3D

3

CD

CL

2
1
0

0.1
0.05
0

-0.05
-0.1

-1
-2
0

2D
3D

1

2

-0.15
0

3

t/T

1

2

3

t/T

Figure 21: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=30,000.

3D

2D

u/Um
2.00
1.75
1.38
1.00
0.62
0.25
-0.12
-0.50

Figure 22: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations of
the flow with oscillating AoA and Re=30,000.

37

4

0.25
0.2
0.15

2D
3D

3

CD

CL

2
1
0

0.1
0.05
0

-0.05
-0.1

-1
-2
0

2D
3D

1

2

-0.15
0

3

t/T

1

2

3

t/T

Figure 23: Lift and drag coefficients obtained by 2D and 3D simulations of the flow with
oscillating AoA and Re=60,000.

3D

2D

u/Um
2.00
1.75
1.38
1.00
0.62
0.25
-0.12
-0.50

Figure 24: Iso-surfaces of z-vorticity ( z C U m  10 ) obtained by 2D and 3D simulations of
the flow with oscillating AoA and Re=60,000.

38

The spanwise domain size Z/C=0.2 is selected based on Galbraith and Visbal findings [20,
23]. They compared solutions with spanwise extent ranging from 0.1C to 0.3C and concluded
that Z/C≥0.1 is sufficient for capturing flow details in the spanwise direction. To further verify
this, the size of flow vortices in planes perpendicular to the freestream direction was examined in
the most “extreme” case with Re=60,000 and AoA=16°which involves relatively large vortical
flow structures in comparison to cases with smaller AoAs. Contours of x-vorticity obtained from
the LES data for AoA=16°case at several x/C locations are plotted in Figure 26, which shows
that the size of almost all large-scale vortices in the spanwise direction is less than 1/5 of the
spanwise extent. Therefore, flow properties in the spanwise direction are not affected by the
chosen spanwise size.
Based on the above discussion, the 3D mesh with spanwise extent Z/C=0.2 will be used for
all computations in the rest of this thesis.

39

80
2D
3D

78
76
0

0.2

0.4

t0
0.6

0.8

1

5

t0+T/4

-5

C

P

0
0

0.2

0.4

0.6

0.8

1

-50
-52

t0+T/2

-54
0

0.2

0.4

0.6

0.8

1

1

t0+3T/4

0
-1

0

0.2

0.4

0.6

0.8

1

X/C

Figure 25: Pressure coefficient obtained by 2D and 3D simulations of the flow with oscillating
AoA and Re=60,000.

Spanwise Direction

x C  0.2

x C  0.5

x C  0.8

height C  0.2

width C  0.2

width C  0.2

width C  0.2

Figure 26: Instantaneous contours of x-vorticity ( x C U m ) obtained by 3D simulations of the
flow with steady freestream and Re=60,000, AoA=16°
.

40

4.

Large Eddy Simulations of Steady Freestream Flows over SD7003 Airfoil

To validate the computational method and to establish reference cases for comparison with
unsteady flows considered in the next chapter, flows around SD7003 airfoil with steady
freestream are computed in this chapter and compared to available experimental and
computational data. The z-vorticity iso-surfaces of flows with Re=60,000 and AoA=4° 8° 12°
, ,
and 16°as obtained by LES in Figure 27 show the general expected trend that at higher AoA the
size of recirculation region over the airfoil is bigger and so are the sizes of the vortices formed by
the separation. For steady flow conditions, high resolution velocity and Reynolds stress
measurements have been reported by Radespiel et al. [8] whose experiments were conducted in a
water channel and a low-noise wind tunnel at Technical University of Braunschweig (TU-BS).
Detailed particle image velocimetry (PIV) measurements for the SD7003 airfoil were also
conducted by Ol et al. [9] at the Air Force Research Laboratory (AFRL) water channel. Katz et
al. [10-13] at Turbulent Mixing and Unsteady Aerodynamics Laboratory (TMUAL) at Michigan
State University also measured the flow over SD7003 airfoil in a water channel with the
molecular tagging velocimetry (MTV) method. Galbraith et al. [7] simulated steady flows over
SD7003 airfoil with LES and a high-order low-pass filter as a SGS model. Results obtained by
two compact differencing schemes (2nd-order and 6th-order) and two meshes (base-line mesh
and overset mesh) were reported. The base-line mesh used in reference [7] is similar to the one
used in our study, while the overset mesh has a refined region around the upper surface of the
airfoil.
For the steady flow with Re=60,000 and AoA=4° the predicted separation, transition, and
,
reattachment points and lift and drag coefficients are compared to the numerical and
experimental results in Table 3. These variables are defined in Section 3.5.

41

AoA=4°

AoA=8°
u/Um
2.00
1.75
AoA=12°
1.38
1.00
0.62

AoA=16°

0.25
-0.12
-0.50
Figure 27: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for the flow with steady
freestream, Re=60,000, AoA=4° 8° 12°and 16°
, ,
.

Separation Transition Reattachment Lift coeff. Drag Coeff.
Xs C
Xt C
Xr C
CL
CD
LES
0.25
0.50
0.64
0.59
0.020
Galbraith 2nd-order [20]
0.25
0.46
0.66
Galbraith 6th-order [20]
0.23
0.55
0.65
0.59
0.021
TU-BS [8]
0.30
0.53
0.62
AFRL [9]
0.18
0.47
0.58
0.56
Date Set

Table 3: Computed and measured time-averaged data for the flow with steady freestream,
Re=60,000 and AoA=4° LES refers to our simulations with WALE SGS model.
.

42

There are some quantitative differences between the computed data and the AFRL
experimental data. According to Galbraith and Ol [20], some of these differences may be
attributed to the fact that the tested airfoil had a tendency to vibrate at certain AoAs, promoting
premature flow transition.
The LES predictions are generally in agreement with those measured in TU-BS wind tunnel,
but seem to predict earlier separation and latter reattachment, thus a larger LSB. Similar trends
are reported by Radespiel et al. [14], Yuan et al. [16], and Galbraith [20]. Part of the differences
may be caused by the freestream turbulence intensity in the experimental facility [12, 29, 52].
Another possible reason is the reported discrepancy between the nominal and effective angles of
attack due to interference effects in the wind tunnel.
Figures 28 and 29 compare the surface pressure coefficient CP and skin friction coefficient
CF obtained from LES to those by Galbraith [20]. Our LES results obtained with a 2nd-order
scheme are expected to be comparable to those of Galbraith obtained with the 2nd-order scheme.
The comparison shows overall good agreement between the two LES results. However, our LES
predictions seem to be closer to Galbraith baseline 6th-order and overset 6th-order predictions,
instead of their baseline 2nd-order results.

43

1

-CP

0.5

0
LES
Galbraith baseline 2nd-order
Galbraith baseline 6th-order
Galbraith overset 6th-order

-0.5

-1
0

0.2

0.4

0.6
0.8
1
X/C
Figure 28: Time-averaged surface pressure coefficient for the flow with steady freestream,
Re=60,000 and AoA=4°
.

0.03
0.025

LES
Galbraith baseline 2nd-order
Galbraith baseline 6th-order
Galbraith overset 6th-order

0.02

CF

0.015
0.01
0.005
0
-0.005
-0.01
0

0.2

0.4

0.6

0.8

1

X/C
Figure 29: Time-averaged skin friction coefficient on the upper airfoil surface for the flow with
steady freestream, Re=60,000 and AoA=4°
.

44

Table 4 and Figures 30 and 31 compare LES results for higher AoA of 8°and Reynolds
number of 60,000. The trends and conclusions are similar to those shown above for the case with
AoA=4°and Re=60,000. There seems to be a slight shift in the CP and CF curves predicted by
the two LES models which is due to differences in the SGS model, mesh, and discretization
method.

Separation Transition Reattachment Lift Coeff. Drag Coeff.
Xs C
Xt C
Xr C
CL
CD
LES
0.042
0.16
0.25
0.93
0.040
Galbraith 6th-order
0.04
0.18
0.28
0.92
0.043
Date Set

Table 4: Computed time-averaged data for the flow with steady freestream, Re=60,000 and
AoA=8°
.

3

LES
Galbraith overset 6th-order

2.5
2

-CP

1.5
1
0.5
0
-0.5
-1
0

0.2

0.4

0.6
0.8
1
X/C
Figure 30: Time-averaged surface pressure coefficient for the flow with steady freestream,
Re=60,000 and AoA=8°
.

45

0.03
LES
Galbraith overset 6th-order
0.02

CF

0.01

0

-0.01

-0.02
0

0.2

0.4

0.6

0.8

1

X/C
Figure 31: Time-averaged skin friction coefficient on the airfoil upper surface for the flow with
steady freestream, Re=60,000 and AoA=8°
.

To consider a much more extreme case in terms of flow separation, the 3D flow over the
airfoil for AoA=16°and Re=60,000 is also computed. The time- and spanwise-averaged data for
this case are presented in Table 5. For comparison, the results for AoA=4° 8° and 12°are also
, ,
listed. The instantaneous vorticity iso-surfaces for AoA=16°in Figure 27 clearly show the rapid
separation of the flow at the leading edge. At this steep AoA, the flow is not able to recover the
static pressure (as shown in the CP plot in Figure 32) and thus cannot reattach to the airfoil
surface. The AFRL measurements [9] show that the maximum CL for SD7003 airfoil at this
Reynolds number is around 1.05, and CL drops to 0.96 at AoA=14° Following this trend, the
.
predicted CL=0.91 at AoA=16°seems reasonable. Since experimental and numerical data for the
simulated Re and AoA are not available for direct comparison, the time-averaged pressure
coefficient is compared with Galbraith [7] results at AoA=14° Re=60,000 in Figure 32. Also
,
shown in this figure are the LES results for the same Reynolds number but lower AoAs of 4°and

46

8° The difference between LES results at different AoA is much more significant on the upper
.
side of the airfoil as expected. The overall trends observed in current LES results are similar to
that of Galbraith, though minor differences exist.

AoA (°
)
4
8
12
16

Separation Transition
Xs C
Xt C
0.25
0.50
0.042
0.16
0.008
0.006
-

Reattachment
Xr C
0.64
0.25
0.92
-

Lift Coeff.
CL
0.59
0.93
1.14
0.91

Drag Coeff.
CD
0.020
0.040
0.124
0.277

Table 5: Time-averaged data for the flow with steady freestream, Re=60,000, and AoA=4° 8°
, ,
12° 16°
,
.

3
AoA=4, LES
AoA=8, LES
AoA=16, LES
AoA=14, Galbraith overset 6th-order

2.5
2

-CP

1.5
1
0.5
0
-0.5
-1
0

0.2

0.4

0.6
0.8
1
X/C
Figure 32: Time-averaged surface pressure coefficient for the flow with steady freestream at
different AoAs and Re=60,000.

Katz et al. [11] at MSU’s TMUAL conducted a series of water channel measurements for
steady freestream flows around a fixed SD7003 airfoil using Molecular Tagging Velocimetry

47

(MTV) method. The reported freestream turbulence intensity in the water channel was about
2.8% of the mean freestream velocity. Also, the Reynolds number in this experiment varied from
20,000 to 40,000 and AoA varied from 3.75°to 10.44° These experiments are simulated here
.
via LES. Table 6 compares the LES predictions and the experimental data. There are some
differences between the numerical and experimental results. At relatively small AoAs, LES
predicts earlier separation in comparison to experiment. The prediction of reattachment point is
generally not as good as that of the separation point. Some of these differences are inevitably due
to modeling and numerical issues. There could be some issues related to the experiment. One
possible reason is the relatively high turbulence intensity in the water channel. Recent
measurements by Olson [12] show that freestream turbulence will move the separation point
downstream and move the reattachment point upstream. More discussion regarding the effect of
freestream unsteadiness on mean separation and reattachment points is given in Section 5.3. Katz
et al. [11] also reported uncertainty of ±
0.4° in the AoA measurement. This may be another
reason for the discrepancy between numerical and experimental results.

Reynolds Number

20,000 30,000

40,000

AoA (°
)

3.75

7.39

5.49

10.44

Separation
LES
( Xs )
TMUAL
Reattachment
LES
( Xr )
TMUAL

0.303
0.402
0.981
0.994

0.063
0.054
0.661
0.438

0.121
0.223
0.652
0.526

0.017
0.024
0.102
0.144

Table 6: Separation and reattachment points obtained by LES and TMUAL experiments.

The lift and drag coefficients by LES for all steady freestream cases are compared in Figure
33. CL and CD for the cases with Re=60,000 have been compared to the available experimental
and numerical results. For other Reynolds numbers and AoAs, we found no numerical or
48

experimental data for CL and CD in the literature. A rough comparison at data of different
Reynolds numbers suggest that, for the same AoA, lift coefficient increases as Reynolds number
increases, and drag coefficient decreases as Reynolds number increases. The same trend was
observed in simulation results of Galbraith [20] at AoA=4°and 8°and Re=10,000 to 90,000.

CL, Re=60,000
CL, Re=40,000

1.3

0.35

CL, Re=30,000

1.2

CL, Re=20,000

0.3

CD, Re=60,000

1.1

CD, Re=40,000

0.25

CD, Re=30,000

1

CD, Re=20,000

0.2

0.8

0.15

C

CL

D

0.9

0.7
0.1
0.6
0.05

0.5
0.4

2

4

6

8

10
12
AoA (  )

14

16

0
18

Figure 33: Time-averaged lift and drag coefficients obtained by LES for flows with steady
freestream condition.

49

5.

Large Eddy Simulation of Unsteady Freestream Flows over SD7003 Airfoil

In this chapter, 3D flows with unsteady freestream condition over SD7003 are computed by
LES with the WALE model. The effects of flow unsteadiness on the aerodynamic force,
boundary layer separation/reattachment, vorticity and velocity fluctuations are systematically
studied. The unsteadiness in freestream flow is assumed to be simple and sinusoidal. Two sets of
unsteady flows are considered. In Section 5.1, the freestream flow direction or AoA is fixed,
while the freestream velocity magnitude or speed is varied harmonically in time. In Section 5.2,
the freestream velocity magnitude is fixed, while the freestream flow direction or AoA is
changed harmonically in time. The effect of freestream flow unsteadiness on boundary layer
separation and reattachment is discussed in Section 5.3. Section 5.4 considers cases involving a
wind gust.
5.1

Flows with Oscillating Freestream Velocity Magnitude
In this set of unsteady flows, the airfoil is fixed at 4° relative to the freestream flow

direction or x-axis, so the AoA is fixed at 4°but the velocity magnitude oscillates sinusoidally in
time according to Equation (37) (Figure 13). In this equation, U m is the mean velocity
magnitude, σ denotes the amplitude of freestream velocity oscillation normalized by the mean
velocity, and   2 f is the angular frequency. The reduced frequency is k   f C U m and the
mean Reynolds number based on U m is 60,000.
 4


 U  t   U m 1   sin  t  

 V W  0


50

(37)

Case

σ

k

CL
CL
CD
CD
Mean Amplitude Mean Amplitude
-3

Steady
mean
freestream

0

-

U1

0.183

U2

XS C
Mean

XR C
Mean

0.25

0.64

-4

0.59

2.6×10
(standard
deviation)

0.020

2.5×10
(standard
deviation)

2

0.56

0.43

0.021

0.485

0.33

0.53

0.183



0.57

0.25

0.020

0.243

0.33

0.51

U3

0.183

 2

0.57

0.18

0.020

0.123

0.33

0.46

U4

0.183

 4

0.57

0.17

0.020

0.063

0.33

0.44

U5

0.183

 8

0.57

0.16

0.020

0.033

0.32

0.47

U6

0.366

2

0.60

0.86

0.021

0.968

0.30

0.53

U7

0.366



0.60

0.49

0.020

0.489

0.32

0.49

U8

0.366

 2

0.60

0.36

0.020

0.248

0.31

0.35

U9

0.366

 4

0.60

0.34

0.020

0.126

0.31

0.37

U10

0.366

 8

0.60

0.33

0.020

0.066

0.32

0.40

Table 7: Mean quantities for cases with oscillating velocity magnitude.

Table 7 lists the flow parameters for ten simulated cases with unsteady freestream speed
(labeled U1 to U10). In cases U1-U5 and cases U6-U10 the amplitude of oscillations is the same
but the frequencies are different. The amplitudes and frequencies chosen cover a range of these
parameters limited by the LES grid. The force coefficients discussed below are defined in
Section 3.5.
5.1.1 Aerodynamic forces
Figures 34-37 show the variation of CL and CD during a cycle after the flows have fully
developed under the periodic freestream condition. In these figures, a negative CD value implies
a net thrust force on the airfoil. Both CL and CD oscillate periodically with the frequency of the
freestream velocity, and a phase shift with respect to freestream velocity which is a function of
the reduced frequency and freestream velocity amplitude. The mean and amplitude of lift and
51

drag coefficients are listed in Table 7 along with data for the steady freestream flow with
AoA=4°and Re=60,000. Low values of the standard deviation of CL and CD in Table 7 for the
steady flow indicate that lift and drag forces are constant when there is no oscillation in
freestream velocity. Interestingly, Table 7 also show that the freestream velocity oscillation has
little or no effect on the mean CL and CD which stay close to their steady freestream values. A
comparison between cases U1~U5 and cases U6~U10 (five pairs: U1-U6, U2-U7, …, U5-U10)
shows that the mean lift increases slightly (around 5%) when the freestream amplitude is doubled.
CL and CD in Figures 34-37 are normalized by the mean freestream velocity U m as
explained in Section 3.5. If the instantaneous freestream velocity U(t) is used for the
normalization,

a new set of lift and drag coefficients can be defined/computed as

*
*
CL  L ( U (t )2 C 2) and CD  D ( U (t )2 C 2) which are plotted in Figures 38-39 for cases

U6-10. Evidentially, the time profiles of CL* and CD* are significantly different from those of CL
and CD and are no longer “sinusoidal-like”.
Figure 40 shows amplitudes of fluctuations in CL and CD as functions of reduced frequency
k and the amplitude of freestream velocity oscillation σ. Evidently, CD amplitude increases
linearly with the reduced frequency. However, CL is not fully linear with respect to the reduced
frequency. At relatively low reduced frequency, CL amplitude increases slowly with increasing
reduced frequency. At higher reduced frequency, CL amplitude increases faster and its relation
with the reduced frequency seems to become closer to being linear. At any reduced frequencies
considered in this study, doubling the amplitude of freestream velocity oscillation causes both CL
and CD amplitudes to double, suggesting a linear relation between the amplitude of freestream
velocity oscillation and the aerodynamic forces. Figure 40 also shows comparison with the
Greenberg theory [53] which will be discussed later.

52

1.6

U/Um, =0.183

1.5

CL, =0.183, k=

1.2

CL, =0.183, k=/4

0.8

1.3

CL, =0.183, k=/2

1.2

1.4

1.1

CL

CL, =0.183, k=/8

0.4

U/Um

CL, =0.183, k=2

1
0.9
0.8

0

0.7
0.6

0.5
180
225
270
315
360
t (deg.)
Figure 34: Time variations CL for flows with oscillating freestream speed (cases U1-5).
0

45

90

135

U/Um, =0.183

1.2

1.5

CD, =0.183, k=2
CD, =0.183, k=

1.2

CD, =0.183, k=/4

0.4

1.3

CD, =0.183, k=/2

0.8

1.4

1.1

C

D

CD, =0.183, k=/8

0

1
0.9
0.8

-0.4

0.7
0.6

-0.8
-1

0.5
180
225
270
315
360
t (deg.)
Figure 35: Time variations of CD for flows with oscillating freestream speed (cases U1-5).
0

45

90

135

53

U/Um

-0.4

U/Um, =0.366

1.4

CL, =0.366, k=

1.3

CL, =0.366, k=/2

1.2

1.5

CL, =0.366, k=2

CL phase shift, case U6

1.2

CL, =0.366, k=/4

0.8
CL

CL, =0.366, k=/8

1.1
U/Um

1.6

1

0.4

0.9
0.8

0

0.7
0.6

0.5
180
225
270
315
360
t (deg.)
Figure 36: Time variations of CL for flows with oscillating freestream speed (cases U6-10). The
phase shift for one of the cases is shown as an example.
0

45

90

135

U/Um, =0.366

1.5

CD, =0.366, k=2

1.2

1.4

CD, =0.366, k=
CD, =0.366, k=/2

CD

1.2

CD, =0.366, k=/8

0.4

1.3

CD, =0.366, k=/4

0.8

1.1
1

0

0.9
0.8

-0.4

0.7
0.6

-0.8
-1

0.5
180
225
270
315
360
t (deg.)
Figure 37: Time variations of CD for flows with oscillating freestream speed (cases U6-10).
0

45

90

135

54

U/Um

-0.4

U/Um, =0.366

3

C* , =0.366, k=2
L
C* , =0.366, k=
L

1.1

*

1
1

U/U

CL

m

1.2

C* , =0.366, k=/8
L

1.5

1.3

C* , =0.366, k=/4
L

2

1.4

C* , =0.366, k=/2
L

2.5

1.5

0.9
0.8

0.5

0.7
0
-0.5

0.6

0.5
180
225
270
315
360
t (deg.)
*
Figure 38: Time variations of CL for flows with oscillating freestream speed (cases U6-10).
0

45

90

135

U/Um, =0.366
C* , =0.366, k=2
D
C* , =0.366, k=
D

1.3

C* , =0.366, k=/4
D

1

1.4

C* , =0.366, k=/2
D

1.6

1.5

1.2

C* ,
D

0.4

=0.366, k=/8

1.1

*

CD

1
0.9

-0.2
0.8
0.7

-0.8

0.6
-1.4

0.5
180
225
270
315
360
t (deg.)
*
Figure 39: Time variations of CD for flows with oscillating freestream speed (cases U6-10).
0

45

90

135

55

U/Um

2

1

CL,  = 0.183

0.9

CD,  = 0.183
CL,  = 0.366

0.7

Amplitude

0.8

CD,  = 0.366

0.6

CL,  = 0.183 (Greenberg)
CL,  = 0.366 (Greenberg)

0.5
0.4
0.3
0.2
0.1
0

0

1

2

3

4

5

6

7

k

Figure 40: Amplitude of CL and CD for flows with oscillating freestream velocity magnitude.
Figures 34-37 show that the phase of aerodynamic forces leads those of the freestream
velocity. This is shown better in Figure 41 which plots the phase shift between the freestream
velocity oscillations and lift/drag oscillations for different reduced frequency and amplitude. The
phase shift is defined as the phase difference between the maximum lift (or drag) and the
maximum freestream velocity. Unlike amplitudes of CL and CD which are significantly affected
by σ, Figure 41 shows that the phase shift is slightly affected by σ and is mainly dependent on
the reduced frequency. For both lift and drag, the phase shift increases as reduced frequency
increases, more so at lower frequencies. For the considered range of reduced frequency (π/8 ~
2π), CD phase shift with respect to freestream velocity varies between 70 and 90 degrees. At high
reduced frequencies, CD phase shift seems to asymptotically approach 90 degrees, indicating that
the maximum drag coincides with the maximum acceleration of freestream velocity (instead of
the maximum velocity magnitude). For the same range of reduced frequency, CL phase shift has
56

a much wider range with the maximum value being around 70 degrees and a minimum close to
zero. However, the rate of increase in phase shift decreases at higher reduced frequencies. At low
freestream frequency, CL phase shift quickly drops to zero.

100
90

Phase Shift (deg.)

80
70
60
50

CL,  = 0.183

40

CD,  = 0.183

30

CL,  = 0.366

20

CD,  = 0.366

10

CL,  = 0.183 (Greenberg)

0

CL,  = 0.366 (Greenberg)

0

1

2

3

4

5

6

7

k

Figure 41: Phase shift between force coefficients and freestream velocity for flows with
oscillating freestream velocity magnitude.

The LES predictions of the CL amplitude and its phase shift are compared with the inviscid
Greenberg’s theory [53] in Figures 40 and 41. According to Greenberg’s theory, the lift
coefficient is obtained from Equation (38), where CL, m is the lift coefficient of the steady
freestream flow obtained by the thin airfoil theory [54] as CL, m  2 .
  2F 
CL
k

 1 
     G  cos  t   1  F  sin  t


CL, m 
2 
2



 F
2

2

cos 2 t 

 G
2

2

57

sin 2 t

(38)

Variables F and G in Equation (38) are the real and imaginary parts of the Theodorsen’s function
[26] both functions of the reduced frequency k. In Figure 40, the two Greenberg CL amplitude
curves have different rates at the same reduced frequency since σ appears as an independent
variable in Equation (38). For CL amplitude, LES values are higher than those predicted by
Greenberg’s theory. However, the CL phase shift obtained from LES data compares much more
favorably with Greenberg’s theory. The inviscid theory predicts that the harmonic freestream
velocity oscillations do not affect the mean CL which still takes the steady freestream flow value
of 0.44. Indeed, the mean CL computed by LES for all ten unsteady cases match the steady flow
CL, even though the LES values are generally higher than those of the inviscid theory.
Greenberg’s theory is based on the assumptions that the flow is two-dimensional, fully attached
over the airfoil, and the wake is planar. These assumptions are not completely true for the 3D
flows studied here, suggesting why the results predicted by the 3D LES do not match very well
with the inviscid theory.
At high reduced frequencies, Figure 41 indicates that the maximum CL values coincide with
the maximum freestream velocity acceleration. To explain this behavior, the airfoil surface
pressure coefficient CP at the time of maximum acceleration of freestream velocity for two
“high” frequency cases U6 and U7 is shown in Figure 42. The pressure component of the lift
coefficient is obtained directly from the integral of CP over the airfoil surface. For the two cases
considered in Figure 42, the pressure component makes up almost the entire lift and the viscous
component is negligible. Our LES results (Figure 43), as expected, indicate that the surface
pressure is strongly affected by the freestream pressure. In the freestream, both the y- and zvelocities are zero and the x-velocity is spatially uniform. When gravity is further neglected, the
momentum equations are reduced to Equation (39) and then Equation (40) for the x-velocity.

58

Equation (40) indicates that the freestream pressure gradient in x-direction is proportional to the
freestream acceleration and at higher frequencies higher pressure gradient around the airfoil is
expected to develop. This will lead to higher pressure difference between the lower and upper
airfoil surfaces, and thus high amplitude of the lift coefficient.

5
Lower airfoil surface

k=2 (Case U6)
k= (Case U7)

0

CP

-5
-10
-15
-20

Upper airfoil surface
0

0.2

0.4

0.6

0.8

1

X/C
Figure 42: Pressure coefficient on the airfoil surface at the time of maximum acceleration of
freestream velocity for case U6 with k  2 and   0.366 and case U7 with k   and
  0.366 .

CP
8
4
0
-4
-8
-12
-16
-20
Figure 43: Pressure coefficient around the airfoil at the time of maximum acceleration of
freestream velocity for case U6 with k  2 and   0.366 .
59

  2u  2u  2u 
 u
u
u
u 
p
 u  v  w       2  2  2    gx
 x
x
y
z 
x
x
x 
 t





p
u
 
x
t

(39)

(40)

5.1.2 Vorticity field
The instantaneous vorticity fields for flows with oscillating freestream velocity magnitude
are compared to those of the steady freestream flow in Figures 45-59. For each flow, the 3D isosurfaces of the z-vorticity component, the 2D contours of spanwise-averaged z-vorticity
component, and the line integrals of the spanwise-averaged vorticity magnitude are shown. All
iso-surfaces are set to have the same normalized z-vorticity value of z C U m  10 . A line is
drawn across the four snapshots for a better comparison. It should be noted that the vortices
formed far beyond the trailing edge are not accurately computed because of the relatively coarse
mesh used in regions away from the airfoil. To quantify the vorticity field, the normalized
spanwise-averaged vorticity magnitude  

C
2
2
2
x   y  z is integrated along the airfoil
Um

surface normal and tangent direction (n and s axes in Figure 44, n and s are normalized by the
airfoil chord C). The spanwise-averaged vorticity magnitude which is an indication of vorticity
strength is integrated along the n axis from the airfoil surface to the far-field boundary of the
computational domain, i.e.

  dn . The spike in   dn close to the trailing edge is due to rapid

and local formation of the wake at the trailing edge. The profile of

  dn

oscillates

considerably in the recirculation region, which makes it difficult to compare vorticity values at
different times. To get a better picture of how vorticity magnitude changes over time,

60

  dn is

further integrated along the s axis from the leading edge to the trailing edge on the upper surface
as

   dnds .
Figures 45-47 show that the vorticity field over the airfoil does not significantly change in

time in the steady freestream flow, despite significant variations of the vorticity and vorticity
magnitude over the airfoil. Even though there are considerable fluctuations in vorticity
magnitude in the recirculation region, the line integral of the spanwise-averaged vorticity
magnitude stays almost constant in time. For the flows with oscillating freestream velocity
magnitude, the vorticity field and magnitude vary much more significantly in time. It appears in
the iso-surface plots of vorticity that the separation point moves back and forth during each cycle
depending on the reduced frequency and amplitude of the freestream velocity oscillations. At
high reduced frequencies (Figures 49 and 55), the fast oscillating freestream flow seems to
promote the formation of distinct vortices which convect downstream and break down to smaller
vortices near the trailing edge (as shown by arrows in Figure 55). This is not seen in flows with
low reduced frequencies (Figures 52 and 58). The integrated vorticity magnitude (   dn or

   dnds ) clearly vary in time much more than that observed for steady freestream flow. For
all cases with oscillating freestream velocity magnitude,

   dnds values correlate well with

the freestream velocity magnitude, and for the time instants shown in Figures 50, 53, 56 and 59,
the maximum (minimum) value of

   dnds

coincides with the maximum (minimum)

freestream velocity magnitude. At t0 and t0+T/2, the freestream velocity magnitudes are the
same, and values of the integral

   dnds are relatively close; the small difference could be due

to possible phase shift between the maximum

   dnds and the maximum freestream velocity.

61

A direct comparison between unsteady cases (cases U1 U5, U6 and U10) in Figure 60 indicates
that the integrated vorticity values (   dn ) on the airfoil upper surface for different k and σ
start to show fluctuations around X/C=0.3. It will be shown in Section 5.1.4 that X/C=0.3 is
1

close to the mean separation point for all unsteady cases listed in Table 7. Clearly and expectedly,
the vorticity magnitude and the deviation between different unsteady flows are significant at or
0.8

after the separation point. At later phases, e.g. at t0+3T/4, the vorticity magnitude for different
cases with different k and σ start to deviate from each other much sooner at locations closer to
0.6
the leading edge. This can be perhaps related to the shift of the separation point upstream toward
the leading edge.
0.4

0.2

n
s

0

Figure 44: The n-s coordinate system around the airfoil surface.
-0.2

0

0.2

0.4

0.6

62

0.8

1

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 45: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one residence time
for the flow with steady freestream, Re=60,000, AoA=4° (T is the residence time C U m , same
.
for Figures 46-47).

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36

t0+3T/4

-40

Figure 46: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one
residence time for the flow with steady freestream, Re=60,000, AoA=4°
.

63

3

||dn

2

(a)

1.5
1

1.8
1.7
1.6

0.5
0
0

(b)

1.9

||dnds

2.5

2

t0
t0+T/4
t0+T/2
t0+3T/4

0.2

0.4

0.6

0.8

1.5
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 47: Integrals of spanwise-averaged vorticity magnitude |Ω| for the flow with steady
freestream, Re=60,000, AoA=4° (a) Integrated along n axis; (b) integrated along n and then s
.
axes.

t0

u/Um
2.00

t0+T/4

1.75
1.38

t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 48: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for case
U1 with k  2 and   0.183 . (t0 is the beginning of the cycle when U equals Um and is
increasing. T is the period. Same for Figures 49-60)

64

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

Figure 49: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one cycle
for case U1 with k  2 and   0.183 .

3

2.2
t0
t0+T/4
t0+T/2
t0+3T/4

||dn

2

(a)

(b)
2

||dnds

2.5

1.5
1

1.6
1.4

0.5
0
0

1.8

0.2

0.4

0.6

0.8

1.2
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 50: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U1 with k  2 and
  0.183 . (a) Integrated along n axis; (b) integrated along n and then s axes.

65

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 51: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for case
U5 with k   8 and   0.183 .

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

Figure 52: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one cycle
for case U5 with k   8 and   0.183 .

66

3

2.2
t0
t0+T/4
t0+T/2
t0+3T/4

2
1.5
1

1.8
1.6
1.4

0.5
0
0

(b)

2

||dnds

2.5

||dn

(a)

0.2

0.4

0.6

0.8

1.2
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 53: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U5 with k   8 and
  0.183 . (a) Integrated along n axis; (b) integrated along n and then s axes.

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 54: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for case
U6 with k  2 and   0.366 .

67

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

Figure 55: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one cycle
for case U6 with k  2 and   0.366 .

3

2

(a)

(b)
2

||dnds

2.5

||dn

2.2

t0
t0+T/4
t0+T/2
t0+3T/4

1.5
1

1.6
1.4

0.5
0
0

1.8

0.2

0.4

0.6

0.8

1.2
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 56: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U6 with k  2 and
  0.366 . (a) Integrated along n axis; (b) integrated along n and then s axes.

68

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 57: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) during one cycle for case
U10 with k   8 and   0.366 .

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

Figure 58: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) during one cycle
for case U10 with k   8 and   0.366 .

69

3

2.2
t0
t0+T/4
t0+T/2
t0+3T/4

||dn

2

(a)

(b)
2

||dnds

2.5

1.5
1

1.6
1.4

0.5
0
0

1.8

0.2

0.4

0.6

0.8

1.2
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 59: Integrals of spanwise-averaged vorticity magnitude |Ω| for case U10 with k   8 and
  0.366 . (a) Integrated along n axis; (b) integrated along n and then s axes.

70

3
2

=0.183, k=2
=0.183, k=/8
=0.366, k=2
=0.366, k=/8

t0
1
0
3
2

t0+T/4

||dn

1
0
3
2

t0+T/2
1
0
3
2

t0+3T/4

1
0
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
X/C

1

Figure 60: Integral of spanwise-averaged vorticity magnitude |Ω| along n axis for cases U1, U5,
U6, and U10.
5.1.3 Velocity fluctuations
The effect of freestream velocity unsteadiness on the velocity fluctuations over the airfoil is
studied in this section by comparing the root mean square (RMS) values of three velocity
components and turbulent kinetic energy (TKE) of steady and unsteady cases.

71

Velocity RMS is calculated for unsteady case U6 with σ=0.366 and k=2π, using LES data
for five freestream cycles. This case is chosen due to its short period and large amplitude of
oscillation in freestream velocity magnitude. The larger amplitude of σ=0.366 is chosen because
it is expected that larger amplitude of oscillation in freestream velocity causes larger variation in
velocity RMS. During each cycle, the velocity data are recorded at four different times or phases.
The 3D computational domain has 101 grids in the spanwise (z) direction. Therefore, for a given
point on the x-y plane and five available cycles, there are 505 (=101× instantaneous velocity
5)
samples for spanwise and phase averaging, from which the velocity RMS and the turbulence
kinetic energy (TKE) are calculated. The velocity RMS of the steady freestream flow with
Re=60,000 and AoA=4° is also calculated for comparison, using velocity data for five flow
residence times. Figure 61 shows the integrated values of the RMS of three velocity components
for the steady case. These integrated values are easier to interpret and are obtained by integrating
the velocity RMS and TKE values along the n or n and s directions (Figure 44). These integrated
values are represented by

  RMS dn ,    RMS dnds , and   TKE dnds . The integrated values

of RMS of u, v, w, and TKE along airfoil surface normal and tangent directions (n and s), i.e.

   RMS dnds ,   TKE dnds ,

are 0.0237, 0.0209, 0.0024 and 0.0004, respectively for the

steady case. Figure 61 shows that, for the steady freestream flow,

  RMS dn increases along the

airfoil surface and peaks around the trailing edge.
The integrated values of velocity RMS and TKE for the unsteady case U6 are plotted in
Figures 62-64. Similar to steady case, Figure 62 shows that the integrated RMS values of all
three velocity components are essentially zero at X/C=0, but gradually increases till trailing edge,
where it suddenly increases. Figures 62-63 show that the RMS of velocity fluctuations changes
only slightly with the phase and are similar to those shown in Figure 61 for steady flow. Figure
72

64 shows that the TKE values are slightly higher at time t0+T/4 compared to other times or
phases, suggesting that higher freestream velocity induces stronger “turbulence” and higher
velocity fluctuations. Also similar to steady flow, it is shown that the fluctuations or RMS of xand y-velocities are comparable, but much bigger than that of z-velocity. Naturally, the TKE is
dominated by the x- and y-velocity fluctuations and behaves similar to RMS of these velocity
components.
By comparing the integrals of velocity RMS and TKE along n and s axes, i.e.

   RMS dnds and   TKE dnds in Figures 63 and 64, it is found that the x- and y-velocity
fluctuations of the unsteady case U6 with k = 2π and σ = 0.366 are generally lower than those of
the steady freestream flow while the z-velocity fluctuations are close for these two flows. This is
also shown in Figure 65 for   RMS dn . It seems that the oscillation in freestream velocity
magnitude suppresses the fluctuations of the x- and y-velocities but does not noticeably affect the
spanwise velocity fluctuations.
0.35
0.3

(RMS)dn

0.25

u'/Um
v'/Um
w'/Um

0.2
0.15
0.1
0.05
0
0

0.2

0.4

0.6

0.8

1

X/C
Figure 61: Integrated values of velocity RMS along n axis for the flow with steady freestream,
Re=60,000 and AoA=4°
.

73

(RMS)dn

0.18
0.12
0.06

0.04

u'/Um

0.03

t0

v'/Um

0.02

w'/Um

0.01

(RMS)dn

0
0.18

0
0.04
0.03

t0+T/4

0.12

0.02

0.06

0.01

(RMS)dn

0
0.18

0
0.04
0.03

t0+T/2

0.12

0.02

0.06

0.01

(RMS)dn

0
0.18

0
0.04
0.03

t0+3T/4

0.12

0.02

0.06
0
0

0.01

0.5
X/C

0
0

1

0.5
X/C

1

Figure 62: Integrated values of velocity RMS along n axis for case U6 with k  2 and
  0.366 .

74

(RMS)dnds

0.015

0.01
u'/Um
v'/Um

0.005

w'/Um

0
t0

t0+T/4

t0+T/2
t0+3T/4
t0+T
Time
Figure 63: Integrated values of velocity RMS along n and s axes for case U6 with k  2 and
  0.366 .

-4

4

x 10

(TKE)dnds

3

2

1

0
t0

t0+T/4

t0+T/2
t0+3T/4
t0+T
Time
Figure 64: Integrated values of turbulence kinetic energy (TKE) along n and s axes for case U6
with k  2 and   0.366 .

75

(RMS)dn

0.04
0.03
0.02

steady freestream
t0
(a) u'/Um
t0+T/4
t0+T/2
t0+3T/4

0.01
0
0.04

(RMS)dn

0.03

(b) v'/Um

0.02
0.01
0
0.01

(RMS)dn

0.008
(c) w'/Um

0.006
0.004
0.002
0
0

0.2

0.4

0.6

0.8

1

X/C

Figure 65: Comparison of integrated values of velocity RMS along n axis between case U6 with
k  2 and   0.366 and the flow with steady freestream with Re=60,000 and AoA=4° (a) u
.
RMS; (b) v RMS; (c) w RMS.

5.1.4 Boundary layer separation and reattachment
Effect of oscillating freestream velocity on the flow over airfoil is further examined in this
section by comparing the separation and reattachment points for various unsteady and steady
cases. Similar to the separation and reattachment points obtained from the mean flow field, the
instantaneous separation point can be calculated to be the location where the instantaneous skin
friction coefficient CF first goes from positive to negative, and the instantaneous reattachment
76

point is defined as the location where CF last goes from negative to positive. The mean CF is
calculated from the time- and spanwise-averaged velocity, and the instantaneous CF is based on
the instantaneous spanwise-averaged velocity. For the steady freestream flow, Figure 66 clearly
shows that the instantaneous CF profiles at different times closely follow the mean profile. The
instantaneous profiles are identical to the mean profile before the separation point, and
consequently, the instantaneous separation and reattachment points oscillate around the mean
values with small amplitude. The small oscillations in the CF are due to small perturbations in the
vortical flow over the airfoil and are not due to any freestream turbulence in the flow which is
negligible. When the freestream velocity magnitude oscillates at low reduced frequency of π/8
(Figure 67), the instantaneous separation and reattachment points fluctuate with amplitudes
larger than those of the steady freestream flow. At higher reduced frequency of 2π (Figure 68),
the four instantaneous CF profiles deviate further from the mean profile, so much so that the flow
at time t0 and t0+T/4 becomes almost fully attached. Compared to Figure 68, Figure 69 shows
that higher freestream velocity amplitude causes even wider oscillation in the skin friction
coefficient profile. The CF profile at t0+T/4 in Figure 69 suggests that the flow is fully attached,
while the CF profiles at t0+T/2 and t0+3T/4 show the flow, after early separation, to reattach as
late as the trailing edge. Figures 67-69 generally show that the separation will be delayed or
completely eliminated as freestream velocity is accelerating, and it will happen earlier and closer
to the leading edge when the freestream velocity is decelerating. This is related to the freestream
pressure gradient as Equation (40) suggests. As the freestream velocity is accelerating, the strong
negative pressure gradient in the freestream flow will alleviate or completely eliminate the
positive pressure gradient induced by the airfoil curvature which causes the separation.

77

The mean separation and reattachment points computed based on the mean CF curves for
steady and unsteady cases are listed in Table 7. The mean separation points for all ten unsteady
cases are shown to be quite similar, while the mean reattachment points vary between 0.35 and
0.53 without a clear trend. Compared to the flow with steady freestream velocity, each of
unsteady cases has a delayed mean separation point and an earlier mean reattachment point
closer to the airfoil leading edge.

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 66: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for the flow with steady freestream, Re=60,000 and AoA=4°
.

78

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 67: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for case U5 with k   8 and   0.183 .

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 68: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for case U1 with k  2 and   0.183 .

79

0.03

Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 69: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for case U6 with k  2 and   0.366 .

5.2

Flows with Oscillating Angle of Attack (AoA)
In this section we discuss the results for a set of simulations with variable AoA. The airfoil

is assumed to be fixed but the direction of freestream flow is periodically changed according to
Equation (41) below. The freestream velocity has a constant magnitude of U m with a mean
Reynolds number of 60,000. The mean AoA is fixed at  m  4 , while the amplitude and
reduced frequency,  f and k, are changed within a range similar to that used for velocity
magnitude oscillations in Section 5.1. The ten cases with unsteady AoA considered in this
section are listed in Table 8. The range of variation in the reduced frequency, k, is chosen again
to be the same as that in Table 7. The amplitude of AoA oscillations is 4°for cases A1 through
A5 and 8°for cases A6 through A10.
80

   t    m   f sin t 

 U  t   U m cos    m 

 V  t   U m sin    m 
W 0


Case

f
(deg.)

k

CL
CL
CD
CD
Mean Amplitude Mean Amplitude
-3

Steady
mean
freestream

0

-

A1

4

A2

(41)

XS C
Mean

XR C
Mean

0.25

0.64

-4

0.59

2.6×10
(standard
deviation)

0.020

2.5×10
(standard
deviation)

2

0.57

2.30

0.018

0.158

-

-

4



0.56

1.13

0.017

0.080

0.35

0.42

A3

4

 2

0.56

0.57

0.017

0.044

0.33

0.46

A4

4

 4

0.56

0.31

0.016

0.030

0.30

0.42

A5

4

 8

0.56

0.25

0.015

0.024

0.29

0.37

A6

8

2

0.60

4.62

0.011

0.322

-

-

A7

8



0.58

2.29

0.008

0.160

-

-

A8

8

 2

0.57

1.15

0.007

0.089

-

-

A9

8

 4

0.56

0.62

0.006

0.062

-

-

A10

8

 8

0.56

0.50

0.003

0.051

-

-

Table 8: Time-averaged data for cases with oscillating AoA.

Lift coefficient, drag coefficient, pressure coefficient, and skin friction coefficient are
defined in Section 3.5. As the AoA oscillates in these flows, it should be noted that the lift and
drag forces can be computed by decomposing the aerodynamic force vector along the coordinate
directions or axes constructed based on fixed mean AoA of  m or they can be obtained by
decomposing the force along coordinate directions constructed based on instantaneous AoA or

 (t ) (Figure 13). Both of these will be considered below.

81

5.2.1 Aerodynamic forces
Figures 70-73 show the variations of lift and drag coefficients (CL and CD) for one cycle
after the flows have fully developed under the periodically variable freestream flow AoA
condition. Lift and drag forces and CL and CD in these figures are calculated over the fixed (mean)
flow axial and normal directions. Similar to what was observed in flows with oscillating
freestream velocity magnitude, CL and CD oscillate in time with the freestream frequency, and CL
and CD profiles are generally “sinusoidal-like”.
The lift and drag forces can also be computed by decomposing the aerodynamic force
vector along the coordinate directions constructed based on instantaneous AoA or  (t ) . The new
lift and drag coefficients, denoted as CL*, CD*, can be calculated by a simple conversion from CL
and CD as shown in Equation (42), and the results for cases A6-10 are shown in Figures 74 and
75. Clearly, CL* is similar to CL, but CD* is significantly different from CD, simply because the
conversion affects the smaller drag coefficient much more than it does the lift coefficient. As
shown in Figure 75, CD* profiles are no longer “sinusoidal-like”.
*
 CL  CL cos  (t )   m   CD sin  (t )   m 

 *
CD  CD cos  (t )   m   CL sin  (t )   m 


(42)

The mean and amplitude of CL and CD are listed in Table 8 along with the data for the
steady freestream flow with Re=60,000 and AoA=4° It is clear that reduced frequency and
.
amplitude of AoA oscillation have only slight effect on the mean lift coefficient which stays
close to its steady freestream flow value of 0.59. However, in all cases, the mean drag coefficient
is lower than its steady freestream flow value of 0.020. The reduction in drag of an airfoil when
the freestream direction is oscillating is referred to as the Katzmayr effect [33]. Toussaint et al.

82

[34] verified the existence of Katzmayr effect in a set of wind tunnel experiments with
oscillating freestream flow direction with respect to an airfoil, induced by pitching blades. Based
on these experiments, Ober [55] provided an explanation for the Katzmayr effect, showing that,
on average, the rotated lift vector has a bigger component in negative drag (thrust) than in
positive drag. In the experiments by Toussaint et al., the reduced frequencies are lower than 0.5,
and effects of reduced frequency and amplitude of the AoA oscillation cannot be separated
because these two parameters could not be changed independently in the experiments. The LES
results here confirm the Katzmayr effect to a much bigger range of reduced frequency and also
show that the drag reduction is increased at lower reduced frequency and higher amplitude of
AoA oscillation. For example, CD drops as much as 83% from case A1 to case A10 because  f
is doubled and k is reduced to 1/16 of its value in case A1. The effect of AoA oscillation
amplitude on drag reduction may be examined by considering case A10. In case A10, the
amplitude of oscillation in AoA is 8° allowing the flow AoA to vary between -4°and 12°during
,
one cycle. For this case, the time interval with negative AoA is 1/3 of the cycle. At a negative
AoA, the freestream velocity has a component in the negative y-direction (Figure 13), resulting
in a thrust on the airfoil. By comparing CD profiles of cases A5 (Figure 71) and A10 (Figure 73),
it is clear that the airfoil spends a bigger portion of each period of CD oscillations with negative
drag (thrust) when the freestream AoA amplitude is higher.
The amplitudes of oscillations in CL and CD for the unsteady AoA cases are listed in Table 8
and are also plotted in Figure 76 against the reduced frequency and freestream AoA amplitude.
Both CL and CD amplitudes seem to increase linearly with the freestream AoA amplitude or  f ,
meaning that if AoA amplitude is doubled while keeping the reduced frequency the same, CL and
CD will also double in amplitudes. CL and CD amplitudes seem to be almost linear to reduced

83

frequency too when the reduced frequency is higher than π/2. At low reduced frequencies, CL
and CD amplitudes change much less with the reduced frequency. Understandably, the flow
varies more slowly and becomes quasi-steady as the reduced frequency gets close to zero. In a
quasi-steady flow, enough time is given to the flow to adjust to changes in AoA, and CL and CD
values become close to those of the corresponding steady freestream flow.

,  f =4

14

5

CL,  f =4, k=2

12

4

CL,  f =4, k=

10

3

CL,  f =4, k=/2

8

2

CL,  f =4, k=/4

6

CL



CL,  f =4 , k=/8

1

4

0

2

-1

0

-2

-2

-3

-4

-4

0

45

90

135

180
225
270
315
t (deg.)
Figure 70: Time variations of CL for flows with oscillating AoA, cases A1-5.

84

-6
360

 (deg.)

6

0.5

,  f =4

0.4

CD,  f =4, k=2

14
12



CD,  f =4, k=/2

0.2

CD,  f =4, k=/4

0.1

CD,  f =4, k=/8

10
8
6

 (deg.)

CD,  f =4 , k=

0.3

CD

4
0
2

-0.1

0

-0.2

-2

-0.3

-4

-0.4

0

45

90

135

180
225
270
315
t (deg.)
Figure 71: Time variations of CD for flows with oscillating AoA, cases A1-5.
,  f =8

6

-6
360

14



CL,  f =8 , k=

4



CL,  f =8 , k=/2

3

CL,  f =8, k=/4

CL

2

CL,  f =8, k=/8

12
10
8
6

1

4

0

2

-1

0

-2

-2

-3

-4

-4

0

45

90

135

180
225
270
315
t (deg.)
Figure 72: Time variations of CL for flows with oscillating AoA, cases A6-10.

85

-6
360

 (deg.)

CL,  f =8 , k=2

5

,  f =8

0.5

CD,  f =8, k=2

14

CD,  f =8, k=

12

0.3

CD,  f =8, k=/2

10

0.2

CD,  f =8, k=/4

8

CD,  f =8, k=/8

6

CD

0.1

4
0

 (deg.)

0.4

2
-0.1

0

-0.2

-2

-0.3

-4
0

45

90

135

180
225
270
315
t (deg.)
Figure 73: Time variations of CD for flows with oscillating AoA, cases A6-10.

-6
360

,  f =8

*

CL

10

C* ,  f =8, k=/2
L

4

12

C* ,  f =8, k=
L

6

14

C* ,  f =8, k=2
L

8

8

C* ,
L
C* ,
L

2



 f =8 , k=/4


 f =8 , k=/8

6
4
2

0

0
-2

-2

-4
-4

0

45

90

135

180
225
270
315
t (deg.)
*
Figure 74: Time variations of CL for flows with oscillating AoA, cases A6-10.

86

-6
360

 (deg.)

-0.4

C* ,  f =8, k=2
D

8

C* ,  f =8, k=/4
D

0.2

10

C* ,  f =8, k=/2
D

0.4

12

C* ,  f =8, k=
D

0.6

14

6

C*

D

C* ,
D



 f =8 , k=/8

0

4
2

 (deg.)

,  f =8

0.8

0

-0.2

-2
-0.4
-0.6

-4
0

45

90

135

180
225
270
315
t (deg.)
*
Figure 75: Time variations of CD for flows with oscillating AoA, cases A6-10.

-6
360

5
CL,  f =4

4.5

CD,  f =4

4

CL,  f =8

Amplitude

3.5

CD,  f =8

3

CL,  f =4 (Theodorsen)

2.5

CL,  f =8 (Theodorsen)

2
1.5
1
0.5
0

0

1

2

3

4

5

6

7

k

Figure 76: Amplitude of oscillations in CL and CD for flows with oscillating freestream AoA and
the prediction by Theodorsen’s theory for a pitching airfoil in the steady mean freestream.

87

A notable difference between the results presented in Section 5.1.1 with those shown in this
section is magnitude of CL and CD amplitudes. For flows with oscillating freestream velocity
magnitude, the oscillation amplitudes in CL and CD are comparable to each other. However, for
flows with oscillating freestream AoA, the CL amplitude is one order of magnitude larger than
the CD amplitude. As explained in Section 5.1.1, when the freestream has constant AoA and
oscillating velocity magnitude, the acceleration and deceleration in freestream velocity induce
pressure gradients in the flow direction. However, when the AoA oscillates, the direction of
pressure gradient vector will also change in time with change in AoA. Figure 77 shows contours
of pressure coefficient around the airfoil at the time of maximum acceleration of freestream AoA.
As expected, the direction of pressure gradient in freestream, induced by changes in velocity
vector is not parallel to the mean flow direction and has a significant component perpendicular to
the mean flow direction, and therefore the pressure difference helps to significantly increase the
lift. The pressure coefficient on the airfoil surface is shown in Figure 78. The large magnitudes
of CP in Figures 77 and 78 are mainly due to large pressure gradient in the freestream flow.

Figure 77: Contours of pressure coefficient around the airfoil at the time of maximum
acceleration of freestream AoA for case A1 with k  2 and  f  4 .
88

74
73
Lower surface

Cp

72

Upper surface

71
70
69

0

0.2

0.4

0.6

0.8

1

X/ C
Figure 78: Pressure coefficient on the airfoil surface at the time of maximum acceleration of

freestream AoA for case A1 with k  2 and  f  4 .

Figure 79 shows that the lift and drag forces lead the freestream AoA in phase. Similar to
the phase shift described in Section 5.1.1, the phase shift here is calculated as the phase
difference between the peak of CL (and CD) and the peak of freestream AoA. Similar to flows
with oscillating freestream velocity magnitude, the phase shift is almost not affected by the
amplitude of oscillations in AoA. However, the CD phase shift increases as the reduced
frequency decreases, which is opposite to the flows with oscillating freestream velocity
magnitude. The minor decrease in CD phase shift at low reduced frequencies (compare case A5
to case A4 for example) is due to the CD profiles being distorted from the “sinusoidal” shape. At
low reduced frequency, the CD profile is “flattened” at high values and the peak values are
delayed, and therefore the phase shift is reduced. Similar to flows with oscillating freestream
velocity magnitude, the reduced frequency affects CL phase shift more than CD phase shift. Phase
shifts in CL and CD seem to converge to 90 degrees at high reduced frequencies, indicating that
the maximum values of the aerodynamic forces coincide with the maximum acceleration of the

89

freestream AoA. This is due to the fact that, at high values of reduced frequency, the noncirculatory or apparent mass terms dominate the solution [54].
Like Figure 76, Figures 79 and 80 compare the lift and drag coefficients of an airfoil
pitching around the chord mid-point given by Theodorsen’s theory [26] with those obtained by
LES here for the fixed SD7003 airfoil in flows with oscillatory direction or AoA. Gharali and
Johnson [36] compared CFD results for flows with oscillating freestream direction around a
S809 airfoil to experimental results of a pitching airfoil, and found some level of agreement
between them. Since the computational domain used in this study has a far-field boundary that is
close to a perfect circle centered at the chord mid-point, it is speculated that flow configuration
of these simulations might be equivalent, in terms of force acting on the airfoil, to the flow with
steady freestream velocity U m and the airfoil pitching around its chord mid-point with mean
AoA  m  4 and AoA amplitude  f . For a pitching airfoil, Theodorsen’s theory gives the lift
coefficient in Equation (43) where F and G are Theodorsen’s function as described in Section
5.1.1. For CL, Figures 76 and 79 show that the predictions of Theodorsen’s theory for the phase
shift and amplitude agree with the LES results quite well with only a slight discrepancy in CL
amplitude at high values of reduced frequency, even though some of the assumptions made in the
theory are not necessarily satisfied for the simulated flow. However, Figure 80 shows that there
is a significant discrepancy in mean CL values obtained by LES for flows with oscillatory
freestream AoA with the predictions of Theodorsen’s theory for a pitching airfoil in a steady
mean freestream.



CL   f k cos t  k 2 sin t

 2  F m   f




kG 

F 
 sin t   f
2 


90


 kF

 G  cos t 

 2



(43)

140
120

Phase Shift (deg.)

100
80
CL,  f =4
CD,  f =4

60

CL,  f =8

40

CD,  f =8
CL,  f =4 (Theodorsen)

20
0

CL,  f =8 (Theodorsen)

0

1

2

3

4

5

6

7

k

Figure 79: Phase shift of CL and CD for flows with oscillating AoA and the prediction by
Theodorsen’s theory for a pitching airfoil in the steady mean freestream. The two Theodorsen
curves overlay each other.
0.7

Mean CL

0.6
 =4, LES

0.5

f



 =8 , LES
f

0.4



 =4 , theory Theodorsen
f

 f=8, theory Theodorsen

0.3

0.2

0

1

2

3

4

5

6

7

k

Figure 80: Computed time-averaged lift coefficient for flows with oscillating freestream AoA
and the prediction by Theodorsen’s theory for a pitching airfoil in the steady mean freestream.
The two Theodorsen curves overlay each other.

91

For flows with steady freestream, both CL and CD increase monotonically with AoA until
AoA reaches the stall angle of about 11°for Re=60,000 [20]. The maximum AoA considered in
this section is 12° slightly above the stall angle. Therefore, when the reduced frequency is very
,
low and the flow has a quasi-steady freestream condition, the maximum CL and CD should
coincide with the maximum AoA, i.e. the phase shift between the aerodynamic force and the
freestream AoA should be nearly zero. The profiles of CL phase shift in Figure 79 clearly show
that the lift coefficient is nearly in phase with the AoA at low frequencies. However, the profiles
of CD phase shift seem to follow a different trend and still show phase shifts higher than 100
degrees at low frequencies. This indicates that the minimum reduced frequency considered in
this study might still be high enough to induce significant phase shift in the CD.
5.2.2 Vorticity field
A comparison of vorticity iso-surfaces/contours made between flows with steady freestream
flows at Re=60,000 and AoA=4° 8° 12°(Figures 81 and 82) indicates typical effects of AoA on
, ,
vorticity field. For example, as the AoA increases, the separation point moves towards the
leading edge and the recirculation region grows in size. However, the trends are very different
when the AoA oscillates. This is illustrated in Figures 83-94, where the z-vorticity isosurfaces/contours, and the line integrals of spanwise-averaged vorticity magnitude

  dn ,

   dnds for cases A1, A5, A6, and A10 are shown. The z-vorticity iso-surfaces and contours
show that lower reduced frequency and higher amplitude in the AoA oscillation lead to wider
variations in the vorticity field in terms of vortex size and the starting point of the recirculation
region. This is consistent with

  dn and    dnds plots in Figures 85, 88, 91, and 94. For

example, in case A6 with k=2π the

  dn profiles at four different times as shown in Figure 91
92

are very similar and the difference between the maximum and the minimum
about 0.2, while in case A10 with k=π/8 the

   dnds values is

  dn profile varies significantly over time and the

difference between the maximum and the minimum

   dnds

values is about 0.7. A

comparison between a flow with oscillating freestream AoA and a steady freestream flow with
the same maximum AoA shows that the unsteady flow has a smaller recirculation region. For
example, in the extreme case A10 with k=π/8 and  f  8 (Figures 92 and 93) the freestream
AoA varies from -4° to 12° during one freestream cycle, but the maximum size of the
recirculation region is similar to that of the steady freestream flow with AoA=8° which is much
,
smaller than that of the steady freestream flow with AoA=12°(Figures 81 and 82).
By relating the CD values of case A10 (  f  8 , k=π/8) in Figure 73 with the z-vorticity
iso-surfaces in Figure 92, it can be seen that the boundary layer separation occurs around the
same time (t0+T/2) when the CD value starts to flatten out. Compared to flows oscillating with
high reduced frequency, those with low reduced frequency have more time to develop and thus
cannot sustain the adverse pressure gradient and experience separation near the leading edge.
After the separation, the static pressure in the recirculation region recovers from its low values,
keeping the drag coefficient low and the CD profile flattened.

93

AoA=4°

u/Um
2.00
1.75

AoA=8°

1.38
1.00

AoA=12°

0.62
0.25
-0.12
-0.50
Figure 81: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for flows with steady
freestream, Re=60,000 and AoA=4° 8°and 12°
,
.

z C U m

AoA=4°

-12
-16

AoA=8°

-20
-24
-28
-32

AoA=12°

-36
-40
Figure 82: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for flows with
steady freestream, Re=60,000 and AoA=4° 8°and 12°
,
.

94

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 83: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A1 with k  2
and  f  4 . At t0, α(t) = 4° at t0+T/4, α(t) = 8° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) = 0°
;
;
;
.

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4  (t )

-40

Figure 84: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A1 with

k  2 and  f  4 . At t0, α(t) = 4° at t0+T/4, α(t) = 8° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) =
;
;
;
0°
.

95

4

2.5
(b)

||dnds

3

||dn

(a)

t0
t0+T/4
t0+T/2
t0+3T/4

2

2

1.5

1

0
0

0.2

0.4

0.6

0.8

1
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 85: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A1 with k  2 and

 f  4 . (a) Integrated along n axis; (b) integrated along n and then s axes.

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 86: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A5 with k   8
and  f  4 . At t0, α(t) = 4° at t0+T/4, α(t) = 8° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) = 0°
;
;
;
.

96

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

 (t )

Figure 87: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A5 with

k   8 and  f  4 . At t0, α(t) = 4° at t0+T/4, α(t) = 8° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) =
;
;
;
0°
.

4
(a)

(b)

||dnds

||dn

3

2.5
t0
t0+T/4
t0+T/2
t0+3T/4

2

2

1.5

1

0
0

0.2

0.4

0.6

0.8

1
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 88: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A5 with k   8 and

 f  4 . (a) Integrated along n axis; (b) integrated along n and then s axes.

97

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 89: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A6 with k  2
and  f  8 . At t0, α(t) = 4° at t0+T/4, α(t) = 12° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) = -4°
;
;
;
.

t0

z C U m
-12
t0+T/4

-16
-20
-24

t0+T/2

-28
-32

t0+3T/4

-36

 (t )

-40

Figure 90: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A6 with

k  2 and  f  8 . At t0, α(t) = 4° at t0+T/4, α(t) = 12° at t0+T/2, α(t) = 4° at t0+3T/4, α(t)
;
;
;
= -4°
.

98

30

||dn

20

(b)

(a)

t0
t0+T/4
t0+T/2
t0+3T/4

||dnds

25

5

15
10

4.5

4

5
0
0

0.2

0.4

0.6

0.8

3.5
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 91: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A6 with k  2 and

 f  8 . (a) Integrated along n axis; (b) integrated along n and then s axes.

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 92: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A10 with k   8
and  f  8 . At t0, α(t) = 4° at t0+T/4, α(t) = 12° at t0+T/2, α(t) = 4° at t0+3T/4, α(t) = -4°
;
;
;
.

99

z C U m

t0

-12
-16
t0+T/4

-20
-24
-28

t0+T/2

-32
-36
t0+3T/4

-40

 (t )

Figure 93: Instantaneous contours of spanwise-averaged z-vorticity ( z C U m ) for case A10
with k   8 and  f  8 . At t0, α(t) = 4° at t0+T/4, α(t) = 12° at t0+T/2, α(t) = 4° at t0+3T/4,
;
;
;
α(t) = -4°
.
4

2.5

||dn

3

(a)

(b)

||dnds

t0
t0+T/4
t0+T/2
t0+3T/4

2

2

1.5

1

0
0

0.2

0.4

0.6

0.8

1
t0

1

X/C

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 94: Integrals of spanwise-averaged vorticity magnitude |Ω| for case A10 with k   8 and

 f  8 . (a) Integrated along n axis; (b) integrated along n and then s axes.
5.2.3 Velocity fluctuations
The effect of freestream AoA oscillations on velocity fluctuations over the airfoil is studied
in this section by considering the RMS values of various velocity components and TKE for case
A6. Case A6 is selected since it has the highest frequency of k=2π and the largest amplitude of

 f  8 . Figure 95 shows that the integrated values of velocity RMS along airfoil surface

100

normal direction, (n),

  RMS dn at four different times during a cycle of oscillation in AoA are

nearly the same. This is further confirmed by the integrated values of velocity RMS along airfoil
surface normal and tangent directions, (n, s),

   RMS dnds in Figure 96 and those of TKE in

Figure 97. It is evident that the velocity fluctuations at different phases, as measured by RMS
values are similar when the oscillation in freestream direction has a high reduced frequency. It
can also be seen in Figures 95 and 96 that the RMS of x-velocity and RMS of y-velocity
components (u′ and v′) are comparable and are one order of magnitude higher than the RMS of
z-velocity component (w′). The RMS of y-velocity is slightly higher than the RMS of x-velocity
(Figure 96), which shows that the oscillation in freestream AoA can cause more fluctuation in yvelocity than in x-velocity.
The trends in Figures 96 and 97 are similar to those in Figures 63 and 64, indicating that
qualitatively the effects of freestream oscillations in AoA on the airfoil flow velocity fluctuations
are similar to those caused by oscillations in the freestream velocity magnitude. However, the
velocity fluctuations induced by AoA oscillations are observed to be quantitatively much more
significant than those caused by velocity magnitude oscillations.

101

(RMS)dn

1.2
0.9
0.6

u'/Um

t0

v'/Um
w'/Um

0.3

(RMS)dn

0
1.2
t0+T/4

0.9
0.6
0.3

(RMS)dn

0
1.2
t0+T/2

0.9
0.6
0.3

(RMS)dn

0
1.2
t0+3T/4

0.9
0.6
0.3
0
0

0.2

0.4

0.6

0.8

1

X/C
Figure 95: Integral of velocity RMS along airfoil surface normal direction, n for case A6 with

k  2 and  f  8 .

102

(RMS)dnds

0.3

0.2
u'/Um
v'/Um

0.1

w'/Um

0
t0

t0+T/4

t0+T/2
Time

t0+3T/4

t0+T

Figure 96: Integral of velocity RMS along n and s axes for case A6 with k  2 and  f  8 .

0.012

(TKE)dnds

0.01
0.008
0.006
0.004
0.002
0
t0

t0+T/4

t0+T/2
t0+3T/4
t0+T
Time
Figure 97: Integral of turbulence kinetic energy (TKE) along n and s axes for case A6 with

k  2 and  f  8 .

103

5.2.4 Boundary layer separation and reattachment
In this section, the effects of oscillating freestream direction or AoA on the boundary layer
separation and reattachment are examined. Mean and instantaneous values of the skin friction
coefficient on the airfoil for cases A1, A3, A5, A6, and A10 are plotted in Figures 98-102.
Similar plots for the steady freestream flow with Re=60,000 and AoA=4°were previously shown
in Figure 66 in Section 5.1.4. A comparison between steady and unsteady cases indicates the C F
variation, or profile over the airfoil becomes very different when the freestream flow direction
oscillates. The CF profiles are also very different at different phases or times. Figures 98-100
show that as the reduced frequency of oscillation in AoA increases, the separation point moves
less. This is in agreement with the observation made in the previous section that higher
frequency of oscillation in freestream AoA leads to lower variation in the vorticity magnitude.
Another observation made from CF plots in Figures 98-100 is that, close to the trailing edge, CF
oscillates around the mean value with larger amplitude as the reduced frequency increases. This
is due to vortex shedding from the trailing edge that occurs at high reduced frequencies (Figures
83, 84, 89 and 90). Because of oscillation in the freestream flow direction, the vortex shed from
the airfoil upper surface near the trailing edge can move either upward or downward. For the
vortex to be shed downward, the fluid above the airfoil upper surface and near the trailing edge
will have to impinge the airfoil surface harder in comparison to the flow with the vortex shed
upward, thus the instantaneous CF has to be higher than the averaged CF. A comparison between
Figure 100 (case A1) and Figure 101 (case A6) shows that as the amplitude of oscillation in AoA
increases, the instantaneous separation point moves upstream towards the leading edge, and the
near trailing edge values of CF also oscillate more. A comparison between Figure 101 (case A6)
and Figure 102 (case A10) shows that the oscillation in CF around the trailing edge subsides

104

significantly when the reduced frequency drops from 2π to π/8, which again shows the
significance of vortex shedding at the trailing edge and the effect it has on the CF.
The mean separation and reattachment points reported in Table 8 are obtained based on the
mean skin friction coefficient. No separation (and thus reattachment) is found for cases A1 and
A6-10. For cases A2-5 which have separation and reattachment points, the mean separation point
is delayed by an increase in the reduced frequency. Compared to the flow with steady freestream
velocity, results of the ten cases considered here show that oscillations in freestream flow AoA
generally suppress the mean separation.

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 98: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for case A5 with k   8 and  f  4 .

105

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 99: Time-averaged and instantaneous skin friction coefficient on the airfoil upper surface
for case A3 with k   2 and  f  4 .

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 100: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A1 with k  2 and  f  4 .

106

0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 101: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A6 with k  2 and  f  8 .
0.03
Mean
t0
t0+T/4
t0+T/2
t0+3T/4

0.025
0.02
0.015

C

F

0.01
0.005
0
-0.005
-0.01
-0.015
-0.02

0

0.2

0.4

0.6

0.8

1

X/C

Figure 102: Time-averaged and instantaneous skin friction coefficient on the airfoil upper
surface for case A10 with k   8 and  f  8 .

107

5.3 A Discussion Regarding Effect of Freestream Unsteadiness on Mean Separation and
Reattachment Points
There have been some discussions in the literature regarding the effect of freestream flow
perturbations on the mean separation and reattachment points. Measurements of laminar
separation bubble (LSB) over SD7003 airfoil at Re=60,000 and AoA=4°were conducted at four
different facilities [9, 29] including the tow tank at Institute for Aerospace Research (IAR), the
low-noise wind tunnel at Technical University of Braunschweig (TU-BS), the free-surface water
tunnel at Air Force Research Laboratory (AFRL), and the closed-circuit water tunnel at RWTH
Aachen University. Mean separation and reattachment data from these facilities are listed in
Table 9, along with the estimated freestream turbulence intensity in the facilities. Evidently, IAR
and TU-BS data agree quite well, the slight difference can be partly due to the finite aspect ratio
of the wing in the IAR facility which induces a smaller effective AoA than those in the other
three cases [9]. In the AFRL experiment, the separation bubble forms considerably further
upstream and the boundary layer reattaches also further upstream in comparison to IAR and TUBS data, which is likely caused by a true AoA that may be slightly larger than the nominal [9].
The RWTH data show delayed separation and advanced reattachment (and consequently smaller
LSB) compared to the other three data sets. This seems to be caused mainly by the substantially
higher freestream turbulence intensity in the RWTH facility.

Data Set
IAR
TU-BS
AFRL
RWTH

Freestream Turbulence
Intensity, (%)
0
0.1
~0.1
1

Separation X s C

Reattachment X r C

0.33
0.30
0.18
0.390

0.63
0.62
0.58
0.515

Table 9: Measured separation and reattachment points for AoA=4°and Re=60,000 from different
experimental facilities.

108

Burgmann et al. [29, 52] at RWTH Aachen University measured the flow over SD7003
airfoil with Re=20,000 and AoA=4° at two levels of freestream turbulence. The freestream
turbulence intensity was 1.5% in their earlier experiments, but with additional flow conditioning
devices in their later experiments it was lowered to 1%. The measured separation and
reattachment points from these experiments are reported in Figure 103, which clearly indicates
that the separation point moves downstream at higher turbulence intensity, except in the case
with an AoA of 8° At this high AoA, the strong adverse pressure gradient at the separation
.
region could make the effect of freestream turbulence intensity become insignificant. These data
do not support the discussion about effect of turbulence intensity on reattachment, because the
method of determining the reattachment point in the later experiments is fundamentally different
from that in the earlier experiments. In the earlier experiments, both separation and reattachment
points were determined using time-mean streamlines. In the later experiments, an instantaneous
reattachment point was determined at each time step dividing the flow field into a main
recirculation region and a region of shed vortices, and then the variables were averaged to be the
time-mean reattachment point. It was noted by Burgmann et al. that this new method results in
reattachment points to significantly move upstream in comparison with points obtained based on
mean streamlines.

109

1
separation, Tu=1.5%
separation, Tu=1%
reattachment, Tu=1.5%
reattachment, Tu=1%

0.9
0.8

Location (x/c)

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

3

4

5

6
7
Angle of attack ()

8

9

Figure 103: Comparison of the separation and reattachment points for freestream turbulence
intensity of 1.0% and 1.5% at Re=20,000 [29, 52].

Galbraith [20] used LES to investigate the sensitivity of the separation bubble size and
shape to the leading edge disturbances in steady flows. The disturbance was selected to be a wall
normal zero-net mass-flow blowing/suction slot with a streamwise sinusoidal distribution, and
the slot was positioned near the leading edge. Two disturbance velocity amplitudes, with 0.1%
and 1% of the freestream velocity, were considered. It was found that, for both disturbance
velocity amplitudes, the separation point moves downstream while the reattachment point moves
upstream, resulting in diminished separation bubbles compared to the case without disturbances.
Olson [12] measured the LSB size over a stationary SD7003 airfoil under different levels of
freestream turbulence and verified that the increased freestream turbulence moves the separation
point downstream and the reattachment point upstream. It was further shown that the increase in

110

freestream turbulence can even eliminate the reversed flow region in the mean flow plots (within
the resolution of the measurements).
In the current study, the mean separation and reattachment results in Sections 5.1.4 and
5.2.4 show that the harmonic oscillations of freestream velocity, either in magnitude or direction,
will delay the mean separation, pushing it downstream and trigger earlier mean reattachment
closer to the leading edge up to the point that the mean flow may become fully attached.
Compared to the studies mentioned above, flows in this study have much higher level of
freestream fluctuations. This implies that suppression of LSB could occur under much higher
level of freestream fluctuations than freestream turbulence. Nevertheless, our findings are
consistent with most of experimental data reported in the literature.

111

5.4

Flow with Wind Gust
This section considers the aerodynamic response of the SD7003 airfoil to a wind gust. Wind

gust is difficult to predict or model as it depends on many factors such as the geographical
location and season. This is evident in the measured time variation of wind speed shown in
Figure 104 [56, 57]. As a universal model of wind gust is not realistic, and a simple model can be
very helpful in better understanding of flow over airfoils in harsh wind gust environment.

Wind velocity (m/s)

14
12
10
8
0

10

20

30
Time (min)

600

800

40

50

60

Wind velocity (m/s)

40
30
20
10
0
0

200

400

1000

1200

1400

1600

1800

Time (sec)
Figure 104: Measured wind speed at (upper) a sea mast in Vindeby offshore wind farm [56] and
(lower) a 10 meter tower [57].

Here, the flow over SD7003 airfoil is computed by LES with a wind gust model [58]
incorporated into the freestream flow. The model involves a sudden increase and smoother
112

decrease in wind speed with added oscillations, while the wind direction is kept constant with a
fixed AoA of 4° It takes the following form
.

U  t   U 0  G1  G2  G3 .

(44)

where U 0 is the freestream velocity prior to the wind gust, and G1 , G2 , G3 are defined as

G1  D1t exp( D2t ) ,

(45)

G2  D3erf  D4t  ,

(46)

G3  D5 exp  D6t  sin  D7t  ,

(47)

D1 ~ D7 are the model constants. The “damping” function, G1 causes a sudden increase in the
freestream wind speed, the saturation function G2 sets the speed after the wind gust, and G3
controls the perturbations of the wind gust. The model constants can be adjusted to reflect
different type of “wind gusts”.
The model constants in Table 10 generate the wind gust shown in Figure 105. The Reynolds
number for the steady flow before the wind gust is Re=60,000. The sudden increase in velocity
drives the Reynolds number up to maximum Re=74,000. It takes about 3 residence times
(residence time= C U 0 ) for the freestream velocity to decline from the maximum velocity and
stabilize at Re=63,000.

Constants

D1

D2

D3

D4

D5

D6

D7

Values

700U 0

7 (3TR )

1

1 TR

0.15

7 (10TR )

20 TR

Table 10: Constants of the wind gust model.

113

The time variation of lift and drag coefficients obtained by LES is shown in Figure 105. The
lift and drag coefficients are defined in the same way they were defined in the previous sections,
with the freestream velocity prior to the wind gust U0 as the reference velocity in both
coefficients. Prior to the wind gust, the flow with steady freestream is fully developed with CL
and CD fluctuating slightly around their time-averaged values of 0.59 and 0.020. As the wind
gust moves over the airfoil, CL and CD rapidly increase to maximum values of 0.84 and 0.17,
respectively, and then start to decrease as the freestream velocity declines. The lift coefficient
closely follows the trend in the freestream velocity and stabilizes at a value slightly higher than
that prior to the wind gust since the lift force is normalized by U0 which is lower than U after the
wind gust as shown in Figure 105, while the drag coefficient goes negative before it stabilizes at
a value very close to that prior to the wind gust. The perturbation parameter G3 generates
oscillations in CL and CD with the frequency comparable to that of the freestream flow. The
amplitudes of fluctuations in CL and CD diminish as the amplitude of wind perturbation decays.
The maximum amplitude in G3 is less than 1% of U0, yet it causes fluctuations in CL and CD with
amplitude up to 3% of the maximum CL value and 15% of the maximum CD value. This suggests
that even a weak (low amplitude) wind gust can have a significant effect on the aerodynamic
forces and the flight stability of MAVs. The lift and drag coefficients obtained by normalizing
the forces with the time-dependent U(t) are also plotted in Figure 105, and as expected, are lower
than those obtained by normalizing the forces with U0.
The iso-surfaces of z-vorticity of the flow over the airfoil at various times are shown in
Figure 106. As the freestream velocity suddenly increases due to the wind gust, the separation
point is pushed downstream toward the trailing edge. But when the wind gust starts to diminish,
the separation point begins to move back and even slightly passes its original (pre-wind gust)

114

position. As the freestream velocity continues to stabilize at a velocity slightly higher than the
original velocity, the separation point gets close to steady state flow value.

115

G1/U0

0.2
0.1
0
G2/U0

0.04
0.02

G /U

0

0
0.01

3

0
-0.01
3

U/U0

1.2

2

4
5

1.15

6

1.1

7 8 9 10

1.05

12

11

13

1

1

CL

0.8

normalized by U

0.7

normalized by U(t)

0

0.6
0.5
0.4

C

D

0.15

normalized by U0

0.1

normalized by U(t)

0.05
0
-0.05

0

1

2

3

4

5
6
tU0/C

7

8

9

10

11

Figure 105: Time variations of the wind gust model parameters, velocity and lift and drag
coefficients induced by the wind gust. CL and CD are computed by normalizing the lift and drag
forces with U0 and U(t), both.

116

1

8

2

9

3

10

4

11

5

12

6

13

7

Figure 106: Instantaneous iso-surfaces of z-vorticity ( z C U0  10 ) simulated by LES with
wind gust at different times. The numbers on figures match the times in Figure 105.

117

6.

Summary and Conclusions

Various methods are used to study flows around the SD7003 airfoil at different flow
conditions. First, linear stability methods are used to investigate the local and global instability of
the time- and spanwise-averaged two-dimensional base flow. Second, large eddy simulation
(LES) method is used to simulate steady and unsteady freestream flows around the airfoil.
For the stability analysis, the LSB region of mean (time- and spanwise-averaged) flow field
around the SD7003 airfoil at Reynolds number of 60,000 and AoA of 4°is employed. The base
flow is obtained using a compressible model and an implicit LES solve. The eigenvalue system
is solved in the global stability analysis and the classic Rayleigh equation is solved in the local
stability analysis. The maximum growth rates obtained from the two methods are in agreement,
and both methods correctly predict that the flow is unstable and will transition to turbulence.
To validate the LES method, a number of flows over SD7003 airfoil with steady freestream
are simulated, and the results are compared to available numerical and experimental data. The
LES results are found to compare reasonably well with other LES data in the literature. The LES
results are also consistent with the experimental data, despite discrepancies among the
experimental data obtained from different facilities.
The effect of flow three-dimensionality is studied by comparing two-dimensional (2D) and
3D solutions of flow over SD7003 airfoil. The 2D and 3D solutions are found to be generally
different. While the difference in instantaneous velocity, pressure and vorticity fields are quite
significant, the global quantities like the lift or drag coefficients are not very different unless in
extreme situations in unsteady flows. The 2D solutions fail to capture the 2D to 3D transition due
to lack of vortex stretching term and spanwise flow variables.

118

The main objective of this thesis is to study the aerodynamic response of the SD7003 airfoil
to unsteady freestream flows. Three types of unsteady freestream conditions are considered. In
the first type, AoA is fixed at 4° while the freestream velocity magnitude varies periodically
,
with mean Reynolds number of 60,000, reduced frequencies ranging from π/8 to 2π, and
normalized oscillation amplitudes of 0.183 and 0.366. In the second type, the freestream velocity
magnitude is fixed at Reynolds number of 60,000, while the freestream flow AoA varies
periodically around the mean AoA value of 4°with the reduced frequency ranging from π/8 to
2π and the oscillation amplitudes of 4°and 8° In the third type of unsteady flows considered in
.
this thesis, a wind gust model is applied to the freestream flow, causing sudden increase and
subsequent gradual decrease in velocity magnitude.
For flows with oscillating freestream velocity magnitude, the oscillations of aerodynamic
forces, vorticity magnitude, separation and reattachment points are studied in details. Our results
indicate that the amplitudes of lift and drag coefficients increase as the reduced frequency and
amplitude of freestream velocity magnitude increase. However, the freestream velocity
oscillation does not noticeably change the mean lift and drag coefficients from their steady mean
freestream values. There is a phase shift between the aerodynamic forces and the freestream
oscillation, which is mainly a function of reduced frequency. Compared to the inviscid theory,
the lift coefficient obtained from the LES simulations have significantly higher amplitudes, while
phase shifts are predicted similarly by the theory and simulation. The oscillating freestream
velocity magnitude causes the vorticity field to change slightly during a cycle. For the cases
studied, the changes in vorticity field are not significant enough to show a clear trend. However,
the integral of the vorticity magnitude along the airfoil wall normal and tangent, n and s
(    dnds ) shows strong correlation with the freestream velocity magnitude. The root-mean-

119

square (RMS) values of flow, calculated for the case with the highest reduced frequency and the
largest amplitude indicate that the fluctuations in various velocity components represented by
velocity RMS values are quite the same at different phases of a cycle. The intensity of
fluctuations in x-velocity (u′) is slightly higher than that of y-velocity (v′), while the fluctuation
of z-velocity (w′) has a much lower magnitude compared to the other two velocity components.
Our results also indicate that the separation and reattachment points oscillate over a wider range
as the reduced frequency and freestream amplitude increase.
Compared to the flows with oscillating freestream velocity magnitude, the flows with
oscillating freestream direction have some similar features but are also different in some respects.
The amplitudes of lift and drag coefficient increase as the reduced frequency and amplitude of
oscillations in flow AoA increase. However, the lift amplitude is one order of magnitude larger
than the drag amplitude. The phase shift between aerodynamic forces and freestream AoA is
significantly affected by the reduced frequency. Compared to the flow with steady freestream
condition, the mean lift is almost unchanged; but, the mean drag is reduced due to the Katzmayr
effect. It is shown that lower reduced frequency and higher amplitude of oscillations in AoA lead
to bigger reduction in drag. The amplitude and phase shift of the lift coefficient compare quite
well with the inviscid theory developed for a pitching airfoil in a steady freestream, while there
is significant difference in the mean lift. At low reduced frequencies and high AoA amplitudes,
the vorticity field changes dramatically during a cycle, and this is related to the wide variations
in the separation point. At high reduced frequencies, the vorticity field shows less variation in
time, but strong vortex shedding at the trailing edge. The strong vortex shedding causes a wide
variation in the skin friction coefficient near the trailing edge. For the case with the highest

120

reduced frequency and the largest amplitude of oscillations in flow AoA, the magnitude of
velocity fluctuations is nearly unchanged in different phases.
Based on the mean separation and reattachment points calculated for unsteady flows with
freestream velocity magnitude and direction (or AoA) oscillations, it is shown that the
oscillations in freestream velocity, either in magnitude or direction, delay the mean separation
and promote earlier reattachment of the boundary layer up to the point that the mean flow
becomes fully attached.
Simulations conducted with the wind gust model shows sudden increase in the aerodynamic
load at the beginning of the gust. Generally, the aerodynamic forces behave the same way as the
freestream velocity but they show relatively high sensitivity to the small changes in freestream
flow.
The current work can be expanded in several ways if substantially better computational
resources become available. Stability analysis can be performed on finer mesh for more detailed
analysis of the onset of instability over a bigger region to develop a broader and better picture of
the flow instability. For flows with unsteady freestream, finer meshes should be used to obtain
details of the flow in regions away from the airfoil and in the wake particularly for flows at high
AoA. Much longer simulations are also needed for capturing the response of flow to very low
frequency oscillations in the freestream. These details may provide a clearer picture of the
development of boundary layer and wake vortices under the effect of oscillating freestream.

121

APPENDICES

122

APPENDIX A
The perturbation equations (13)-(16) are solved with the finite element method. For this, the
velocity and pressure eigenfunctions are interpolated as
ˆ
ˆ
v   v
ˆ
ˆ
p   p

  1, 2, , n 

ˆ
ˆ
where v and p are the node values,  are the shape functions, and n is the number of grid
nodes in an element. The coefficient matrices in Equation (19) for a 2D base flow are obtained
via the following equations.
 Auu
A
A   vu
A
 wu

 Apu

Auv
Avv
Awv
Apv

Auw
Avw
Aww
Apw

Aup 
Avp 
,
Awp 

App 



1


Auu    Ac  Aux  R e1 ( S  k 2 M )  Dxx  
3



1


Auv    Auy  R e1 Dxy 
3


1
T
Auw   k R e1 Cx
3
Aup  M a 1CxT
1


Avu    Avx  R e1 Dyx 
3



1


Avv    Ac  Avy  R e1 ( S  k 2 M )  Dyy  
3



1
Avw   R e1 kCT
y
3
Avp  M a 1CT
y
1
Awu   R e1 kCx
3
1
Awv   R e1 kC y
3

123

M
B




M
M




M



1


Aww    Ac  R e1 ( S  k 2 M )  k 2 M  
3




Awp  M a 1kM



Apv    Apy  M a 1C y 
Apu   Apx  M a 1Cx

Apw  kM a 1M

App   Ac

The elemental matrices are obtained from the following equations.
M      d 

 ,   1, 2, , n 





Ac     U   , x   V   , y d 

  1, 2, , n 

Aux       , xU d 
Auy       , yU d 

Avx       , xV d 
Avy       , yV d 

Apx       , x P d 

Apy       , y P d 


Cx    , x  d 
C y      , y d 

Dxx    , x   , x d 
Dxy    , x   , y d 

Dyx    , y   , x d 
Dyy    , y   , y d 

S  Dxx  Dyy

124

APPENDIX B

The plots of z-vorticity iso-surface not shown in Section 5.1.2 for flows with oscillating
freestream speed are presented here.

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 107: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U2 with k  
and   0.183 .

125

t0

u/Um
2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 108: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U3 with k   2
and   0.183 .
u/Um

t0

2.00
1.75
t0+T/4
1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 109: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U4 with k   4
and   0.183 .

126

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 110: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U7 with k  
and   0.366 .
u/Um

t0

2.00
1.75
t0+T/4
1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 111: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U8 with k   2
and   0.366 .

127

u/Um

t0

2.00
1.75
t0+T/4
1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 112: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case U9 with k   4
and   0.366 .

128

APPENDIX C
The plots of z-vorticity iso-surface not shown in Section 5.2.2 for flows with oscillating
freestream AoA are presented here.

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 113: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A2 with k  
and  f  4 .

129

u/Um

t0

2.00
1.75

t0+T/4

1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 114: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A3 with k   2
and  f  4 .
u/Um

t0

2.00
1.75
t0+T/4
1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 115: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A4 with k   4
and  f  4 .

130

u/U

t0

m

2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 116: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A7 with k  
and  f  8 .
u/U

t0

m

2.00
1.75

t0+T/4

1.38
t0+T/2
1.00
t0+3T/4

0.62
0.25
-0.12
-0.50

Figure 117: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A8 with k   2
and  f  8 .
131

u/Um

t0

2.00
1.75
t0+T/4
1.38
t0+T/2

1.00
0.62

t0+3T/4

0.25
-0.12
-0.50
Figure 118: Instantaneous iso-surfaces of z-vorticity ( z C U m  10 ) for case A9 with k   4
and  f  8 .

132

BIBLIOGRAPHY

133

BIBLIOGRAPHY
[1] Jones, B. M., “Stalling,” Journal of the Royal Aeronautical Society, Vol. 38, 1938, pp. 747–
770.
[2] Young, A. D., and Horton, H. P., “Some Results of Investigation of Separation Bubbles,”
AGARD Conference Proceedings, Vol. 4, AGARD, London, 1966, pp. 779–811.
[3] Tani, I., “Low Speed Flows Involving Bubble Separations,” Progress in Aeronautical
Sciences, Pergamon, New York, 1964, pp. 70–103.
[4] Tezuka, A., “Global stability Analysis of Attached or Separated Flows over a NACA.0012
Airfoil,” AIAA Paper 2006-1300.
[5] Hammond, D. A. and Redekopp, “Local and Global Instability Properties of Separation
Bubbles,” Eur. J. Mech. B/Fluids, Vol. 2, 1998, pp. 145-164.
[6] Theofilis, V., Hein, S. and Dallmann, U., “On the Origins of Unsteadiness and ThreeDimensionality in a Laminar Separation Bubble,” Phil. Trans. Roy. Soc. Lond. A. Vol. 358,
2000, pp. 3229-3246.
[7] Galbraith, M. and Visbal, M., “Implicit large Eddy Simulation of Low Reynolds Number
Flow Past the SD7003 Airfoil,” AIAA Paper 2008-225.
[8] Radespiel, R., Windte, J., and Scholz, U., “Numerical and Experimental Flow Analysis of
Moving Airfoils with Laminar Separation Bubbles,” AIAA Paper 2006-501, Jan. 2006.
[9] Ol, M. V., McAuliffe, B. R., Hanff, E. S., Scholz, U., and Kähler, C., “Comparison of
Laminar Separation Bubble Measurements on a Low Reynolds Number Airfoil in Three
Facilities,” AIAA Paper 2005-5149, Jun. 2005.
[10] Katz, A., Naguib, A., Koochesfahani, M., “Investigation of a low Reynolds number airfoil
using molecular tagging velocimetry”, American Physical Society, 62nd Annual Meeting of
the APS Division of Fluid Dynamics, November 22-24, 2009.
[11] Katz, A., “Molecular Tagging Velocimetry Measurements of the Low-Reynolds-Number
Flow Around an SD7003 Airfoil”, Master Thesis, Michigan State University, 2010
[12] Olson, D., “Facility and Flow Dependence Issues Influencing the Experimental
Characterization of a Laminar Separation Bubble at Low Reynolds Number,” Master Thesis,
Michigan State University, 2011.
[13] Olson, D., Katz, A., Naguib, A., Koochesfahani, M., Rizzetta, D., and Visbal, M., “An
Investigation of the Effect of Freestream Turbulence on the Laminar Separation Bubble on

134

an SD7003 Airfoil”, AIAA 2011-395, 49th AIAA Aerospace Sciences Meeting Including
The New Horizons Forum and Aerospace Exposition , 4-7 January 2011, Orlando, Florida..
[14] Radespiel, R., Graage, K., and Brodersen, O., “Transition Predictions Using ReynoldsAveraged Navier–Stokes and Linear Stability Analysis Methods,” AIAA Paper 91-1641,
June 1991.
[15] Stock, H. W., and Haase, W., “A Feasibility Study of eN Transition Prediction in Navier–
Stokes Methods for Airfoils,” AIAA Journal, Vol. 37, No. 10, 1999, pp. 1187–1196.
[16] Yuan, W., Khalid, M., Windte, J., Scholz, U., and Radespiel, R., “An Investigation of low
Reynolds-Number Flows Past Airfoils,” 23rd AIAA Applied Aerodynamics Conference,
Toronto, AIAA Paper 2005-4607.
[17] Lian Y., and Shyy W., “Laminar-turbulent Transition of a Low Reynolds Number Rigid or
Flexible Airfoil,” AIAA Journal, vol. 45, No. 7, 2007.
[18] Sherer, S. and Visbal, M., “Multi-resolution Implicit Large Eddy Simulations Using a HighOrder Overset-Grid Approach”, International Journal for Numerical Methods in Fluids, Vol.
55, pp. 455-482, 2007.
[19] Rizzetta, D., Visbal, M. and Morgan, P., “A High-Order Compact Finite-Difference Scheme
for Large-Eddy Simulation of Active Flow Control”, Progress in Aerospace Science, Vol.
44, pp. 397-426, 2008.
[20] Galbraith, M., “Implicit Large Eddy Simulation of Low-Reynolds-Number Transitional
Flow Past the SD7003 Airfoil”, Master Thesis, University of Cincinnati, 2009.
[21] Visbal, M., Gaitonde, D. and Roy, S., “Control of Transitional and Turbulent Flows Using
Plasma-Based Actuators”, AIAA Paper 2005-5302, 36th AIAA Plasmadynamics and Lasers
Conference, Toronto, Ontario, June 6-9, 2005.
[22] Bohl, D. and Koochesfahani, M., “MTV measurements of the Vortical Field in the Wake of
an Airfoil Oscillating at High Reduced Frequency”, Journal of Fluid Mechanics, Vol. 620,
pp. 63-68, 2009.
[23] Visbal, M., “High-Fidelity Simulation of Transitional Flows Past a Plunging Airfoil”, 47th
AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace
Exposition 5 - 8 January 2009, Orlando, Florida.
[24] McGowan, G., Gopalarathnam, A., OL, M., Edwards, J., and Fredberg, D., “Computation vs.
Experiment for High-Frequency Low-Reynolds-Number Airfoil in Pitch and Plunge,”
AIAA Paper 2008-0653, 2008.
[25] Koochesfahani, M., “Vortical Patterns in the Wake of an Oscillating Airfoil”, AIAA Journal,
Vol. 27, No. 9, 1989.

135

[26] Theodorsen, T., “General Theory of Aerodynamic Instability and the Mechanism of
Flutter,” N.A.C.A. TR No. 496, May 1934.
[27] von Ká n, T. and Sears, W. R., “Airfoil theory for non-uniform motion,” J. Aeronaut. Sci.
rmá
5, 1938, 379–390.
[28] “Unsteady Aerodynamics for Micro Air Vehicles,” RTO Technical report, TR-AVT-149,
2010.
[29] Burgmann, S., Dannemann, J., Schröder, W. “Time-resolved and volumetric PIV
measurements of a transitional separation bubble on an SD7003 airfoil”. Experiments in
Fluids, 44, 2007, pp. 609-622.
[30] Isaacs, R., “Airfoil Theory for Flows of Variable Velocity,” Journal of the Aeronautical
Sciences, Vol. 12, No. 1, 1945, pp. 113-117.
[31] Greenberg, J. M., “Airfoil in Sinusoidal Motion in a Pulsating Stream,” N.A.C.A. Report
No. 1326, June 1947.
[32] Gursul, I., Lin, H., and Ho, C., “Effects of Time Scales on Lift of Airfoils in an Unsteady
Stream”, AIAA Journal, Vol. 32, No. 4, April 1994.
[33] Katzmayr, R., “Effect of Periodic Changes of Angle of Attack on Behavior of Airfoils”,
N.A.C.A. Technical Memorandum No. 147, 1922.
[34] Toussaint, Kerneis, Girault, “Experimental investigation of the effect of an oscillating
airstream (Katzmayr effect) on the characteristics of airfoils,” N.A.C.A. Technical Report
No.202, 1924.
[35] Williams, D., Quach, V., Kerstens, W., Buntain, S., Tadmor, G., Rowley, C., Colonius, T.,
“Low-Reynolds Number Wing Response to an Oscillating Freestream with and without
Feed Forward Control,” AIAA 2009-143, 47th AIAA Aerospace Sciences Meeting
Including The New Horizons Forum and Aerospace Exposition , 5-8 January 2009, Orlando,
Florida.
[36] Gharali, K. and Johnson, DA., “Numerical modeling of an S809 Airfoil under Dynamic
Stall, Erosion and High Reduced Frequencies”, Appl Energy (2011),
doi:10.1016/j.apenergy.2011.04.037.
[37] Ding Y, and Kawahara M., “Linear Stability of Impressible Flow Using a Mixed Finite
Element Method,” J. Comp. Phys, Vol. 139, 1998, pp. 243-273.
[38] Malkus, D., “Eigenproblems associated with the discrete LBB condition for incompressible
finite elements,” International Journal of Engineering Science, Vol. 19, 10, 1981.

136

[39] Frank M. White, “Viscous Fluid Flow (Third Edition),” McGraw-Hill, NY, 2006.
[40] Huerre P., and Monkewitz P. A., “Absolute and Convective Instability in Free Shear
Layers,” J. Fluid Mech, vol 159, 1985, pp. 151-168.
[41] Triantafyllou G. S., Triantafyllou M. S., and Chryssostomidis C., “On the formation of
vortex streets behind stationary cylinders,” J. Fluid Mech. Vol. 170, 1986, pp. 461-477.
[42] Zhuang, M., Dimotakis, P. E. and Kubota, T. [1990], “The Effect of Walls on a Spatial
Growing Supersonic Shear Layer,” Physics of Fluids, 2(4), 599-604.
[43] Zhuang, M., Kubota, T. and Dimotakis, P. E. “On the instability of inviscid, compressible
free shear layers,” AIAA Journal, Vol. 28, No.10, 1990, pp. 1728-1733.
[44] Theofilis, V., Duck, P., and Owen, J., “Viscous Linear Stability Analysis of Rectangular
Duct and Cavity Flows”, J. Fluid Mech. Vol. 505, 2004, pp. 249-268.
[45] ANSYS Inc, “FLUENT 6.3 Manual”, 2008
[46] Pope, S., “Turbulent Flows”, Cambridge University Press, UK, 2000.
[47] BARDINA, J., FERZIER, J. H. and REYNOLDS, W. C., “Improved subgrid-scale models
for large eddy simulation”, AIAA paper, 78, 1980.
[48] Hinze, J., “Turbulence”, McGraw-Hill Publishing Co., New York, 1975.
[49] Nicoud, F. and Ducros, F., “Subgrid-Scale Stress Modeling Based on the Square of the
Velocity Gradient Tensor”, Flow, Turbulence, and Combustion, 62(3):183–200, 1999.
[50] Guleren, K., “Large-Eddy Simulation of Wall-Bounded Flows Subjected to Curvature and
Rotation”, PhD Dissertation, University of Manchester, 2007.
[51] Hodge, J. K., Stone, A. L., and Miller, T. E., “Numerical Solution for Airfoils near Stall in
Optimized Boundary-Fitted Curvilinear Coordinates,” AIAA Journal, Vol. 17, No. 5, 1979,
pp. 458-464.
[52] Burgmann, S., Schrö
der, W. Investigation of the vortex induced unsteadiness of a separation
bubble via time-resolved and scanning PIV measurements”. Experiments in Fluids, 45, 2008,
pp. 675-691.
[53] Greenberg, J.M., “Airfoil in Sinusoidal Motion in a Pulsating Stream,” N.A.C.A. TN-1326,
1947.
[54] Leishman, J. G., “Principles of Helicopter Aerodynamics”, Cambridge University Press,
2000.

137

[55] Ober, S., “Note on the Katzmayr effect on airfoil drag,” N.A.C.A. Technical Report No. 214,
1925.
[56] Sorensen, P., Hansen, A., Andre, P., Rosas, C., “Wind models for simulation of power
fluctuations from wind farms”. Journal of Wind Engineering and Industrial Aerodynamics,
Vol. 90, 2002, pp. 1381-1402.
[57] Gast, K., “A Comparison of Extreme Wind Events as Sampled in the 2002 Thunderstorm
Outflow Experiment”, Master’s Thesis, Texas Tech University, Lubbock.
[58] Hemida, H., Krajnovic, S., “Transient Simulation of the Aerodynamic Response of a
Double-Deck Bus in Gusty Winds”. Journal of Fluids Engineering, Vol. 131, 2009.

138