Grid and Time Step Requirements to Accurately &

Eciently Resolve Flow around a Rigid Flapping

Airfoil using OVERFLOW
Joshua I. Leell

Stanford University, Stanford, CA 94305

Thomas H. Pulliam

NASA Ames Research Center, Moett Field, CA 94305
Micro air vehicle design is an emerging area in aeronautics requiring accurate solutions
of apping wing ight in combination with numerical optimization. The highly unsteady
and complex nature of the ow elds associated with apping wings can translate to massive
computational costs. A large-scale optimization over a signicant design space can result
in hundreds to thousands of objective evaluations making it critical to ensure maximal
eciency while maintaining the requisite level of accuracy. Full-scale optimization can be
used for both the trajectory and geometry of the wings but this work focuses on a reduced
two-dimensional problem with xed geometry. A series of analyses are performed to isolate
the minimum spatial and temporal requirements to accurately and eciently model a two-
dimensional apping rigid airfoil subjected to sinusoidal pitching and plunging motion.
NASAs highly parallelized overset ow solver OVERFLOW is used for all of the computations
on a rigid NACA 0012 airfoil. Grid-spacing and domain size parameters are explored
using spatial analysis. The physical time-step with implicit sub-iterations in the context
of OVERFLOWs dual time-stepping scheme is explored using temporal analysis. Finally, a
method for evaluating periodic convergence in the force and moment time-history signals
is discussed.
Nomenclature
f Flapping frequency, Hz U
ref
Reference velocity
k Reduced frequency,
2fc
U
ref
t Physical time-step
h Nondimensional plunge amplitude N Number of time-steps per cycle,
2
kt
Maximum pitch amplitude NSI Sub-iterations per t
Pitch-plunge phase shift T Flapping period,
2
k
a Pitch axis location Dx Length of level-one grid in x-direction
y(t) Vertical position Dy Length of level-one grid in y-direction
(t) Rotation angle s Level one Cartesian grid spacing
c Chord length Re Reynolds number,
U
ref
c

n Iteration M
ref
Reference Mach number
c
l
Lift coecient c
d
Drag coecient
cm Moment coecient CCF Cross-correlation coecient
cT Time-averaged thrust coecient rms Root mean square
cP Time-averaged power-input coecient p Propulsive eciency,
c
t
cp
fM Fuzzy set for time-mean convergence fS Fuzzy set for convergence of signal shape
fA Fuzzy set for DFT amplitude convergence f

Fuzzy set for DFT phase convergence

Graduate Student, Department of Aeronautics & Astronautics and Research Scientist, STC, NASA Ames Research Center,
AIAA Student Member. jleell@stanford.edu

Senior Researcher, AIAA Associate Fellow. Thomas.H.Pulliam@nasa.gov

1 of 17
American Institute of Aeronautics and Astronautics
49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition
4 - 7 January 2011, Orlando, Florida
AIAA 2011-573
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
I. Introduction
A
signicant body of research has been performed analyzing periodically forced ows using both experimen-
tal and computational methods. Designers of micro air vehicles hope to exploit the highly maneuverable
aerodynamics exhibited by birds and insects in their work. A rst step in understanding the mechanisms that
lead to strong apping performance is reducing the three-dimensional problem with complex and deformable
geometry to a two-dimensional problem using simplied, rigid geometry. A sinusoidally forced airfoil can
produce either thrust or drag depending on the frequency and amplitude at which it plunges and pitches.
Jones, Dohring and Platzer [1] performed experiments with plunging airfoils and published results of a par-
ticular thrust- and lift-producing case. A plunging rigid NACA 0012 airfoil with reduced apping frequency,
k = 12.3, and nondimensional plunge amplitude, h = 0.12 results in the ow eld depicted in Figure 1a. The
reference Mach number for this case is M
ref
= 0.2 with Reynolds number, Re = 1850. Vorticity computed
by an OVERFLOW simulation using the same parameters is presented in Figure 1b. Snapshots of vorticity from
(a) Experimental results from Jones, Dohring & Platzer (b) OVERFLOW solution
Figure 1: Comparison of experimental and computed ow-elds
an OVERFLOW simulation at dierent intervals of the apping cycle are presented in Figure 2. The paired
vortices create thrust by producing a jet of uid in the downstream direction. Vertical position, y, and the
lift and drag signals are plotted in Fig. 3 at intervals corresponding to the labels of Figure 2. At station a
the airfoil is moving upwards with a lift decit. The airfoil has reached the apex of its trajectory by station
b with limited lift production. Maximum lift is attained at station c as the airfoil plunges rapidly down-
wards. Finally, the airfoil reaches the bottom of its stroke with a slight lift decit at station d. The plunging
conguration described and depicted above is used throughout the paper as an example of high-frequency
apping. The reader is directed to references [1-5] for more detailed discussion on the analysis of apping
airfoils.
Modeling ow over apping geometry can lead to signicantly higher computational cost compared to
its xed-body counterpart. Increased resolution in space is required because the body is moving through
its domain and may be shedding ow structures such as vortices downstream. Finer temporal resolution
is needed to maintain consistency with the unsteady physics of the problem and avoid interpolation errors
from excessive grid motion between each time-step. Other costs arise from the need to simulate over a
number of apping cycles until the ow eld converges to a periodic state achieving functional invariance
of the time-averaged performance measures over the period of a apping cycle. The ability to isolate the
minimum discretization parameters and stop the simulation upon convergence will signicantly reduce the
computational burden. This computational work is amplied in the context of unsteady optimization where
hundreds to thousands of ow solutions are required. The goal of this paper is to understand the spatial
and temporal discretization requirements for the reduced two-dimensional problem and establish a reliable
method for measuring convergence of the periodic ow to reduce the computational cost of optimization.
The spatial investigation includes a grid-spacing study in addition to analysis of solution sensitivity to the
domain size. The temporal analysis examines solution dependence on the total number of iterations used per
apping cycle. Finally, a method for analyzing convergence of a periodically forced ow eld is introduced
and applied. Analyses performed in two-dimensions will hopefully aid future work in three-dimensions where
both computational cost and complexity increase signicantly.
2 of 17
American Institute of Aeronautics and Astronautics
(a) t/T =
1
4
(b) t/T =
1
2
(c) t/T =
3
4
(d) t/T = 1
Figure 2: Flow eld visualization of vorticity magnitude
y
t/T
0 0.2 0.4 0.6 0.8 1
-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
d
c
b
a
y
t/T
(a) y
c
l
t/T
0 0.2 0.4 0.6 0.8 1
-80
-60
-40
-20
0
20
40
60
a
b
c
d
c
l
t/T
(b) c
l
c
d
t/T
0 0.2 0.4 0.6 0.8 1
-2.5
-2
-1.5
-1
-0.5
0
0.5
d
c
b
a
c
d
t/T
(c) c
d
Figure 3: y, c
l
and c
d
at corresponding stations
3 of 17
American Institute of Aeronautics and Astronautics
II. Overow
OVERFLOW is a highly parallelized overset nite-dierence Navier-Stokes ow solver capable of handling grid
motion such as apping airfoils. The code oers several time integration methods and spatial schemes but all
simulations presented use second-order implicit dual time-stepping and central dierencing in space providing
a sixth-order accurate convection operator and a fth-order accurate dissipation operator. OVERFLOWs overset
framework uses a structured body-conforming grid embedded in a set of automatically generated Cartesian
o-body grids. X-Rays are used to provide hole-cutting [6] with domain connectivity and interpolation across
grid interfaces [7]. A more detailed explanation of the specics of OVERFLOW can be found in the paper by
Nichols, Tramels and Buning [8].
Figure 4 shows an example of a typical two dimensional OVERFLOW grid system (reduced for the purpose
of illustration). The curvilinear body-conforming mesh (as seen in orange in Fig. 4a), referred to as the
near-body grid, is embedded in a series of Cartesian o-body grids. The rst Cartesian grid is referred to as
the level-one grid with domain size D
x
by D
y
with uniform grid-spacing of s in both the x- and y-directions
(see Fig. 4). The series of Cartesian o-body grids (the second grid is referred to as the level-two grid and
so on) span the domain to the outer boundaries with each level doubling the previous grids spacing. Each
successive grid level has an overlap region of at least two grid cells of the ner spacing as demonstrated in
Fig. 4a between the red and blue hashed lines. The X-Rays corresponding to the near-body grid cut the
level-one grid and dene an interface through which information will be interpolated between the two grids.
Fifth-order spatial discretization requires the use of three fringe, or overlap, points. Figure 4b outlines the
X-Ray hole cutting in the upper surface, leading edge region, with a hashed red line.
(a) Level-one grid domain (b) X-ray hole cutting
Figure 4: Grid denitions
The solution procedure for a given apping conguration requires the near-body grids, the X-Rays,
surface integration information and an OVERFLOW input le that denes the boundary conditions and essential
parameters for the simulation. Two XML les control the grid dynamics. One moves the geometry into its
initial position and provides grid-name associativity and the other prescribes the translational and rotational
velocities of the associated components. Force and moment coecients are easily extracted and used to
integrate the time-averaged thrust and power-input coecients as demonstrated in Section IV. Far-eld
boundaries are placed at 32 chord-lengths in both the positive and negative x- and y-directions for all
simulations presented.
III. Kinematics of a Rigid Pitching & Plunging Airfoil
Equations 1 and 2 dene the kinematic model of a sinusoidally pitching and plunging rigid airfoil with a
frequency of f Hz.
y (t) = hcos (2ft) (1)
(t) = cos (2ft +) (2)
4 of 17
American Institute of Aeronautics and Astronautics
The airfoils vertical position, y(t) and the rotation angle, (t), are fully dened by time, t, and ve parameters
of interest: the reduced frequency, k =
2fc
U
ref
, the nondimensional plunge amplitude, h, the maximum pitch
amplitude, , the phase shift between pitching and plunging, and the pitch-axis along the chord line, a.
(t)
y(t)
a
Figure 5: Kinematics of apping airfoil
These parameters can be used as design variables for optimizing the trajectory of apping airfoil. A
schematic of the kinematics of a rigid apping airfoil is shown in Figure 5. Full-scale optimization may
include the geometry of the apping body but a NACA 0012 airfoil is used for all of the cases presented.
IV. Analysis of Temporal & Spatial Requirements
Aerodynamic performance can be measured in any number of ways. Force and moment coecients work
well for steady-state problems but apping wing solutions produce periodic time histories of these functionals
without a single value to compare directly with other cases. Computing the time-averaged mean thrust and
power-input coecients, c
T
and c
P
, provides the analyst with scalar performance measures over a apping
cycle. The ratio of these scalars is dened as the propulsive eciency of the apping conguration,
P
.
Equations 3-5 dene each of these functionals.
c
T
=
1
T
_
t+T
t
c
d
(t)dt (3)
c
P
=
1
T
_
t+T
t
_
c
l
(t) y(t) +c
m
(t)

(t)
_
dt (4)

P
= c
T
/ c
P
(5)
A tness function to measure apping performance was constructed by Tuncer and Kaya [9] as a weighted
combination of c
T
and
P
. A generalized objective function, O(X, S, C), is the convex combination of c
T
and
P
.
O(X, S, C) = c
T
(X, S, C) + (1 )
P
(X, S, C) where 0 1 (6)
The independent variables X = (k, h, , , a) represents the design variables of a given apping congu-
ration, S = (s, D
x
, D
y
) represents the spatial discretization and C = (t, N
SI
) represents the temporal
discretization. Lift, drag and moment coecients are extracted from the ow solution and y(t) and

(t)
are computed analytically by dierentiating the expressions for y(t) and (t) in Equations 1 and 2. Further
details on optimization of apping airfoils are provided by Oyama, Okabe, Fujii and Shimoyama [10].
While the objective function is evaluated with the time-averaged integrated functionals of c
T
, c
P
and

P
, it is critical to examine how the signal shape of the time-history of the integrated force and moment
coecients dier as the spatial and temporal resolutions are adjusted. Two solutions may share a similar
thrust coecient but exhibit signicant shape, magnitude and phase dierences in the time-history response
of the drag coecient, for example. Two similar signals are guaranteed to have similar time-averaged
functionals but similar time-averaged functionals are not guaranteed to have similar signal shape. Clark
and Grover [11] dene the cross-correlation coecient, CCF(f, g), which qualies how well the shape of one
signal, f, matches another, g.
CCF (f, g) =
N

n=1
f (n) g (n)
_
N

n=1
f
2
(n)
N

n=1
g
2
(n)
_
1
2
(7)
5 of 17
American Institute of Aeronautics and Astronautics
The value of the cross-correlation coecient ranges from, 0 CCF 1, where 0 describes no correlation
and 1 indicates complete correlation identical signals. CCF is used to check similarity of a signal from
period to period, g (n) = f (n +N), or to check how two individual signals, f and g, compare.
A second measure of signal correlation between two responses compares the root mean square, rms, of
each signal over one apping cycle. This proves a better tool than pure averaging because the force and
moment coecients may vary between positive and negative values over a apping cycle. Equation 8 denes
the root mean square for a signal, f.
rms (f) =

_
1
N
N

n=1
f (n)
2
(8)
It is important to verify that the analytical tools successfully perform their tasks. Initially, percentage error
in thrust coecient was used to compare signals but this proved inaccurate because certain frequencies
integrated thrust approached zero. Results to be presented later conrm that the quantitative assessments
of signal correlation using CCF and rms values agree with qualitative observations.
IV.A. Spatial Requirements
Two spatial discretization factors which greatly inuence solutions of apping airfoil ow are: a) the grid
spacing near the body and in the region traversed by shed vortices, and b) the size of most rened region
the level-one grid. A brute-force approach suggests solving on an arbitrarily large domain with exceedingly
ne grid spacing. This naturally leads to prohibitive computational costs especially in the context of
optimization that may demand hundreds to thousands of unsteady ow solutions. Grid selection must isolate
the minimum spatial requirements to obtain solutions with the requisite level of accuracy. The following two
sections explore the sensitivity of apping airfoil solutions to near-body grid-spacing and level-one domain
size over a range of apping frequencies.
IV.A.1. Grid Spacing
Four near-body grids, with decreasing resolution, are used to determine grid independence for apping airfoil
ow solutions. The goal is to select a near-body grid-spacing that successfully resolves solutions for apping
congurations over a spectrum of apping frequencies. An optimizer will require solutions to a variety
of apping congurations and it is unreasonable to determine the spatial requirements for each objective
function evaluation a priori. The spacing of the level-one grid is chosen to match the spacing of the near-
body grid in the overlap region to ensure high-quality interpolation between the two grids. Table 1 provides
information about each of the grids where index j runs tangential to the airfoil and index k runs normal to
the airfoil in a structured c-mesh. The rst point normal to the airfoil in k is 0.0036c for each grid.
Grid # n
j
n
k
n
Airfoil
L1 Pts s
1 1601 201 1025 6,352,020 0.0025c
2 801 151 513 1,603,305 0.005c
3 401 101 257 407,721 0.01c
4 201 51 129 71,145 0.025c
Table 1: Near-body grid points in j, in k, along the airfoil, level-one (L1) grid-points and grid-spacing
The objective functional for pure plunging is computed based on the lift and drag coecients. Of these,
the drag coecient is the slowest to converge and as such it will be used to determine grid independence.
Figure 6 shows the c
d
signal of the 24th (re: converged) apping cycle for the four test case frequencies. The
highest frequency cases use the conguration outlined in Section I. The frequency of k = 12.3 is halved three
times arriving at the four test frequencies: k = 1.5375, 3.075, 6.15 and 12.3. All cases are pure plunging
with h = 0.12. The level-one grid domain size is D
x
D
y
= 13c 3c with the airfoil placed such that
3
4
of the level-one domain remains downstream. Each case uses 2000 time-steps per cycle (N) and 10 implicit
sub-iterations (N
SI
) as dened in Section IV.B.
6 of 17
American Institute of Aeronautics and Astronautics
t/T
c
d
Grid 1
Grid 2
Grid 3
Grid 4
0.5 0.6 0.7 0.8 0.9 1
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.15
0.16
0.17
(a) k = 1.5375
t/T
c
d
Grid 1
Grid 2
Grid 3
Grid 4
0.5 0.6 0.7 0.8 0.9 1
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
(b) k = 3.075
t/T
c
d
Grid 1
Grid 2
Grid 3
Grid 4
0.5 0.6 0.7 0.8 0.9 1
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
(c) k = 6.15
t/T
c
d
Grid 1
Grid 2
Grid 3
Grid 4
0.5 0.6 0.7 0.8 0.9 1
-2.5
-2
-1.5
-1
-0.5
0
(d) k = 12.3
Figure 6: c
d
signal for dierent apping frequencies, k, over half a apping period, T
The cross-correlation coecient and percent error in root mean square of the drag signal between grids
2-4 and the nest grid (Grid 1) are presented in Table 2. Error in rms, e
rms
, provides more meaningful results
than error in thrust coecient because c
T
is an integral of positive and negative values. If the integration
results in a small number (as is the case for the for the k = 3.075 conguration) then the quantitative error
calculation can dier greatly from the observed qualitative error. Equation 9 denes percent error in rms
for a signal f.
%e
rms
= 100

rms (f) rms (f

exact
)
rms (f
exact
)

(9)
In Table 2, CCF and e
rms
are computed for each grid taking the c
d
signal of the nest grid as exact.
7 of 17
American Institute of Aeronautics and Astronautics
Grid 2 Grid 3 Grid 4
k CCF e
rms
CCF e
rms
CCF e
rms
1.5375 0.9999 0.171% 0.9999 0.407% 0.9999 1.75%
3.075 0.9999 0.0834% 0.9995 0.472% 0.9986 0.720%
6.15 0.9999 0.427% 0.9985 2.81% 0.9986 3.74%
12.3 0.9999 0.280% 0.9973 3.13% 0.8284 7.60%
Table 2: c
d
signal comparison for each test frequency
The signal quality of c
d
is satisfactory for the three coarser grids at the lower frequencies of k = 1.5375,
3.075 and 6.15. However, the correlation of c
d
signal for k = 12.3 on the coarsest grid (Grid 4) is unacceptable
because of substantial magnitude and phase errors. The error estimates agree with observed results in Fig. 6d
where CCF = 0.8284 and e
rms
= 7.60%. Grid 3 suciently resolves the drag signal for each test frequency
and will thus be used for much of the remaining analysis to save considerable computational cost over the
two ner grids.
IV.A.2. Domain Size
A second task in grid selection is determining an appropriate level-one domain size. The solution at every
grid point in the entirely subsonic domain inuences the solution at every other point. Thus, the solution
at every point in the domain contributes to the force and moment functional values. The signicance of the
contribution, however, is unknown but it is clear that points near the boundary are important warranting
an investigation into the solution sensitivity of level-one grid size.
Label C1 C2 C3 C4 U1 U2 U3 U4 Control
D
x
5c 9c 17c 33c 5c 9c 17c 33c 33c
D
y
4c 4c 4c 4c 4c 4c 4c 4c 16c
Points 200k 360k 680k 1320k 200k 360k 680k 1320k 5280k
Table 3: Level-one domain size information
Two apping congurations are solved on a series of level one domains to investigate the solution sensi-
tivity to domain size and the location of the airfoil within that domain. A low-frequency case of combined
pitching and plunging and a high-frequency case of pure plunging are included to derive insight into the
domain requirements throughout a reasonable design space of apping parameters, X. Additionally, the
location of the airfoil within the domain is included in the analysis. In the cases examined, large vortical
structures produced by the apping motion are convected downstream of the airfoil. This suggests placing
the airfoil upstream of the domains center to resolve more of these structures on a ner mesh for a longer
period of time. The sizes of the level-one domain used are dened in Table 3. Grid labels beginning with a
C- prex indicate central placement of the airfoil whereas a U- prex refers to upstream airfoil placement. A
control domain is used to compute the CCF and e
rms
. It extends 33 chord-lengths in x and 16 chord-lengths
in y with the airfoil in the center resulting in a level-one grid with approximately 5.3 million grid-points. For
upstreamed-centered cases, the airfoil is placed such that
3
4
of the level-one domain remains downstream.
Grid 3 from the previous section is used as the near-body grid with s = 0.01, N = 2000 physical time-steps
per cycle and N
SI
= 10 implicit sub-iterations per time-step. Signals are extracted from the 10th apping
period.
8 of 17
American Institute of Aeronautics and Astronautics
t/T
c
d
C1
C2
C3
C4
U1
U2
U3
U4
Control
0.5 0.6 0.7 0.8 0.9 1
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
(a) Low-frequency pitching & plunging
k = 1, h = 0.5, = 5

, = 30

t/T
c
d
C1
C2
C3
C4
U1
U2
U3
U4
Control
0.5 0.6 0.7 0.8 0.9 1
-2
-1.8
-1.6
-1.4
-1.2
-1
-0.8
-0.6
-0.4
-0.2
(b) High-frequency plunging
k = 12.3, h = 0.12, = 0

, = 0

Figure 7: c
d
dependence on level-one grid domain size and airfoil location
Figure 7 demonstrates the increased sensitivity to level-one domain size for the high-frequency case. The
smallest domain (D
x
= 5c) was sucient for resolving the low frequency case but not the high frequency
plunging case. Tables 4 and 5 quantitatively demonstrate the increased domain size dependence for the high
frequency case. CCF and rms values compare each grid to the control.
Center Upstream
D
x
D
y
CCF e
rms
CCF e
rms
5c 4c 0.9996 0.0905% 0.9998 0.172%
9c 4c 0.9998 0.441% 0.9999 0.365%
17c 4c 0.9998 0.0111% 0.9999 0.324%
33c 4c 0.9999 0.0455% 0.9999 0.125%
Table 4: Signal comparison for k = 1, h = 0.5, = 5

, = 30

Center Upstream
D
x
D
y
CCF e
rms
CCF e
rms
5c 4c 0.9989 3.78% 0.9986 1.66%
9c 4c 0.9998 1.56% 0.9999 0.561%
17c 4c 0.9999 0.586% 0.9999 0.647%
33c 4c 0.9999 0.598% 0.9999 0.594%
Table 5: Signal comparison for k = 12.3, h = 0.12, = 0

, = 0

The smallest domain, with the airfoil centered or located upstream, is satisfactory for the lower-frequency
case while the high-frequency case demands a domain that spans at least 9 chord lengths in the x-direction.
Placing the airfoil upstream improves the solution for the smallest and second to smallest domains as it
allows for ner resolution of vortices in the proximity of the airfoil. There is still much to be explored in
establishing grid requirements for forced periodic ows but proximity- and feature-based adaption will assist
in automating the process by dynamically generating new grids to resolve desired features such as vortices
convecting throughout the domain.
9 of 17
American Institute of Aeronautics and Astronautics
IV.B. Temporal Requirements
OVERFLOWs use of implicit dual time-stepping requires selection of a physical time-step, t, and the number
of implicit sub-iterations, N
SI
, to be performed at each time-step. The combination of t and N
SI
must
suciently model the physics of the problem while reducing the residual at each time-step by a reasonable
amount to ensure a commensurate level of accuracy. There are additional temporal parameters built into the
code such as maximum CFL and local time-step but discussion here is limited to t and N
SI
. OVERFLOWs
implicit dual time-stepping scheme, described in Eq. 10, uses both physical time, t, and computational time,
, to iterate the solution forward.
dQ
d
+
dQ
dt
+F (Q(t)) = 0 (10)
where Q is a vector of the ow-variables and F represents the inviscid and viscous ux terms of the Navier-
Stokes equations. Applying second-order implicit dierencing in t and rst-order implicit dierencing in
,
Q
s+1
Q
s

+
_
3Q
s+1
4Q
n
+Q
n1
2t
+F
_
Q
s+1
_
_
= 0 (11)
Linearizing in s results in Eq. 12 where L(Q
s
) is the implicit linear driver with sub-iterations advancing in
s, i.e. s = 1, . . . , N
SI
.
L(Q
s
)
_
Q
s+1
Q
s
_
=
_
3Q
s
4Q
n
+Q
n1
2t
+F (Q
n
)
_
= R(Q
s
, Q
n
) (12)
R(Q
s
, Q
n
) in Eq. 12 is the residual at time-step n and sub-iteration s. The norm of R, ||R||, measures
sub-iteration convergence. In the context of optimization, temporal resolution must be selected such that
computational resources are not wasted by solving more time-steps per apping cycle than necessary or
driving the residual beyond a necessary level with superuous implicit sub-iterations. A parametric study is
performed to determine appropriate values for t and N
SI
. The low- and high-frequency apping congu-
rations from the domain size study in Sec. IV.A.2 are used. For consistency among the dierent frequencies
t is dened as:
t =
2
kN
(13)
Results using a range of N and N
SI
values are computed to demonstrate solution trends as the number of
total iterations per apping cycle, N N
SI
, increases to a maximum value of 128,000. Figures 8 and 9 show
the root mean square of the drag signals for the low- and high-frequency cases, respectively. The root mean
square of the drag coecient is plotted as a function of total iterations per cycle. Figures 8a and 9a plot the
root mean square of drag as a function of total iterations per apping cycle for dierent values of N where
N
SI
increases for each data-point moving from left to right. Figures 8b and 9b plot the root mean square
of drag as a function of total iterations per apping cycle for dierent values of N
SI
where N increases for
each data-point moving from left to right. Signals are extracted from the 10th apping period.
10 of 17
American Institute of Aeronautics and Astronautics
N x NSI
r
m
s
(
c
d
)
N = 250
N = 500
N = 1000
N = 2000
N = 4000
10
3
10
4
10
5
10
6
0.24
0.245
0.25
0.255
0.26
0.265
0.27
= 1.6%
(a) N vs N N
SI
N x NSI
r
m
s
(
c
d
)
NSI = 2
NSI = 4
NSI = 8
NSI = 16
NSI = 32
10
3
10
4
10
5
10
6
0.24
0.245
0.25
0.255
0.26
0.265
0.27
(b) N
SI
vs N N
SI
Figure 8: rms(c
d
) for k = 1, h = 0.5, = 5

, = 30

Figures 8a and 8b show that using most values of N and N

SI
will converge as the total number of
iterations per cycle increases. Figure 8a shows that combinations of total iterations per cycle using N = 250
fail to completely converge in the limit of N N
SI
. The percentage error in rms for N = 250 versus
N = 4000 for the maximum number of total iterations is 1.6%. Comparatively, the error in rms for N
SI
= 2
versus N
SI
= 32 is only 0.298% for the maximum number of total iterations. These plots show that once a
certain number of iterations is reached, N N
SI
2 10
4
(see dotted red lines), the specic values of N
and N
SI
are free to change contingent upon N 500. Using N = 250 results in a physical time-step that
may be too large to capture all of the unsteady physics of the ow.
N x NSI
r
m
s
(
c
d
)
N = 250
N = 500
N = 1000
N = 2000
N = 4000
10
2
10
3
10
4
10
5
10
6
1.05
1.06
1.07
1.08
1.09
1.1
1.11
1.12
1.13
1.14
= 2.3%
(a) N vs N N
SI
N x NSI
r
m
s
(
c
d
)
NSI = 2
NSI = 4
NSI = 8
NSI = 16
NSI = 32
10
2
10
3
10
4
10
5
10
6
1.05
1.06
1.07
1.08
1.09
1.1
1.11
1.12
1.13
1.14
(b) N
SI
vs N N
SI
Figure 9: rms(c
d
) for k = 12.3, h = 0.12, = 0

, = 0

Results for the high-frequency case presented in Fig. 9 oers similar temporal guidelines for the low-
frequency case. Figures 9a and 9b show that solutions using N = 250 will not converge to a consistent
value in the limit of total iterations per cycle. It appears that this apping conguration also requires
approximately 2 10
4
total iterations for all solutions to fully converge. The error in rms for N = 250
for the maximum of total iterations per cycle is 2.3% whereas the error in rms for N
SI
for the maximum
number of total iterations is only 0.0056%. Again, a simulation using too large a physical time-step will
never converge in the limit of temporal resolution. Using a small number of implicit sub-iterations with a
large number of time-steps per cycle will converge for both apping congurations at which point the scheme
essentially reverts back to a single time-stepping method. These examples only demonstrate convergence in
11 of 17
American Institute of Aeronautics and Astronautics
the limit of temporal resolution and do not address overall accuracy of the solution.
While determination of a threshold for iterations required per cycle can serve as a guideline it should not
be followed dogmatically as temporal convergence can be achieved with signicantly less resolution.
t/T
c
d
N = 250, NSI = 8
N = 500, NSI = 8
N = 1000, NSI = 8
N = 2000, NSI = 8
N = 4000, NSI = 8
0 0.2 0.4 0.6 0.8 1
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
(a) c
d
for constant N
SI
= 8
t/T
c
d
N = 250, NSI = 512
N = 500, NSI = 256
N = 1000, NSI = 128
N = 2000, NSI = 64
N = 4000, NSI = 32
0 0.2 0.4 0.6 0.8 1
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
(b) c
d
for maximum N N
SI
= 128, 000
Figure 10: Comparison of c
d
signal between N
SI
= 8 and maximum iterations for k = 1, h = 0.5, = 5

,
= 30

Figure 10 shows that the drag signal for the low-frequency case fails to converge to if fewer than 1000
time-steps per cycle are combined with only 8 sub-iterations but that all combinations converge in the limit
of temporal resolution with limited exception to the N = 250 case. Thus, at least 8000 total iterations per
apping cycle are required to resolve the low-frequency case. Figure 11b shows that the drag signals for all
but N = 250 converge in the limit of iterations per apping cycle for the high-frequency case but Fig. 11a
shows that convergence in the drag signal is achieved for N = 500 with only eight sub-iterations totaling
only 4000 iterations per apping cycle and less than the low-frequency case. Pitching and increased plunging
amplitude in the low-frequency case necessitates increased temporal resolution.
t/T
c
d
N = 250, NSI = 8
N = 500, NSI = 8
N = 1000, NSI = 8
N = 2000, NSI = 8
N = 4000, NSI = 8
0 0.2 0.4 0.6 0.8 1
-2.5
-2
-1.5
-1
-0.5
0
(a) c
d
for constant N
SI
= 8
t/T
c
d
N = 250, NSI = 512
N = 500, NSI = 256
N = 1000, NSI = 128
N = 2000, NSI = 64
N = 4000, NSI = 32
0 0.2 0.4 0.6 0.8 1
-2.5
-2
-1.5
-1
-0.5
0
(b) c
d
for maximum N N
SI
= 128, 000
Figure 11: Comparison of c
d
signal between N
SI
= 8 and maximum iterations for k = 12.3, h = 0.12
The total iterations per cycle required for these cases are signicantly less than the threshold discussed
above and a small fraction of the maximum iterations used. The threshold can be used as an initial estimate
but less iterations have proven to be sucient for both test cases.
12 of 17
American Institute of Aeronautics and Astronautics
It would be benecial to monitor the solution during the computation to determine if adjusting temporal
parameters could reduce the wall-time necessary to reach a periodic unsteady convergence. One way to
achieve this would be reducing the number of time-steps per cycle or implicit sub-iterations. Another
approach would increase the number of time-steps per cycle with the goal of reducing the number of apping
cycles for the solution to enter a periodic state.
V. Periodic Unsteady Convergence
Solutions of apping airfoils require the computation of multiple apping cycles before the ow eld
reaches a periodic state. The work required to reach this periodic state is dependent upon the kinematic
parameters of the apping conguration in addition to spatial and temporal discretization. Because the ow
is periodically forced by apping the airfoil, the force and moment functionals will converge to a periodic
signal as opposed to xed values. Clark and Grover [11] refer to this type of convergence as periodic unsteady
convergence. Terminating the computation upon rst reaching unsteady periodic convergence can provide
signicant computational savings especially in the context of optimization.
The objective function, dened in Eq. 6, is a convex combination of the time-averaged thrust and
propulsive eciency coecients. Monitoring convergence of these two functionals, or the objective function
directly, is straightforward. As previously mentioned, similar values of the thrust and power coecients can
be generated from dissimilar force and moment time-history signals. A more robust method of convergence
estimation must be employed to evaluate the correlation of functional signal shape from one cycle to the
next to guarantee a periodically converged ow eld.
Clark and Grover [11] present one such framework for assessing periodic unsteady convergence. A method-
ology is developed for computing a global measure of periodic unsteady convergence by creating a series of
fuzzy sets that each describe membership into unique elements of periodic unsteady convergence. Clark and
Grover [11] use pressure signals from an array of sample-points in the domain to monitor time-averaged
uctuations, Discrete Fourier Transform (DFT) magnitudes and phases, cross-correlation coecients and
power spectral densities at each apping cycle. Values from the current and previous cycle are compared to
determine the degree of convergence in each measure after every apping period. Clark and Grover [11] em-
ploy multi-valued logic to produce an additional fuzzy set represented by a single value, f
c
, describing overall
convergence of the unsteady periodic solution. Fuzzy sets for mean-level, amplitude, phase angle and overall
signal shape are employed. Integrated pressure on the surface of the airfoil provides time-history signals for
the force and moment coecients. Because the force and moment time-history signals vary between positive
and negative values, the root-mean-square is used in the time-mean fuzzy sets instead of pure averaging.
CCF and rms are dened previously in Eqs. 7 and 8, respectively. Equations 14-15 dene the magnitude,
A, and phase, , for a DFT represented by P.
A = 2
|P|
N
(14)
= tan
1
_
(P)
(P)
_
(15)
Membership grade in a fuzzy set varies between 0 and 1, 0 f 1, where 0 suggests no membership in the
fuzzy set f and 1 indicates complete membership. Unlike boolean logic, partial membership in fuzzy sets
is permitted. More information on the theory and application of fuzzy sets can be found in Klir and Yuan
[12]. Equations 16-19 dene the fuzzy sets for mean-value, DFT amplitude, DFT phase and signal shape,
respectively.
f
M
= 1

1
rms (c
2
)
rms (c
1
)

(16)
f
A
= 1

1
A
2
A
1

(17)
f

= 1

(18)
f
S
= CCF (c
1
, c
2
) (19)
The subscripts of 2 and 1 in the fuzzy set denitions in Eqs. 16-19 indicate the current and previous apping
cycle, respectively. The c refers to a force or moment coecient signal. In practice, three fuzzy sets, f
M,c
l
,
13 of 17
American Institute of Aeronautics and Astronautics
f
M,c
d
and f
M,cm
, are dened to monitor the root-mean-square convergence for the lift, drag and pitching
moment coecients at each cycle. Three additional fuzzy sets, f
S,c
l
, f
S,c
d
and f
S,cm
, are constructed to
monitor the signal shape of the three force and moment coecients using the cross-correlation coecient.
The spectral-based fuzzy sets of f
A
and f

are computed for the fundamental frequency of each force and

moment coecient. In summary, there are a total of three instances (corresponding to c
l
, c
d
and c
m
) of each
of the four fuzzy set categories, f
M
, f
A
, f

and f
S
.
The values of Eqs. 16-19 dene the degree of membership in each fuzzy set. Global convergence, f
c
, is
treated as a fuzzy set characterized by the intersection of the other sets.
f
c
= f
M
f
A
f

f
S
(20)
The intersection of fuzzy sets are dened as
f
c
= min (f
M
, f
A
, f

, f
S
) (21)
Thus, overall convergence is evaluated as the weakest membership grade in any sub-global set after each
cycle. Clark and Grover [11] suggest that two consecutive cycles with f
c
0.95 implies convergence of the
periodic unsteady ow eld.
A time-history of the membership grade for each of the 12 non-global fuzzy sets for the high-frequency
case is presented in Figure 12a. The coarsest grid, Grid 4, from Sec. IV.A.1 was used with N = 2000 and
N
SI
= 10. Fuzzy-sets corresponding to the drag signal are the slowest to reach the tolerance level of 0.95
designated by the dotted black line. The corresponding drag coecient time-history is plotted in Fig. 12b
for comparison.
f
t/T
0 5 10 15 20 25
0
0.5
1
(a) Non-global fuzzy sets
fM,cl
fS,cl
fA,cl
f,cl
fM,cd
fS,cd
fA,cd
f,cd
fM,cm
fS,cm
fA,cm
f,cm
20
c
d
t/T
0 5 10 15 20 25
-3
-2
-1
0
(b) c
d
Figure 12: Time history of drag coecient, c
d
, and non-global fuzzy sets for k = 12.3, h = 0.12
Each of the four test-cases described in Sec. IV.A.1 are analyzed for periodic unsteady convergence over
a range of time-steps per cycle: N = (250, 500, 1000, 2000, 4000). Ten implicit sub-iterations, N
SI
, are used
14 of 17
American Institute of Aeronautics and Astronautics
in each case. The three lower-frequency congurations converge rapidly but the high-frequency conguration
requires signicantly more cycles to converge. The iterations required for each of the test-cases to converge
is provided in Table 6. The time-histories of the thrust coecients and global fuzzy sets, f
c
, using N = 2000
are presented in Figures 13-16 to demonstrate typical evolution of the drag signal for each case.
k N Grid 1 Grid 2 Grid 3 Grid 4
1.5375
250 3 3 3 3
500 3 3 3 3
1000 3 3 3 3
2000 3 3 3 3
4000 3 3 3 3
3.075
250 3 3 3 3
500 3 3 3 3
1000 3 3 3 3
2000 3 3 3 3
4000 3 3 3 3
6.15
250 3 3 3 3
500 3 3 3 3
1000 3 3 3 3
2000 3 3 3 3
4000 3 3 3 3
12.3
250 8 8 8 19
500 8 9 7 11
1000 8 11 8 13
2000 8 11 8 13
4000 8 11 8 13
Table 6: Number of apping cycles computed before reaching periodic unsteady convergence for each grid
at each frequency and for each number of time-steps per cycle
t/T
c
T
Grid 1
Grid 2
Grid 3
Grid 4
0 5 10 15 20 25
-0.16
-0.155
-0.15
-0.145
-0.14
-0.135
-0.13
-0.125
-0.12
-0.115
(a) c
T
t/T
f
c
Grid 1
Grid 2
Grid 3
Grid 4
Tolerance = 0.95
0 5 10 15 20 25
0.4
0.5
0.6
0.7
0.8
0.9
1
(b) fc
Figure 13: Time history of thrust coecient, c
T
, and global fuzzy set, f
c
for k = 1.5375, h = 0.12
15 of 17
American Institute of Aeronautics and Astronautics
t/T
c
T
Grid 1
Grid 2
Grid 3
Grid 4
0 5 10 15 20 25
-0.08
-0.07
-0.06
-0.05
-0.04
-0.03
-0.02
-0.01
(a) c
T
t/T
f
c
Grid 1
Grid 2
Grid 3
Grid 4
Tolerance = 0.95
0 5 10 15 20 25
0.4
0.5
0.6
0.7
0.8
0.9
1
(b) fc
Figure 14: Time history of thrust coecient, c
T
, and global fuzzy set, f
c
k = 3.075, h = 0.12
t/T
c
T
Grid 1
Grid 2
Grid 3
Grid 4
0 5 10 15 20 25
0.3
0.32
0.34
0.36
0.38
0.4
0.42
(a) c
T
t/T
f
c
Grid 1
Grid 2
Grid 3
Grid 4
Tolerance = 0.95
0 5 10 15 20 25
0.7
0.75
0.8
0.85
0.9
0.95
1
(b) fc
Figure 15: Time history of thrust coecient, c
T
, and global fuzzy set, f
c
for k = 6.15, h = 0.12
t/T
c
T
Grid 1
Grid 2
Grid 3
Grid 4
0 5 10 15 20 25
0.9
0.95
1
1.05
1.1
1.15
1.2
1.25
1.3
(a) c
T
t/T
f
c
Grid 1
Grid 2
Grid 3
Grid 4
Tolerance = 0.95
0 5 10 15 20 25
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(b) fc
Figure 16: Time history of thrust coecient, c
T
, and global fuzzy set, f
c
for k = 12.3, h = 0.12
The lower-frequency cases shown in Figures 13-15 converge rapidly for all values of N and on each level
of grid resolution while f
c
for the higher frequency case uctuates for each grid resolution. The spectral
16 of 17
American Institute of Aeronautics and Astronautics
fuzzy sets, f
A
and f

, for drag are the most resistant to convergence in this case.

Periodic unsteady convergence analysis provides key insights into the behavior of solutions with dierent
apping parameters and spatial and temporal discretizations but it does not account for accuracy. Accuracy
is controlled by the spatial and temporal discretization discussed in the previous section, whereas, convergence
is primarily controlled by the apping conguration. This is exemplied by the slow convergence for the
high-frequency case for all temporal and spatial discretizations. The overall goal is to balance convergence
rate with sucient accuracy.
VI. Summary
This work focused on understanding how spatial and temporal discretization aect solutions to apping
airfoil simulations over a range of apping parameters in the context of OVERFLOWs overset grid and dual
time-stepping frameworks. Informative signal processing error estimators were introduced and veried. Two
types of spatial analysis were performed that focused on grid-spacing in near-body grid and domain size of
the level-one gird. Temporal analysis investigated the interplay between time-steps per cycle and implicit
sub-iterations per time-step and suggests that solutions converge in the limit of iterations per cycle for all
but the largest physical time-steps. Insight gained studying the two-dimensional process will be applied
to three-dimensional cases where selection of controlling parameters becomes less obvious. This was by no
means an exhaustive study but key trends in solution sensitivity to spatial and temporal discretizations were
uncovered.
A method to measure periodic convergence was applied based on the work of Clark and Grover [11]. This
demonstrated that convergence rate is highly dependent on the apping conguration. Future plans include
development of a monitoring process to analyze solution convergence during the simulation. An additional
process aimed at controlling the convergence of a simulation by dynamically adjusting the quantity of sub-
iterations or the time-step size will also be addressed. All of these eorts are directed at the development of
an ecient framework for optimization of forced periodic ows.
Acknowledgments
This work is funded by the Army High Performance Computing Research Center (AHPCRC). The authors
would like to thank NASA Advanced Supercomputing for use of their facilities.
References
1
Jones, K.D., Dohring, C. M., and Platzer, M. F. 1998. Experimental and Computational Investigation of the Knoller-Betz
Eect, AIAA Journal, 36(7): 1240-1246
2
Young, J. and Lai, J. C. S. 2007. Vortex Lock-In Phenomenon in the Wake of a Plunging Airfoil, AIAA Journal, 45(2):
485-490
3
Lai, J. C. S., and Platzer, M. F. 1999. Jet Characteristics of a Plunging Airfoil, AIAA Journal, 37(12):1529-1537
4
Platzer, M.F, and Jones, K.D. 2006. Flapping Wing Aerodynamics - Progress and Challenges, 44th AIAA Aerospace
Sciences Meeting and Exhibit, Reno, NV, paper no. AIAA-2006-500
5
Ashraf, M. A., Young, J., and Lai, J. C. S. 2009. Eect of Airfoil Thickness, Camber and Reynolds Number on Plunging
Airfoil Propulsion, 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition,
Orlando, FL, paper no. AIAA-2009-1274
6
Chan, W. M. 2009. Overset Grid Technology Development at NASA Ames Research Center, Computers and Fluids,
38(3): 496-503.
7
Meakin, R. L. 2001. Automatic O-body Grid Generation for Domains of Arbitrary Size, AIAA 15th Computational
Fluid Dynamics Conference, Anaheim, CA, paper no. AIAA-2001-2536.
8
Nichols, R. H., Tramel, R. W., and Buning, P. G. 2006. Solver and Turbulence Model Upgrades to OVERFLOW 2
for Unsteady and High-speed Applications, 24th AIAA Applied Aerodynamics Conference, San Francisco, CA, paper no.
AIAA-2006-2824.
9
Tuncer, I. H., and Kaya, M. 2005. Optimization of Flapping Airfoils for Maximum Thrust and Propulsive Eciency,
AIAA Journal, 43(11): 2329-2336
10
Oyama, A., Okabe, Y., Fujii, K., and Shimoyama, K. 2007. A Study on Flapping Motion for MAV Design Using Design
Exploration, AIAA 2007 Aerospace Sciences Conference, Rohnert Park, CA, paper no. AIAA-2007-2878.
11
Clark, J.P., and Grover, E.A. 2006. Assessing Convergence in Predictions of Periodic-Unsteady Flowelds,ASME
TURBO EXPO 2006: Power for Land, Sea and Air, Barcelona, Spain, 1-11
12
Klir, G. J. and Yuan, B., 1995. Fuzzy Sets and Fuzzy Logic: Theory and Applications, Prentice Hall PTR, Upper
Saddle River, NJ.
17 of 17
American Institute of Aeronautics and Astronautics