ABSTRACT
Title of dissertation: THE FLUID DYNAMICS OF
MAYFLY NAIADS
Khaled Abdelaziz, Doctor of Philosophy, 2014
Dissertation directed by: Professor Kenneth Kiger
The Mechanical Engineering Department
During their aquatic phase, mayflies Centroptilum triangulifer use a series of
tracheal gills to facilitate gas exchange. Recent experimental studies on nymphal
mayflies have identified two coupled features associated with the ontogenic progres-
sion of their ventilatory kinematics: 1) there is an abrupt shift from a rowing mech-
anism in small instars to a flapping mechanism in larger instars, and 2) the flapping
mechanism is associated with the development of a flexural hinge that permits the
passive movement of a distal flap. The primary role of the tracheal gills is tied to
ventilation rather than locomotion. As such, it is not yet understood why such a tran-
sition happens and which performance metric is improved, if any. Hence, the goal of
the current research is to investigate both features using numerical simulations.
First, a computational model of the mayfly is built from a dissected animal.
Then, a 3-level prescribed kinematic chain is introduced to a previously in-house
developed and in-house validated explicit parallel Navier-Stokes solver where both
the advective and diffusive terms are advanced explicitly using a third-order, low-
storage, Runge-Kutta scheme. Finally, an immersed boundary method based on a
moving least squares reconstruction is implemented to enforce the correct moving
boundary conditions.
Two different parametric spaces are constructed. The first one investigates the
transition from rowing to flapping kinematics, the morphological effects on the flow
field and on the proposed performance parameters while the second aims to provide
an explanation for the hinge development.
Two metrics based on control volume analysis are proposed to quantify the per-
formance of each numerical case. The first metric is simply the mechanical efficiency
of the energy transfer from the moving gills to the surrounding flow field while the
second incorporates the mass flow rate across the control surface. The second metric
is promising because it is able to provide a plausible explanation of both features by
showing that the rate of work done by the mayfly is diminished throughout ontogeny
with respect to the induced mass flow rate.
THE FLUID DYNAMICS OF
MAYFLY NAIADS
by
Khaled Abdelaziz
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2014
Advisory Committee:
Professor Kenneth Kiger, Chair
Professor Elias Balaras
Professor Peter Bernard
Professor James Duncan
Professor Jeffrey Shultz, Dean’s Representative
c© Copyright by
Khaled Abdelaziz
2014
Dedication
To GOD, the Watchful.
To Mankind, in peace.
To my deceased father, from whom I learned my first key strokes.
To my mother, for going beyond limits in raising her children.
To my brother, for his support from the moment I stepped foot in JFK.
To my sister, for her continuous encouragement.
To Mariam and Jannah, my adorable nieces.
To Ryaan, my nephew and the newest family member.
ii
Acknowledgments
I would like to express my sincere appreciation to my advisors, Professor Elias
Balaras and Professor Kenneth Kiger, for their continuous support during the course
of my graduate studies at the University of Maryland. Professor Balaras provided me
with invaluable scientific expertise and gave me the chance to work with his CFD team
and with state-of-the-art computational resources. My deepest gratitude to Profes-
sor Kenneth Kiger for his never-ending motivation, insight, inspiration, constructive
feedback, and for effectively and enthusiastically teaching me courses in viscous flows
and experimental techniques.
I take this opportunity to express my gratitude to Professors James Duncan,
Peter Bernard and Jeffrey Shultz for agreeing to serve on my PhD dissertation com-
mittee and for providing helpful feedback on my research proposal. Professors James
Duncan and Peter Bernard have also inspired me throughout many fluid mechanics
courses. I also wish to acknowledge the help provided by Professor Jeffrey Shultz for
supplying the high-resolution images and segments of the mayfly naiad.
I would like to thank my lab mates Marcos Vanella and Nikolaos Beratlis, for
their help and support at the beginning of my graduate studies and for providing me
with useful numerical codes. I also thank my fellow lab mates Grigorios Panagakos,
Don Daniel, and Clarence Baney for their helpful discussions and moral support.
Last but not least, I would like to acknowledge the technical and prompt help
of Melvin Fields our IT coordinator. I also thank Dylan Hazelwood and Dan Wysling
for their IT support.
iii
Table of Contents
List of Tables vi
List of Figures vii
List of Abbreviations ix
1 Introduction 1
1.1 What is a Mayfly? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Literature Review and Prior Work . . . . . . . . . . . . . . . . . . . . 5
1.4.1 Rowing and Flapping . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2 The Metachronal Wave . . . . . . . . . . . . . . . . . . . . . . 14
1.4.3 Advances in Mayfly Larvae . . . . . . . . . . . . . . . . . . . . 18
1.5 Research Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2 Mathematical Modeling 21
2.1 Surface Modeling of the Mayfly Nymph . . . . . . . . . . . . . . . . . 21
2.2 The Prescribed Kinematics . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.1 The Kinematic Chain . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.2 The Flapping Kinematics . . . . . . . . . . . . . . . . . . . . 41
2.2.3 The Rowing Kinematics . . . . . . . . . . . . . . . . . . . . . 42
2.3 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.4 Performance Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4.1 Mechanical Efficiency . . . . . . . . . . . . . . . . . . . . . . . 51
2.4.2 Specific Mass Flow Rate . . . . . . . . . . . . . . . . . . . . . 52
3 Numerical Methodologies and Validation 54
3.1 Spatial and Temporal Discretization . . . . . . . . . . . . . . . . . . . 54
3.2 Treatment of Immersed Boundaries . . . . . . . . . . . . . . . . . . . 57
3.2.1 MLS Reconstruction . . . . . . . . . . . . . . . . . . . . . . . 59
3.2.2 Evaluation of Surface Forces . . . . . . . . . . . . . . . . . . . 62
3.3 Computational Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
iv
3.4 Comparison with Experiments . . . . . . . . . . . . . . . . . . . . . . 66
3.4.1 Comparison with in-vivo Experiments . . . . . . . . . . . . . 67
3.4.2 Comparison with Robotic Experiments . . . . . . . . . . . . . 80
3.5 Parametric Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4 Results I: The Rowing-Flapping Transition 88
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.2 General Description of the Time-averaged Flow Fields . . . . . . . . . 89
4.3 Dynamics of Rowing and Flapping . . . . . . . . . . . . . . . . . . . 92
4.3.1 Flapping Dynamics at Reω = 21.6 . . . . . . . . . . . . . . . . 93
4.3.2 Rowing Dynamics at Reω = 21.6 . . . . . . . . . . . . . . . . 101
4.4 Viscous and Inertial Effects . . . . . . . . . . . . . . . . . . . . . . . 104
4.5 Planform Morphology . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.6 Micro-pumping and Performance Metrics . . . . . . . . . . . . . . . . 122
4.6.1 Mechanical Efficiency . . . . . . . . . . . . . . . . . . . . . . . 125
4.6.2 Specific Mass Flow Rate . . . . . . . . . . . . . . . . . . . . . 132
5 Results II: Flexibility and the Role of the Hinge 150
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5.2 Qualitative Investigation . . . . . . . . . . . . . . . . . . . . . . . . . 152
5.3 Quantitative Assessment . . . . . . . . . . . . . . . . . . . . . . . . . 159
5.4 Last Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
6 Conclusions and Future Work Suggestions 170
A Published Article: Experiments in Fluids 175
Bibliography 188
v
List of Tables
3.1 Parametric space used to explore the rowing-flapping transition . . . 86
4.1 Mass flow rates for flapping at Reω = 21.6 and rowing at Reω = 2.3 . 91
4.2 Normalized specific mass flow rate for control volume 2 . . . . . . . . 138
vi
List of Figures
1.1 Two different phases of mayfly C. triangulifer. . . . . . . . . . . . . . 2
1.2 Typical motion of cilium. Photomicrograph of a live copepod . . . . . 15
2.1 High-resolution raster images of C. triangulifer and cross sections . . 23
2.2 Computational model of C. triangulifer . . . . . . . . . . . . . . . . . 24
2.3 Gill morphology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Illustration of the kinematic chain . . . . . . . . . . . . . . . . . . . . 31
2.5 Illustration of the forward transformation between levels 1 and 2 . . . 35
2.6 Schematic of position vectors in different levels . . . . . . . . . . . . . 38
2.7 Prescribed mayfly flapping kinematics for gill 4 . . . . . . . . . . . . 43
2.8 Eulerian angles for mayfly flapping kinematics . . . . . . . . . . . . . 44
2.9 Prescribed mayfly rowing kinematics for gill 4 . . . . . . . . . . . . . 45
2.10 Eulerian angles for mayfly rowing kinematics . . . . . . . . . . . . . . 46
2.11 Schematic of control volume . . . . . . . . . . . . . . . . . . . . . . . 50
3.1 Schematic of computational Domain . . . . . . . . . . . . . . . . . . 66
3.2 Comparison of the velocity magnitude for flapping at Reω = 21.6 . . . 69
3.3 Comparison of u for flapping at Reω = 21.6 . . . . . . . . . . . . . . . 71
3.4 Comparison of v for flapping at Reω = 21.6 . . . . . . . . . . . . . . . 72
3.5 Comparison of ωz for flapping at Reω = 21.6 and t/τ = 0.102 . . . . . 74
3.6 Comparison of ωz for flapping at Reω = 21.6 and t/τ = 0.326 . . . . . 75
3.7 Comparison of the velocity magnitude for rowing at Reω = 2.3 . . . . 79
3.8 Comparison of the instantaneous fields for rowing at Reω = 2.3 . . . . 81
3.9 Illustration of computational setup and slice at z/L = 1.3 . . . . . . . 82
3.10 Comparison with robotic experiments for rowing at Reω = 17.6 . . . . 84
4.1 Time-averaged flow fields for rowing and flapping . . . . . . . . . . . 90
4.2 Instantaneous flapping dynamics at Reω = 21.6 . . . . . . . . . . . . 95
4.3 Flapping dynamics at Reω = 21.6 and t/τ = 0.65 . . . . . . . . . . . 97
4.4 Pathlines for flapping kinematics at Reω = 21.6 . . . . . . . . . . . . 100
4.5 Instantaneous rowing dynamics at Reω = 21.6 . . . . . . . . . . . . . 102
4.6 Effect of Reω on ωz for flapping kinematics . . . . . . . . . . . . . . . 107
4.7 Effect of Reω on ωz for rowing kinematics . . . . . . . . . . . . . . . . 108
vii
4.8 Pathlines for flapping kinematics at Reω = 21.6 and Reω = 1.0 . . . . 111
4.9 Pathlines for rowing kinematics at Reω = 21.6 and Reω = 1.0 . . . . . 113
4.10 Effects of planform morphology on flow fields at Reω = 21.6 . . . . . 118
4.11 Effects of planform morphology on flow fields at Reω = 21.6 . . . . . 119
4.12 Mechanical efficiency, η, for 3 control volumes . . . . . . . . . . . . . 127
4.13 Viscous and pressure work rate on the ventral plane . . . . . . . . . . 131
4.14 Specific mass flow rate, ξ, for CV1 . . . . . . . . . . . . . . . . . . . . 135
4.15 Specific mass flow rate, ξ, for CV1 (normalized axis) . . . . . . . . . 139
4.16 Sensitivity of specific mass flow rate, ξ, to the control volume size . . 141
4.17 Quantitative mass flow rate analysis around the mayfly . . . . . . . . 143
4.18 Quantitative mass flow rate analysis around the mayfly (%) . . . . . 144
5.1 Time-averaged vorticity component ωz for Reω = 21.6 . . . . . . . . . 153
5.2 Time-averaged normalized velocity magnitude for Reω = 21.6 . . . . . 154
5.3 Pathlines integration for gill pair 4 at Reω = 21.6 . . . . . . . . . . . 157
5.4 Specific mass flow rate, ξ, for the parametric space (normalized axis) 160
5.5 Rate of work done for flapping using gill pair 4 at Reω = 21.6 . . . . 163
5.6 Viscous and pressure work rate on left gill 4 . . . . . . . . . . . . . . 165
viii
List of Abbreviations
CV : Control Volume
CS : Control Surface
FFT : Fast Fourier Transform
DFFT : Discrete Fast Fourier Transform
FOR : Frame of Reference
IB : Immersed Boundary
LHS : Left Hand Side
MLS : Moving Least Squares
MPI : Message Passing Interface
NSE : Navier-Stokes Equations
PDE : Partial Differential Equation
RHS : Right Hand Side
RK3 : Runge-Kutta 3rd order
SP : Stroke Plane
wrt : with respect to
+ve : Positive
-ve : Negative
Anterior : Of or near the head end of the mayfly
Coronal : Divides the body into ventral and dorsal sections
Dorsal : Belonging to or near the upper surface of the mayfly
Medial : Dividing the mayfly into right and left halves
Posterior : Located at or towards the rear of the mayfly
Sagittal : Located in a plane that is parallel to the medial plane
Transverse : Divides the body into superior and inferior parts
Ventral : Belonging to or near the belly surface of the mayfly
C. triangulifer : Centroptilum triangulifer
c : cosine function
s : sine function
A : Peak-to-peak amplitude
L : Gill length and reference length in NSE
f : Frequency
NrM : Position vector from frame N to point M
t? : Physical time
t : Dimensionless time
U : Reference velocity in NSE
u : Dimensionless velocity vector
ix
x : Dimensionless Cartesian position vector
XαYZ : angular acceleration of frame Y relative to frame X in Z frame
β : Stroke plane inclination angle
γ : Hinge flexion angle
 : Stroke plane lateral offset angle
η : Mechanical efficiency
θ : Stroke plane deviation angle
λ : Proximal-distal inclination angle
µ : Dynamic viscosity
ν : Kinematic viscosity
ξ : Specific mass flow rate
ρ : Fluid mass density
τ : Beating Period
φ : Angle that modifies the pitch
ψ : In-plane beating angle
XωYZ : angular velocity of frame Y relative to frame X in Z frame
Reω = fL2/ν : Oscillating Reynolds number
Ref = fAL/ν : Flapping Reynolds number
ReL = ULbody/ν : Translation Reynolds number
St = fL/U : Strouhal number
P e´ = LU/D : Pe´clet Number
x
Chapter 1: Introduction
1.1 What is a Mayfly?
Mayflies are classified in the insect order Ephemeroptera, meaning in Greek
“lasting but one day”. The name refers to the fact that the winged adults live only
one to two days. Mayfly adults are winged, terrestrial and are also referred to as
imagoes (sexually mature). Mayflies are unique in having a winged subimago stage
(sexually immature). Early in their life -referred to as larvae, naiads or nymphs-,
mayflies are aquatic (Fig. 1.1). Most of a mayfly’s life cycle is spent as a naiad.
Naiads are found in wetlands, ponds, lakes, slow or swift flowing streams and rivers.
There are many species of mayflies, and they are very common and diverse in river
habitats. Many species are relatively intolerant of pollution. The chief importance
of mayflies in the ecosystem lies in their value in breaking down organic matter
and as food for fish, although they are also food for many other animals, including
amphibians, birds, spiders and insects.
Mayfly naiads possess gills, which can be divided into two types: the plate-like
abdominal gills and the accessory gills. From a biological standpoint, the abdominal
gills have been studied in detail, and a variety of functions have been documented,
among which are respiration, osmoregulation, locomotion, water circulation, protec-
1
(a) (b)
Figure 1.1: Two different phases of mayfly C. triangulifer. (a) Aquatic phase, photo-
graph courtesy of David Funk. (b) Terrestrial phase.
tion and maintenance of the position of the body [84]. Accessory gills have other
characteristics that distinguishes them from the abdominal gills. For example, they
are thread-like or finger-like filaments or tufts, never plate-like. Also, they are sit-
uated on the head and/or thorax, never on the abdomen. Most genera of Baetidae
(small minnow mayflies) lack these accessory gills and the species under considera-
tion in this research (Centroptilum triangulifer) is not any different. According to
the Integrated Taxonomic Information System (www.itis.gov), the taxonomic hier-
archy of the species under consideration is: Animalia (Kingdom); Arthropoda or
jointed-foot (Phylum); Hexapoda (Subphylum); Insecta (Class); Pterygota or winged
insects (Subclass); Palaeoptera or ancient winged insects (Infraclass); Ephemeroptera
or mayflies (Order); Pisciforma (Suborder); Baetidae (Family); Centroptilum (Eaton
1869) (Genus), Centroptilum triangulifer (McDunnough, 1931) (Species).
2
1.2 Motivation
A typical feature of the mayfly naiad is the presence of 7 pairs of abdominal
gill plates on the lateral dorsal region of the abdomen. In many species, these plates
actively beat in a metachronal fashion to produce an external current that aids in the
circulation of fresh oxygenated water and allows mayflies to tolerate fluctuations in
oxygen concentrations [3]. More recent studies documented the ventilatory kinemat-
ics of the nymphal mayfly C. triangulifer [57–59]. In these studies, two distinctive
and coupled features associated with the ontogenic progression of the ventilatory
kinematics of this specific species (C. triangulifer) have been identified: 1) there is
an abrupt shift from a rowing mechanism in small instars to a flapping mechanism in
larger instars, and 2) the flapping mechanism exhibited the development of a flexural
hinge that allowed for the passive movement of a distal flap. The reported oscillating
Reω ranged from 2.0 to 22.0. This was the first study to reveal that a switch from
a rowing mechanism to a flapping mechanism can occur systematically in the same
species as it molts or grows and therefore highlighting the effect of increasing the Re,
and possibly morphology, on beating kinematics.
Centroptilum triangulifer therefore presents a very interesting problem for both
areas of fluid mechanics and biology. For the former, C. triangulifer operates in
the transitional regime where a demarcation from a viscosity-dominated flow to an
inertia-dominated flow is known to occur. Examining the behavior of C. triangulifer
may shed more light on how important gait changes are, in order to escape viscosity,
increase mechanical efficiency or thrust, or enhance water circulation. For the latter,
3
C. triangulifer traverses the biologically intermediate regime (1 < ReL < 100) where
it is known that below this regime (ReL < 1), rowing is exclusively used by the
biological taxa while above this regime (ReL > 100), flapping is predominantly used.
Indeed much research has been conducted that aims to link rowing mechanisms to
viscous flows and flapping mechanisms to inertial flows, but it is unwarranted to
generalize previous hypotheses especially that rowing is often used at ReL > 100 and,
in some cases, is associated with better maneuverability, higher thrust or intermittent
burst-and-coast. Mayfly larvae might be optimizing for respiration rather than thrust
or mechanical efficiency and therefore presents a different type of a problem.
What makes the problem even more appealing is its applied engineering perspec-
tive, where understanding the micro-circulation dynamics around mayfly nymphs can
help improve future generations of distributed and autonomous micro-scale chemical
sensor networks, which often require an external convective current to achieve rapid
and reliable performance. Such requirement was concluded by a recent report by the
National Academy of Sciences [71]. Micro-scale sensors are usually limited by diffu-
sion which calls for the use of a convective sampling system to enhance performance,
and indeed, use of such systems has shown over an order of magnitude decrease in
response time when implemented [48, 56, 63, 70]. Bio-inspired sensors have exciting
applications in micro-fluidic devices, which was recently demonstrated by fabricated
artificial cilia consisting of electro-statically actuated polymer structures. The device
showed substantial fluid flow experimentally [13], which was also modeled numeri-
cally [5]. Worth mentioning is that conventional centrifugal pumps, for example, are
typically ineffective at low Re and have only been miniaturized to a limited extent [45],
4
which paves the road towards unconventional designs for micro-pumping.
1.3 Objectives
The objective of the current work is to determine why the transition from a
rowing-like mechanism to a flapping-like mechanism occurs in C. triangulifer, to pos-
sibly corroborate previous experimental results [57–59] and present new ideas or hy-
potheses. The research also aims to shed more light on the more general question as
to why rowing is exclusively used as one approaches the viscosity-dominated flows.
Using solely numerical techniques, the goal is to simulate hypothetical scenarios that
were impossible to obtain using in-vivo experiments such as rowing at high Reω or
flapping at low Reω or flapping without a flexural hinge, and interpret the numerical
results in order to explain the systematic behavior of growing nymphs. These hypo-
thetical scenarios could be obtained using robotic mayflies, but given the complicated
nature of the kinematics, the metachronal wave, and the three-dimensionality of the
flow field, this approach presents a different set of challenges.
1.4 Literature Review and Prior Work
This section provides a comprehensive overview of prior work related to my
research, presents the problem in greater depth and familiarizes the reader with the
different aspects related to it as, for example, the highly viscous nature of creeping
flows. The multidisciplinary problem is tied to what could be optionally classified into
three different categories of research. The first research area is basically the rowing
5
and flapping axis in the biological taxa, which has been under investigation for many
decades. The rowing-flapping axis is largely, if not entirely, an aquatic phenomenon
and is rarely discussed in the flight literature [74] since most flight literature deals
with flapping flight. Rowing and flapping, however, in the context of enhancing micro-
mixing or micro-pumping efficiency is almost nonexistent and therefore is the primary
contribution from my work. Second, a need to assess the performance of the micro-
mixing of the mayflies rose, and given the possible applicability towards bio-inspired
and autonomous sensors, the second research area is basically articles related to micro-
scale sensors and lab-on-a-chip devices aiming to look primarily into how these sensors
are being evaluated (e.g. response time, effectiveness, energy consumption, etc.).
Both these areas have been synthesized in order to present a coherent description
of the rowing-flapping dichotomy and the metachronal wave in sections 1.4.1 and
1.4.2 respectively. The last research area pertains to the numerical methodologies
used to perform the simulations, which relies primarily on the immersed boundary
technique [4, 69], which will be presented in Chapter §3.
The following subsections are divided as follows. First, prior work on rowing and
flapping is reviewed in §1.4.1. Second, a description of the metachronal wave and how
it may affect the flow field is presented in §1.4.2. Third and last, advances pertaining
to mayfly larvae specifically is presented, and its relation to rowing, flapping and the
metachronal wave is discussed in §1.4.3.
6
1.4.1 Rowing and Flapping
The primary distinction between rowing and flapping is made by comparing the
direction of the net force generation (thrust) with respect to the stroke plane. Flap-
ping strokes are characterized by an appendage motion perpendicular to the direction
of the net force generation, while rowing strokes are characterized by predominant
motion in the direction of net force generation. With the primary objective of the
mayfly larva under consideration being different than locomotion, the difference be-
tween rowing and flapping is tweaked by simply comparing the direction of the mean
flow (instead of the thrust) with the motion of the gills. If the mean flow is perpendic-
ular to the motion of the gills, then the beating kinematics are classified as a flapping
mechanism. On the contrary, if the direction of the mean flow is almost parallel to
the motion of the gills, then the kinematics are classified as a rowing mechanism. In
order to remove any confusion that may arise from the distinction between rowing
and flapping, it is emphasized that rowing and flapping are idealized extremes of a
continuum [74]. In other words, there exist many sets of kinematics that produce
mean flow (or thrust) in directions that are not necessarily parallel or perpendicu-
lar to the appendages motion as for example in mayfly naiads in their intermediate
regime [59]. Flapping strokes are usually symmetric while rowing strokes introduce
some asymmetry. The asymmetry in rowing can be classified into temporal and spa-
tial asymmetry. Temporal asymmetry is introduced by using a power stroke and a
recovery stroke where the recovery stroke is slower in order to reduce drag. Temporal
asymmetry is known to fail in the Stokesian regime (scallop theorem [53]). Spatial
7
asymmetry is introduced by changing the shape or profile of the appendage during
the recovery stroke in an attempt to reduce drag. This is accomplished through ar-
ticulation of the appendage (remipede), actuated flexion (cilia in the ctenophore) or
retraction of setae on the limb (Artemia, remipede) [57]. Spatial symmetry may also
be broken by the surrounding fluid domain as for example when beating near a wall
(see for example [46]).
The subject of rowing and flapping is rarely discussed in the flight literature
and most research on insect flight has been done on flapping flight. There has been
a mention that the smallest flying insects swim, or row, in air but then the geometry
of rowing would be quite different from that occurring in aquatic flight [74]. In his
review of unsteady mechanisms of aquatic and aerial locomotion (intermediate and
high Reynolds number), Dickinson [14] pointed out that the separate histories of the
two fields (aerial and aquatic locomotion) serve as an obstruction to a useful and
powerful synthesis of the knowledge gained throughout the previous decades. For ex-
ample the concept of added mass is variably called ’added mass’, ’added mass inertia‘,
’virtual mass’, ’acceleration reaction’ within the biological literature. Nevertheless,
the abundant literature about insect flapping flight provides profound interpretations
to the inertia-dominated regime and something to look for in the transitional regime
encountered in this research. One important lesson to learn from insect flight litera-
ture is that non steady-state (instead of quasi-steady) approaches are currently much
more common in analysis, due in large part to the influence of Torkel Weis-Fogh in
the 1970s [76] and more recently of Charlie Ellington and his colleagues [21–25]. This
is due to the fact that quasi-steady analysis failed to explain many insect flights, for
8
example by frequently underestimating the lift coefficient. A numerical solution of
the full unsteady Navier-Stokes equations is the best approach given the availability
of higher computing power.
Another nomenclature is repeatedly associated with rowing and flapping kine-
matics, which is the “drag-based” and “lift-based“ propulsion respectively. I speculate
that the terms loosely originated to indicate or imply that the low Reynolds num-
ber regime is dominated by viscous drag and therefore a “drag-based” propulsion or
mechanism, while the high Reynolds number regime is dominated by inertia forces
producing lift and therefore the term “lift-based” propulsion or mechanism. This
nomenclature may be misleading. After all, how can a drag-based mechanism pro-
duce lift (needed at least to sustain weight in non-neutrally buoyant animals or flying
insects) if by formal definition drag and lift are orthogonal to each other? On the
other hand, at high Reynolds numbers, part of the pressure forces contributes to the
drag, known as form drag which is clearly not implied by the term “lift-based”. At
best, the nomenclature may be a good approximation in a heuristic model of heav-
ing and pitching plates but numerous mechanisms, including viscous stresses, bound
circulation, flow separation, fin-vortex interactions, fin inertia, added mass (accelera-
tion reaction) forces, and jet-like (squeeze) forces, potentially contribute to the force
balance in both rowing and flapping foils that makes it unwise to distinguish beating
or locomotion mechanisms using only these two terms. Again, the “drag-based” and
“lift-based” distinction seem to work only for idealized extremes. As Walker put it:
“The drag vs. lift dichotomy, then, while useful as an elementary model, both con-
founds and obscures the different sources of forces on real oscillating fins.” [74], which
9
was emphasized earlier by Dickinson [14]. A detailed presentation of such unsteady
forces associated with inertia-dominated flows, may be found in [11,55,75].
In one of the more recent and rich technical articles that examines the rowing-
flapping dichotomy, Walker presented a comprehensive comparative study of rowing
and flapping [72]. Based on the available data, Walker showed that rowing occurs
throughout the range of biologically intermediate Re (1 < Re < 100). On the other
hand, flapping is restricted to Re > 10. In addition, for neutrally buoyant animals,
the data suggested that flapping is restricted to Re > 100. Walker pointed out that
the minimum Re for flapping appendages should be accepted cautiously because of the
scarcity of the data on the most relevant animals. In other words, the critical Reynolds
number shall not be extended to all biological entities in general and certainly not
to C. triangulifer. Using the blade-element method and despite the simplicity of his
model, Walker concluded that the results predict (or retrodict) both the exclusive use
of a rowing stroke below a cut-off Re of about 20 and the reliance on a reduced area
and span during the recovery stroke of small rowing animals.
The problem has been also approached by Childress and Dudley [12] using ex-
periments and several simplified theoretical models. They performed experiments on
the swimming of the shell-less pteropod mollusc Clione antarctica, which achieves
locomotion using ciliated surfaces alone (rowing) or using a pair of flapping wings.
It was proposed that forward, reciprocal flapping flight is impossible below some fi-
nite critical value of the Reynolds number provided the motion of the insect be fully
determined by that dimensionless number. They substantiated their findings using
theoretical approaches, namely a modified Oseenlet model and a venetian blind flap-
10
ping model. The former showed a critical Reynolds number of 33 for an approximate
two-dimensional wing model, while the latter showed it to be 9.23 for a spatial period
of 1.
In yet another simple and elegant experimental approach, and employing a
rotational horizontal flat rectangular wing, Vandenberghe [66, 67] was able to show
that flapping flight occurs as a symmetry-breaking bifurcation from a pure flapping
state with no prior horizontal motion. A sharp bifurcation to flapping flight for a
frictionless system in the range 20 < Ref < 55 was indicated. Visualization of the
flow field around the heaving and plunging foil showed a symmetric pattern below
transition, while above threshold, an inverted von Ka´rma´n vortex street was observed
in the wake of the wing.
Another numerical study shows that flapping has a higher mechanical efficiency
for all biologically relevant swimming speeds and therefore it was suggested that flap-
ping is usually employed by animals that require energy conservation such as during
migration (e.g. adult green turtles and loggerheads), while rowing maximizes maneu-
verability such as acceleration, turning and braking such as during escape swimming
or feeding [73]. The existence of both flapping and rowing in the intermediate regime
and the occasional rowing in the inertial regime might be explained by the fact that
the optimized parameter or parameters are not known. After all, locomotion is not
only about getting efficiently from one point to another, but also a complex and
flexible system of behavior required for feeding, courtship, escape [14] and migration.
[1] used in-silico experiments and in particular the vorticity stream function to
study a 2D flapping ellipsoidal body. By introducing a perturbation normal to the
11
flapping direction, they showed that at sufficiently large Ref , unidirectional locomo-
tion emerges as an attracting state. The onset of the instability (steady locomotion)
depended on the aspect ratio of the body and appeared to be independent of the
density ratio. They also showed that for sufficiently small Ref , horizontal motions
imposed transiently on the flapping body rapidly dissipate and the fluid resumes a
qualitative left-right symmetry.
Now that it is established that a transition to rowing at low Reynolds numbers
must exist [12,67,72,73], the question becomes why is flapping flight “prohibited” at
low Reynolds numbers? To answer that, one must first look into how really “viscous”
a fluid can get as one approaches the Stokes regime and how bad it is for living
organisms. First, an organism is confronted with a large viscous drag and its inertia
plays no role. To give an example, an organism with one micron order that is moving
with a typical biological speed of 30µm/s will coast for 0.1A˚ in 0.6µs to a complete
stop. This is a tremendous deceleration of 450m/s2. In other words, the organism
stops instantaneously once the external force is removed (the distance coasted is 5
orders of magnitude less than that of the organism size). Simply speaking if an
organism needs to move in a Stokesian realm, it has to keep propelling itself all
the time. Coasting is just not possible. Viscous flows are not easily intuitively
grasped simply because they are not encountered on a daily basis, as opposed to
other inertia-dominated flows (swimming, driving a car, insect flight, water faucet,
etc.). A second characteristic of the Stokesian realm is the low mechanical efficiency
associated with moving organisms. For example, a sphere driven by a helical propeller
in a creeping flow has a mechanical efficiency of only 1%. Nevertheless, it was pointed
12
out that the energy exerted in locomotion still represents a minute amount of an
organism’s metabolism and that efficiency may not be of primary concern [53]. The
third important aspect of creeping flows is the kinematical reversibility resulting from
the linearity of its basic differential equations. In other words, any creeping flow
solution is indistinguishable from its solution if the flow direction is reversed [77]. This
paves the road to why symmetrical flapping would not work in the Stokesian realm and
was further solidified by the famous scallop theorem: A neutrally buoyant organism
exhibiting time-reversal symmetry is a non-swimmer in the Stokesian realm [11]. The
scallop theorem could be traced back to the famous lecture by Purcell [53] where it
was also shown that with only one degree of freedom in configuration space, one can
only have reciprocal motion and therefore swimming is not possible at low Re.
So how can an organism propel itself in the Stokesian realm? One way is to
introduce some asymmetry in the oscillating appendage other than time asymmetry
(scallop theorem). Such asymmetry can be accomplished by changing the appendage
shape or profile during the recovery stroke, which, in turn, might help reduce the
viscous drag. This mechanism is sometimes referred to as non-reciprocal drag-based
propulsion. One abundant example is cilia (Figure 1.2a). A cilium (Latin for eyelash)
has a diameter of order 0.2µm and a typical length that varies from 2µm to 15µm.
Cilia are found in humans (lining of the trachea to sweep mucus, Fallopian tubes to
move ovum). They are also found in a diversity of aquatic animals like ctenophores
(comb jellies) [6] and molluscs [12]. To achieve asymmetry, a cilium moves approxi-
mately as a straight rod during the effective stroke, while during the recovery stroke,
it rolls close to the surface in a tangential motion. The mechanical work done of
13
cilia during the effective stroke was modeled computationally and was shown to be
approximately five times the amount of work done during the recovery stroke [31].
Another way to introduce spatial asymmetry is achieved by articulation of the
appendages as in remipede (Remipedia: Latin for oar-footed) [40], or retraction of
setae (Latin for bristle) on the limb as in Artemia (brine shrimp) [78, 79] and remi-
pede. The most recent study on C. triangulifer showed that at its low Re end, the
gills are fairly stiff offering little chance to change its shape. Therefore, the spatial
asymmetry is introduced by changing the gills pitch during the recovery stroke in
order to “feather“ them into the incoming flow and minimize the drag [57].
A recent article has classified asymmetry into temporal, spatial and orienta-
tional asymmetry [38], where temporal asymmetry means that one of the half strokes
lasts longer than the other, spatial asymmetry is when the shape of the appendage
changes between power and recovery strokes (e.g. cilia), and orientational asymme-
try is the asymmetry introduced with respect to the channel boundaries for example.
The temporal asymmetry is known to be ineffective in the Stokesian regime (scallop
theorem) since both half strokes occur in a highly viscous regime. A system achieving
“configurational” symmetry (the three aforementioned symmetries) will not generate
a net flow or thrust in the Stokesian realm [38].
1.4.2 The Metachronal Wave
An animal having multiple appendages has one more degree of freedom, which is
the phase lag between its beating appendages. On some occasions, appendages beat
14
(a) (b)
Figure 1.2: (a) Typical asymmetric motion of a cilium. The dashed lines represent
the trajectory of the tip of an individual cilium [39]. (b) Photomicrograph of a live
copepod, Euchaeta norvegica [83].
in phase, but in most scenarios, there is a phase lag between the beating appendages
that introduces a spatial wave, which travels along the animal in what is known
among biologists as a metachronal wave. A metachronal wave that travels in the
same direction as the power stroke is called a symplectic metachronal wave. On the
contrary, a metachronal wave that travels in the opposite direction of the effective
stroke is termed an antiplectic metachronal wave. An ametachronal wave is different
from the typical metachronal wave in that many consecutive legs beat in the same
phase. Many research articles investigated the benefits of metachronal waves using
either a numerical approach or an experimental one.
On the modeling side, Gueron modeled multiciliary configurations. It was found
that when two adjacent cilia start beating at different phases, they became synchro-
nized within several beating periods which was also observed in experiments when
two flagella are brought into close proximity. The models also showed that an ap-
15
proximately antiplectic wave pattern evolved autonomously and it was conjectured
that metachronism in cilia may occur, at least partially, as an organized phenomenon
due to the hydrodynamic interactions between neighboring cilia [30, 32]. In other
words, the metachronal wave occurs as a result of the fluid-structure interaction and
the proximity of the cilia to each other. Later it was shown that having neighboring
cilia beat metachronally is energetically advantageous [31]. Using a coupled magneto-
mechanical solid-fluid computational model, Khaderi et al. [39] showed that artificial
cilia beating metachronally enhance the fluid flow and that antiplectic metachrony
leads to a considerable flow enhancement compared to symplectic metachrony, but
only when the cilia spacing is small.
On the experimental side, the flow of water around the comb plates of Pleuro-
brachia pileus (phylum Ctenophora, class Tentaculata (comb jellies)) was studied [6].
Typical ctenophores swim by means of eight rows of comb plates where each comb
plate resembles a rectangular paddle and is composed of a hundred thousand or more
cilia. A typical power stroke of the comb plate has a Reynolds number of 9. Using sus-
pended plastic beads and digitized high-speed films, it was shown that the antiplectic
motion of the comb plates smooths out the otherwise intermittent flow enabling the
jelly to impart greater momentum to the flow. Otherwise said, the comb plates work
together by capturing particles accelerated by the preceding comb plate. There was
no quantitative approach to compute the efficiency, but qualitative investigation of
the flow field suggested that antiplectic metachronism may be responsible for the pro-
gressive acceleration of water to a much higher speed than would be possible with a
single paddle thereby increasing the mechanical efficiency. It is interesting to mention
16
that the beating frequency of Pleurobrachia pileus could be controlled externally by
a fine glass needle depressing against the shaft of the first comb plate at the aboral
(back) end.
Studies on individual species that transit the biologically intermediate regime
have shown that different types of metachronal waves may be used depending on the
life situation. For example, remipede crustaceans beat with one set of appendages
using a synchronized ametachronal stroke pattern during escape, while they beat with
another set of appendages using a metachronal stroke during cruising and swimming.
The copepod Temora longicornis forages using a metachronal rowing stroke with
the first three of its five feeding appendages while the other two row synchronously,
but rows with its swimming legs during escape with a metachronal power stroke
and synchronized recovery stroke [65]. Copepods have been studied extensively to
investigate the different types of propulsion and metachronal waves [34–37,65,83].
Many mayfly species use a metachronal wave or rhythm to beat their gills
[3,15–18,20,57–59], but only recent studies have revealed that C. triangulifer not only
uses an antiplectic metachronal wave, but also switches systematically throughout
ontogeny [57–59] from a rowing mechanism to a flapping mechanism thereby bringing
an interesting problem into light where transition from antiplectic metachronal rowing
to metachronal flapping is occurring. In the following section, previous work on mayfly
larvae is reviewed.
17
1.4.3 Advances in Mayfly Larvae
The rhythmical movements of tracheal gills could be traced all the way back to
Baba´k and Foustka [3] where they performed experiments on many arthropods in-
cluding Ephemeroptera larvae. Baba´k and Foustka showed that the different beating
frequencies were dependent on the oxygen content of the water. They observed a
frequency of 200bpm (beat per minute) in previously boiled water, 100bpm in normal
running water and as low as 3.5bpm in well-aerated water. It was also indicated that
the lower oxygen content affects the beating amplitude.
In 1932, Eastham [15] presented in detail the water currents generated by
nymphal gill movements for five different mayfly species. He speculated that tracheal
gills moving in a metachronal rhythm create differences of pressure in the inter-gill
spaces leading to a period of suction followed by a period of compression. Using
suspended mud, cultures of ciliate protozoa, and a stroboscope, he was also able
to provide crude visualizations of the currents produced by three different species
(Caenis horaria in 1934 [16], Leptophlebia marginata in 1936 [17] and Ecdyonurus
venosus in 1937 [18]). In 1938, Eastham [19] also studied a species that is quite dif-
ferent from those previously investigated [15–18], namely, the nymphs of Ephemera
danica which lives in burrows in the sandy mud of flowing waters. He showed that
gills move in metachronal rhythm that starts backwards, and set up currents, which
are symmetrical with the body axis. In 1958, Eastham made another contribution
to mayflies [20] when he presented the nymph of Chloeon dipterum L., described its
bilamellate gills, its metachronal rhythm, the currents produced by the gills over the
18
body, and its muscular mechanisms, which makes these and swimming possible. East-
ham illustrated again the periods of suction followed by periods of squeezing between
the posterior lamella of one gill and the anterior lamella of the next. He also showed
that the normal musculature is modified to provide direct muscles to gill bases and
to enable swimming using rhythmical movements of the abdomen. In his literature
about nymphal mayflies, Eastham never mentioned an abrupt switch or change in gill
kinematics or shape within the same species.
It was not until recently, however, that more quantitative information started
to emerge. Sensenig, Kiger & Shultz [58] used Particle Image Velocimetry (PIV) to
investigate the fluid dynamics of C. triangulifer. In particular, the velocity fields as
well as the gill kinematics were measured for 13 mayfly naiads with ages ranging from
20 to 42 days, and gill lengths ranging from 0.25mm to 0.77mm. The corresponding
oscillating Re ranged from 2.1 to 21.6. The results are presented in detail in [57]
and it is this quantification of the biological model system that makes up one of the
cornerstones of the upcoming computational study.
In a follow up study, Sensenig, Kiger & Shultz [59], investigated the hydrody-
namic pumping mechanism, also termed “phased vortex pump”. They hypothesized
that rowing should be favored when the radius of the vortices shed by the gills ex-
ceeds the inter-gill spacing, while flapping should be favored when the vortices radius
is smaller than the inter-gill spacing. Their results showed that the vortices radius
decreases relative to the inter-gill spacing during mayfly growth, and therefore a tran-
sition point exists during ontogeny. The vortices were defined within each PIV frame
as those spatial regions where vorticity exceeded 25% of the peak value observed
19
during the complete beating cycle.
1.5 Research Outline
The goal of the current research is to investigate the transition from a rowing
type of kinematics to a flapping type of kinematics in C. triangulifer, where the
flapping kinematics is coupled with an increased gill flexibility and the development
of a distal flap. The outline of the thesis is as follows. §2 describes the 3D modeling of
the mayfly nymph, the 3-level kinematic chain used to prescribe the gill kinematics,
the flapping and rowing kinematics in C. triangulifer, the mathematical modeling and
the performance metrics used. §3 will briefly describe the numerical methodologies
along with the immersed boundary formulation. This will be followed by comparisons
with previous experimental results. §4 will present the results for the first parametric
space that was designed to investigate the rowing-flapping transition, while §5 will
investigate in more detail the effect of the flexural hinge. Finally, §6 summarizes the
conclusions along with possible future extensions.
20
Chapter 2: Mathematical Modeling
The goal of this work is to model the laminar micro-circulation of water around
mayfly larvae using mass and momentum conservation laws. Therefore, four major
milestones were accomplished. First, an accurate mayfly 3D model was created from
a dissected animal. Second, the gill kinematics have been translated into a set of
Eulerian angles, which were implemented into the mathematical model using a 3-level
kinematic chain. The first and second milestones basically enable the introduction of
the time-varying boundary conditions into the governing equations. Therefore, the
third milestone was implementing the 3D model and the prescribed kinematics into
the mathematical modeling, where both mass and momentum conservation laws are
satisfied. The last milestone deals with the quantitative assessment of the micro-
circulation around the mayfly nymphs. Sections 1 through 4 present milestones 1
through 4 respectively. The numerical methodologies will be introduced in §3.
2.1 Surface Modeling of the Mayfly Nymph
In order to implement the mayfly nymph into the numerical solver, it was nec-
essary to build an accurate 3D model. The groundwork for the solid modeling was 4
high-resolution raster images of the mayfly nymph, which were provided by Professor
21
J. Shultz 1. These images represent the ventral (bottom), dorsal (top), lateral (side)
perspectives of the mayfly, magnified under the microscope with a magnification fac-
tor of 32, 25 and 32 respectively (see Figures 2.1b and 2.1c). A high-resolution image
was also provided for the caudal filaments 1. In addition to these, 10 posterior and
10 anterior sections were provided for 10 segments along the mayfly’s abdomen (Fig-
ure 2.1a) which were obtained by patiently and artfully dissecting one nymph under
the microscope using a drawing tube1.
The raster images and segments were then imported to a CAD software (Au-
toCAD 2008). Splines were visually fitted in order to vectorize the 20 sections and
the main body longitudinal axis. Only one half of each segment was vectorized then
mirrored by the medial plane in order to guarantee symmetry (otherwise, undesirable
asymmetry will be introduced into the numerical solver). Care has been given to
insure the smoothness of the sections at the medial plane by having tangents that are
as close as possible to being normal to the medial plane. Using the lateral view of
the nymph (Figure 2.1c), the location of each section has been identified and the 20
section have been positioned on the longitudinal axis. Then, each of the 20 sections
has been rotated so that its plane lies normal to the longitudinal axis. Five more
sections for the thorax and head have been created using the raster images making
the total number of sections 25. The end result of these CAD operations is shown in
Figure 2.2a. Finally, all the 25 sections have been input to the AutoCAD command
“loft”, which reconstructed the surface of the mayfly using interpolation techniques.
After some considerable processing time, the complicated surface of the mayfly ab-
1Professor Jeffrey Shultz, Entomology Department, University of Maryland, College Park.
22
Figure 2.1: High-resolution raster images of C. triangulifer and cross sections.
(a) Posterior and anterior margins of 10 abdomen segments of the mayfly nymph.
Sections where traced using a drawing tube attached to a high-magnification dissect-
ing microscope.(b) Dorsal view of the mayfly nymph at 25X (ruler in mm). (c) Lateral
view of the mayfly nymph at 32X (ruler in mm).
domen, thorax and head emerged and the tapered nature of each segment and the
“bulges” at the gills roots conformed to what had been seen under the microscope 1,2.
The final rendered version of the mayfly nymph is shown in Figure 2.2b. Noteworthy
is that mayfly antenna, legs and caudal filaments were not included in the model since
23
Figure 2.2: Computational model of C. triangulifer nymph. (a) 25 sections created
in modeling software. Black: posterior sections. Grey: anterior sections. Red: thorax
and head sections. Gills shown are decreasing in size towards the head and the tail
(shown for illustration only). (b) Rendered version of C. triangulifer nymph showing
different landmarks on the gills: (R) Root, (V) Ventral, (D) Dorsal and (T1) Tip.
Blue: proximal plate. Red: distal plate.
their contribution to the fluid mechanics of the problem are negligible.
In order to create the mayfly gills, the same but simpler procedure has been
followed. First, the raster images of the tracheal gills 1 are imported to AutoCAD,
then splines are fitted visually in order to vectorize the gill shape. Figure 2.3 shows
both raster and vectorized pictures of the different gill planforms used throughout
the study.
The variation in the gill planform from younger nymphs to older nymphs was
thought to have a marginal effect on the flow field. While this turned out to be gener-
ally true, its effect on the proposed performance parameters could not be overlooked.
In particular, the mean flow direction and vortices configuration were not affected by
the planform variation, but the work done by the gills as well as the induced mass
1Professor Jeffrey Shultz, Entomology Department, University of Maryland, College Park, MD.
2Professor Andrew T. Sensenig, PhD, Tabor College, Hillsboro, KS
24
Figure 2.3: Gill morphology. (L)Gill length. (R)Root. (V)Ventral. (D)Dorsal.
(M)Midchord. (T1)Gill Tip. (T2)Intersection of gill edge with
−−→
RM . (a) Computa-
tional model of planform A (older nymphs). (b,c) Computational models of planforms
B and C (younger nymphs). (d,e) Actual planforms A and B [58]. (f) Computational
model of all planforms at mean position for rowing (Gill 4). Planforms A, B and C
have an area of 0.638, 0.449, and 0.356 (length squared) and an aspect ratio of 1.57,
2.23, and 2.81 respectively.
flow rate was strongly dependent on the planform. As a result, the gill planform
was treated as an independent variable and three different planforms were selected to
perform the simulations. The first planform corresponds to adult nymphs, while the
other two relate to younger nymphs.
Figure 2.3d shows Planform A, which is the actual gill planform of adult nymphs,
which utilize flapping kinematics. The branching trachea structure emerging at the
25
gill root is clear and interestingly presents a good analogy with vascular systems and
branched water networks. The curved flexural hinge splits the gill into proximal and
distal plates where the proximal plate is the plate in proximity of the gill root. The
vectorized or computational version of Planform A is shown in Figure 2.3a. Point ’R’
is the root of the gill, where it attaches to the mayfly abdomen. Point ’D’ is termed
the dorsal point, i.e. the point of intersection of the dorsal edge with the flexural hinge
V D. Likewise, point ’V’ is termed the ventral point, i.e. the point of intersection of
the ventral edge with the flexural hinge V D. Point ’T1’ is the tip of the gill, defined
as the furthest point from the root ’R’. Point ’M’ is the “mid-chord”, while point ‘T2’
is the extension of segment RM . The importance of all these points arises from the
fact that four of them (’R’, ‘D’, ’V’ and ’T1’) have been tracked experimentally to
extract the gill kinematics, which will be discussed in more detail in the next section.
These points are also shown in Figure 2.2b, in order to give a clear perspective of the
mean orientation of the flapping gills with the respect to the abdomen. Segment V D
is the approximation of the flexural hinge that develops in the gills. The flexural line
again splits the gill into its proximal (defined by points ’R’, ’V’, and ’D’) and distal
(defined by ’V’, ’T1’, and ’D’) plates. The in-plane curvature of the flexural hinge
has been idealized into a straight line owing primarily to the rigidity assumption of
the gill plates.
The actual and vectorized versions of Planform B are shown in Figure 2.3e
and 2.3b respectively. Planform B is representative of younger naiads, which use a
rowing mechanism. Planform B is narrower compared to planform A and does not
possess the bulge seen in Planform A near the root. In other words, the dorsal edge
26
(specifically between point ’R’ and point ’D’) is the major contributing factor in the
planform change. Again, in extracting the gill kinematics, the four points ’R’, ‘D’, ’V’
and ’T1’ were use to track the gill. For rowing kinematics, the flexural hinge angle is
negligibly small and therefore was neglected in rowing simulations.
During the course of study, and after implementing planforms A and B, the
effect of the gill shape on the proposed performance parameters became increasingly
pronounced. A mayfly naiad molts more than a few times, and as it switches from
rowing kinematics to beating kinematics, the gill planform also changes from planform
A to B. While only two gill planforms were supplied by my collaborators in the
Entomology Department, it is not valid to assume that the mayfly abdomen and gill
shapes are ”statistically convergent”. After all, only one mayfly has been dissected
and only one gill per planform has been drawn. When extracting the gill plate
kinematics (Eulerian angles), I had to assume that the gill shape (but not the length)
in Figure 2.3a-b remains similar for all six pairs of moving gills and for any mayfly.
By looking at the average lengths of segments RV , RD and V D in the experiments,
it was clear that this is not true. In other words, there is an amount of uncertainty
in the shape of the gill. In fact, a gill with a higher aspect ratio would yield a better
representation of the rowing kinematics. In light of that, planform C was devised.
Planform C was obtained by scaling down planform B by 25% but only in a direction
perpendicular to RT1. This effectively increased the gill aspect ratio further, which
is illustrated in Figure 2.3c.
The three planforms are superimposed and shown in Figure 2.3f for gill 4 at its
mean location for the rowing kinematics. Note that the local gill axis defined by RT2
27
coincides for all 3 planforms, i.e midpoints of segment DV are collinear. All the gills
shown in Figure 2.3 have been scaled to have the same length in order to highlight
their morphological or numerical differences. However, a typical gill length of older
nymphs (planform A) is of order 0.7mm, while that of younger nymphs (planform B
or C) is of order 0.3mm. Planforms A, B and C have areas of 0.638, 0.449, and 0.356
(length squared) respectively. The aspect ratio of planforms A,B and C, defined as
the gill surface area divided by the square of the gill length, is 1.57, 2.23, and 2.81
respectively.
The following step was to create a surface mesh (Lagrangian mesh) for the
mayfly nymph and the different gills. In other words, the vectorized version of the
mayfly and the gills need to be transformed into the discrete space. In order to ac-
complish this, the AutoCAD drawings have been exported as Standard ACIS Text
and imported into a meshing software (Gambit 2.4.6) where the geometry has been
checked, cleaned then meshed. Details on the numerical methodologies will be pre-
sented in the §3.
2.2 The Prescribed Kinematics
In this section, the kinematic chain used to prescribe the different kinematics is
presented. A fluid-structure interaction model would have required parameters that
were not readily available at the beginning of the study as for example the mass
distribution of the gill plates or elasticity parameters that describe the distal plate
flexing. The interaction between the different gills becomes increasingly important at
28
lower Re and the structural properties will have more pronounced effects. Also, it is
not known that these properties will be the same for all 7 pairs of gills. More interest-
ingly, it is known if these structural properties remain constant throughout ontogeny.
The gills of older nymphs (planform A) are known to be more flexible than those of
younger nymph (planform B). The unavailability and the introduced uncertainty in
the structural properties, if it were available, coupled with the uncertainty in the gill
shape, kinematics and the variability of mayflies would further complicate the analy-
sis of the results. In addition to its computational cost, introducing a fluid-structure
interaction model was neither desirable nor feasible. Therefore, prescribed kinematics
present a much simpler and more direct approach especially with the availability of
the experimental data (kinematics and flow fields).
2.2.1 The Kinematic Chain
In this section, the kinematic chain is discussed in four steps. First a description
of the kinematic chain itself, the different rigid bodies involved and the transformation
matrices resulting from the multi-body system. Second, a short derivation of the
absolute angular velocities of the different bodies. This is followed by a derivation of
the absolute angular accelerations. The last step introduces the equations required to
compute the absolute linear velocity and acceleration of any point in the kinematic
chain.
The kinematic chain is illustrated in Figure 2.4 and consists of 3 levels, where
the 1st level describes the orientation of the mayfly body ’M’ (abdomen, thorax and
29
head) with respect to the inertial Frame of Reference (FOR) ’N’ (0th level), the 2nd
level describes the orientation of the gill proximal plates ’G’ with respect to the mayfly
body (1st level) and the 3rd level describes the orientation of the distal plates ’D’ with
respect to the proximal plate (2nd level). For the flapping kinematics, both the 2nd
and 3rd levels introduce time-dependent angles and are required to fully describe the
gills motion. In the rowing case however, the 3rd level is not required since the entire
gill is treated as one single rigid plate. This was easily executed by setting the hinge
angle γ(t) = 0 when simulating the rowing kinematics.
In the kinematic chain, three body-attached coordinate systems or FORs are
used. These are ’M’, ’G’, and ’D’ for levels 1, 2 and 3 respectively. A series or a
sequence of Eulerian angles is performed in order to compute the location of a child
body in level (l) with respect to its parent body in level (l-1) along with its angular
velocities and accelerations. The first coordinate system M{xm, ym, zm} is attached
to the mayfly body and its origin lies at the symmetry (medial) plane on the midpoint
of the segment connecting the attachment points of gill pair 4. The second coordinate
system is attached to the root any of the proximal plates 1 through 7 and is denoted
by G{xg, yg, zg}. The third coordinate system is attached to any of the distal plates 1
through 7 with its origin at the midpoint of DV and is denoted by D{xd, yd, zd} (see
Figure 2.4). Although a different amount of translation is required to reach each gill
root from 1 to 7 from level 1, the transformation matrix for the rotation sequence is
exactly the same. The same applies for the distal plates. In other words, there is only
one transformation matrix in going from level 1 to level 2 and only one transformation
matrix in going form level 2 to level 3. Only the numerical values of these matrices
30
Figure 2.4: Illustration of kinematic chain. Level 0 {xn, yn, zn}: Inertial FOR. Level 1
{xm, ym, zm}: Mayfly abdomen, thorax and head. Level 2 {xg, yg, zg}: Proximal plate
of gill. Level 3 {xd, yd, zd}: Distal plate of gill. Rotation sequence shown is for the
forward transformation from the inertial FOR down to the distal plates. Level 1 is
useful for orienting the mayfly inside the Eulerian domain and is not time-dependent.
Level 2 has 3 time-dependent angles ψ(t), θ(t) and φ(t) signifying the beating angle,
the stroke plane deviation angle and the angle that modifies the pitch. Level 3 has
only one time-dependent angle γ(t), which signifies the hinge flexion angle. Right
hand rule is strictly followed to assign signs to different angles.
31
will differ depending on the type of kinematics used and the gill number. Worth
mentioning is that gill pair 7 is not moving and that gills 1 through 6 have similar
kinematics but different values for , β, ψ(t), θ(t), φ(t) and γ(t). FOR ’M’ is used
to describe the location of the mayfly body with respect to the inertial FOR ’N’.
The transformation to level ’M’ adds more flexibility because it allows for different
orientations of the mayfly body within the Eulerian domain. The transformation
was implemented to also allow future extensions of the code (insect flight). The
transformation from the inertial FOR ’N’ to the mayfly FOR ’M’ is achieved by a 180◦
rotation with respect to the 1 axis, followed by a 3-2-1 Eulerian angle sequence with
ψm(t)-θm(t)-φm(t) . This is the typical flight mechanics coordinate transformation
(yaw, pitch and roll respectively) and yields a singularity for pitch angles multiple of
pi/2. In the current work, the values of ψm(t), θm(t), and φm(t) were kept constant
at −pi/2, 0, and pi/2. The transformation matrix from ’N’ to ’M’ is given by:
[TMN ] =








cθcψ −cθsψ sθ
sφsθcψ − cφsψ −sφsθsψ − cφcψ −sφcθ
cφsθcψ + sφsψ −cφsθsψ + sφcψ −cφcθ








=








0 1 0
0 0 −1
−1 0 0








(2.1)
where s() and c() are the sine() and cosine() functions respectively. Level 1 may be
interpreted as follows: 1) a positive rotation about the 1 axis in level 1 translates to
a positive rotation about the inertial 2 axis. 2) a positive rotation about the 2 axis
in level 1 translates to a negative rotation about the inertial 3 axis. 3) a positive
rotation about the 3 axis in level 1 translates to a negative rotation about the inertial
1 axis.
FOR ’G’ is used to describe the location of any of the proximal plates 1 through
32
7 G{xg, yg, zg} with respect to the mayfly body M{xm, ym, zm}. The transformation
from the mayfly system ’M’ to the proximal plates system ’G’ is achieved in 3 steps
with a total of 6 rotations as illustrated in Figure 2.5. First, a rotation of pi about 1.
This was primarily done to restore the axis orientation used in previous experimental
work. Second, a 3-2 rotation of  (stroke plane lateral offset angle) and β (stroke
plane inclination angle), which describe the location of the stroke plane (SP). Finally
a 3-2-1 Eulerian rotation of ψ(t), θ(t), and φ(t), which signifies the in-plane beating
angle, the stroke plane deviation angle and the pitch respectively. The transformation
matrix of this 1-3-2-3-2-1 sequence is:
[TGM ] =








T11 T12 T13
T21 T22 T23
T31 T32 T33








(2.2)
where:
T11 = −cθ(ssψ − cβccψ)− csβsθ
T12 = sβssθ − cθ(csψ + cβcψs)
T13 = cβsθ + cψsβcθ
T21 = −cφ(cψs+ cβcsψ)− sφ(sθ(ssψ − cβccψ)− csβcθ)
T22 = −sφ(sθ(csψ + cβcψs) + sβcθs)− cφ(ccψ − cβssψ)
T23 = −sφ(cβcθ − cψsβsθ)− cφsβsψ
T31 = sφ(cψs+ cβcsψ)− cφ(sθ(ssψ − cβccψ)− csβcθ)
T32 = sφ(ccψ − cβssψ)− cφ(sθ(csψ + cβcψs) + sβcθs)
T33 = sβsφsψ − cφ(cβcθ − cψsβsθ)
33
FOR ’D’ is used to describe the location of any of the distal plates 1 through
7 D{xd, yd, zd} with respect to its parent proximal plate G{xg, yg, zg}. Level ’D’ is
simply a 3-2 transformation with λ and γ(t), where λ is the angle which aligns the 2
axis with the flexural hinge line, DV , and γ(t) is the time-dependent hinge flexural
angle, which defines the deflection or the flexion of the distal plate with respect to
its parent proximal plate. The angles are illustrated in Figure 2.4. The resulting
transformation from the ’G’ FOR to the ’D’ FOR is:
[TDG] =








cγcλ cγsλ −sγ
−sλ cλ 0
cλsγ sγsλ cγ








(2.3)
Equations (2.1), (2.2), and (2.3) define the transformation matrices between any two
consecutive levels and therefore enable the transformation of any vector quantity
between the inertial FOR ’N’ and any other FOR (’M’, ’G’, or ’D’) using equations
similar to the following:
[TNM ] = [TMN ]T (2.4)
[TNG] = [TNM ][TMG] = [TGN ]T = [TMN ]T [TGM ]T (2.5)
[TND] = [TNM ][TMG][TGD] = [TDN ]T = [TMN ]T [TGM ]T [TDG]T (2.6)
The next step in the kinematic analysis of this mechanical system is to compute
the angular velocities of the different bodies in terms of the time derivatives of the
Eulerian angles. In general, the angular velocity of some frame ’Y’ with respect to
another frame ’X’ may be expressed in terms of the time derivatives of those Eulerian
34
Figure 2.5: Illustration of the forward transformation from the mayfly local FOR
M{xm, ym, zm} to any local left gill FOR G{xg, yg, zg}. The sequence is a 1 − 3 −
2 − 3 − 2 − 1 Eulerian sequence with pi −  − β − ψ(t) − θ(t) − φ(t), where  is the
lateral offset angle, β is the stroke plane inclination angle, ψ(t) is the beating angle,
θ(t) is the stroke plane deviation angle and φ(t) is then angle that modifies the pitch.
The illustration shown is for  = 100◦, β = −60◦, ψ = 30◦, θ = −15◦, φ = −45◦.
(a) Isometric view of mayfly naiad with its local FOR, left row of gills omitted.
(b) Isometric view (matching view in (a)) of the stroke plane {xsp, ysp, zsp} after the
1−3−2 rotation of pi,  = 100◦, β = −60◦. (c) Perpendicular view of the stroke plane
showing the in-plane gill (ψ = 0◦, θ = 0◦, φ = 0◦). (d) Gill location {x2, y2, z2} after
beating of ψ = 30◦. (e) Gill location before {x2, y2, z2} and after {x3, y3, z3} a stroke
plane deviation of θ = −15◦. (f) Gill location before {x3, y3, z3} and after {xg, yg, zg}
a pitch of φ = −45◦. Gill at {x2, y2, z2} is also shown (axis omitted for clarity).
35
angles used to rotate between ’Y’ and ’X’ in some frame ’B’ in the form XωYB = MΦ˙
where the definition frame ’B’ is usually one of the frames ’X’ or ’Y’. This is obtained
by expressing the components of the angular velocity vector ω in terms of only one
FOR rather than in terms of all the FORs involved in the rotation sequence, which
normally results from an Eulerian sequence. The angular velocity of the proximal gill
plate ’G’ with respect to the mayfly ’M’ expressed in the mayfly frame ’M’ is denoted
MωGM = M(θ,ψ)Φ˙(φ˙,θ˙,ψ˙) and reads:
MωGM =








cβccθcψ − scθsψ − sβcsθ −cβcsψ − scψ sβc
−cβscθcψ − ccθsψ + sβssθ cβssψ − ccψ −sβs
sβcθcψ + cβsθ −sβsψ −cβ











φ˙
θ˙
ψ˙



(2.7)
Similarly, the angular velocity of the distal gill plate ’D’ with respect to the proximal
gill plate ’G’ expressed in the proximal gill frame ’G’ is denoted GωDG = M(λ)Φ˙(γ˙) and
reads:
GωDG =








0 −sλ 0
0 cλ 0
0 0 0











0
γ˙
0



=



−sλγ˙
cλγ˙
0



(2.8)
Equations (2.7) and (2.8) fully describe the angular velocities of the proximal and
distal plates of the gill respectively. Specifically it describes the relative angular
velocity between a child (level l) and its parent (level l − 1). Although the angular
velocities are expressed in different FORs, the transformation matrices introduced in
(2.1), (2.2) and (2.3) along with the relations in (2.4) allows the expression of the
angular velocities in any FOR and in the inertial FOR ’N’ in particular. The absolute
36
angular velocities are calculated by use of the addition theorem [28,47]:
NωlN =
Nωl−1N +
l−1ωlN 1 ≤ l ≤ 3 (2.9)
where l is the level number.
The next step in the kinematic analysis of this mechanical system is to compute
the absolute angular accelerations of the different bodies in terms. With the angular
velocity of the different frames known, it is possible to calculate the angular accel-
eration with the respect to the inertial FOR ’N’ using the fundamental composition
rule for the absolute angular acceleration [7]:
NαlN =
Nαl−1N +
l−1αlN +
Nωl−1N ×
l−1ωlN 1 ≤ l ≤ 3 (2.10)
where l is the level. The angular acceleration of a child with respect to its parent is
computed using Coriolis theorem as follows:
NαMN =
N d
dt
NωMN = 0 (2.11)
MαGM =
N d
dt
MωGM
=
M d
dt
MωGM +
NωMM ×
MωGM
= M(θ,ψ)Φ¨(φ˙,θ˙,ψ˙) + M˙(θ,ψ)Φ˙(φ˙,θ˙,ψ˙) (2.12)
GαDG =
N d
dt
GωDG
=
G d
dt
GωDG +
NωGG ×
GωDG
= M(λ)Φ¨(γ˙) + M˙(λ)Φ¨(γ˙) +
NωGG ×
GωDG (2.13)
Equations (2.10) through (2.13) completely define the absolute angular accelerations
37
Figure 2.6: Schematic of position vectors in different levels. Pl: Point in level l. N,
M, G and D: Inertial, Mayfly, Gill and Distal plate frames of reference. irPl: Position
vector from FOR i to point Pl. ivPl, iaPl: Velocity and acceleration of point Pl with
respect to frame i respectively.
of the entire kinematic chain with the help of the transformation relations in (2.4).
With the absolute angular velocity and acceleration of all bodies (proximal ’G’
and distal ’D’ plates) determined, the last step is to compute the absolute position,
linear velocity and linear acceleration of any point in the chain. With reference to
Figure 2.6, the position vector for points ’P1’, ’P2’ and ’P3’ lying on the mayfly ’M’,
38
proximal plate ’G’ and the distal plate ’D’ respectively can be expressed as:
NrP1 = NrM + MrP1 (2.14)
NrP2 = NrM + MrG + GrP2
= NrG + GrP2 (2.15)
NrP3 = NrM + MrG + GrD + DrP3
= NrD + DrP3 (2.16)
The corresponding linear velocity of any point lying on these rigid bodies may be
respectively expressed as:
NvP1 = NvM + NωM × MrP1 (2.17)
NvP2 = NvG + NωG × GrP2 (2.18)
NvP3 = NvD + NωD × DrP3 (2.19)
Similarly, the corresponding linear acceleration of any point lying on these rigid bodies
may be respectively expressed as:
NaP1 = NaM + NαM × MrP1 + NωM × (NωM × MrP1) (2.20)
NaP2 = NaG + NαG × GrP2 + NωG × (NωG × GrP2) (2.21)
NaP3 = NaD + NαD × DrP3 + NωD × (NωD × DrP3) (2.22)
In investigating the fluid dynamics of mayfly naiads, two sets of mayfly kine-
matics representative of flapping and rowing were under consideration, both extracted
from the experimental measurements [57]. A description of both sets follows. In both
sets, the beating cycle is split into two parts: 1) retraction, where the proximal plate
39
is moving towards the rear of the insect, and 2) protraction, where the proximal plate
is moving towards the head. This notion is most comparable to the down-stroke and
the up-stroke half-strokes in flapping flight literature because it implies a direction of
motion rather than a function. The notion of a power stroke and a recovery stroke is
not useful in flapping kinematics because both half-strokes are equally important and
both half-strokes contribute to the induced micro-pumping. The power and recovery
half-strokes notion, however, can be applied to the rowing kinematics since each half
stroke has a distinguishable function.
In both sets, the kinematics of gill pairs 1 to 6 in the array are similar in terms of
magnitude and frequency, and only differ in terms of their mean positions with respect
to the body; gills towards the head of the insect are positioned closer to the abdomen
(higher lateral offset angle, , for SP definition). In both sets, gill pair number 7 is
not moving and is, therefore, believed to act as a buffer to enhance performance. The
phase lag between any two consecutive gills is approximately 90◦, with the gills at the
back having the lead. Whenever the mayfly naiad starts to beat its gills, it starts by
beating gill pair 6. The metachronal wave starts posteriorly at gill pair number 6 and
progresses forward towards the head of the animal. The actual frequency of the gills
varies from 20Hz to 37Hz and no correlation could be made with the mayfly size or
age [58]. Noteworthy is that experimental results from an idealized rowing robotic
model of the mayfly naiad at Reω = 17.4 showed that the 90◦ phase lag produces the
highest efficiency for the rowing case when compared to phase lags of 0◦, 180◦ and
270◦ [44]. In other words, it is possible that the 90◦ phase lag is nature’s choice to
optimize the induced micro-pumping. Both sets of kinematics are presented next.
40
2.2.2 The Flapping Kinematics
The flapping kinematics of gill 4 is illustrated in Figure 2.7. Part (a) of the figure
shows an entire period of the gill as projected onto an xy-plane (the median plane
or the symmetry plane), while parts (b) and (c) show gill 4 during retraction and
protraction respectively. The projection onto the xy-plane is a good representation
of the actual kinematics, i.e. antero-posterior (along y) and dorso-ventral (along x)
motion represents the dominant motion. Frames 1 and 3 depict the extreme positions
of the proximal plate, while frames 4 and 2 depict mid-retraction and mid-protraction
respectively. The first thing to note about the flapping kinematics is that the motion
of the proximal plate is almost symmetric. It takes 48% of a cycle to protract and 52%
of to retract. In addition, looking at any two time frames in Figure 2.7a that lie in the
same vertical column (e.g. 0.55 and 0.65) the segments DV (red) and RM (black) are
always parallel, which again highlights the symmetry of motion of the proximal plate
around its mean position. The segments in lighter dashed shades (gray and light red),
which represent the position of the gill 0.025 earlier in the period, demonstrate that
the acceleration of the proximal plate is the same during retraction and protraction.
For this particular gill, the mean pitch angle is −92 ◦ with a range of only 6 ◦. The
beating angle has a mean of −22 ◦ and an amplitude of 22 ◦. The range of the stroke
plane deviation is 3 ◦. It is evident from these numbers that the motion is dominated
by the beating angle with relatively little pitch variation or stroke plane deviation.
Finally, the motion of the distal plate relative to the hinge can also be seen in the
figure. The hinge flexion angle varies from 10 ◦ at t/τ = 0.45 to 65 ◦ at t/τ = 0.8. As
41
a result, the mean position of the distal plate in not parallel to the mean position of
the proximal plate. A plot of the different angles describing the motion of all the 6
gills is shown Figure 2.8.
2.2.3 The Rowing Kinematics
The rowing kinematics for the same gill (number 4) are illustrated in Figure 2.9.
The stroke plane in this case is inclined −53 ◦ from the yz plane. Compared to
flapping, there is no symmetry between the protraction and retraction part of the
cycle. It takes 38% of the cycle for the gill to retract and 62% of the cycle to protract.
The beating angle has a range of 47 ◦ over the entire period, which is more than twice
the flapping beating range (22 ◦). The presence of much higher accelerations during
retraction is evident, as manifested by the larger lag between the black and gray
lines (i.e. t = 0.3 − 0.4). In addition, the pitch angle has a range of 37 ◦, which is
more than six times that of the flapping. The mean pitch is −90 ◦ during retraction,
−70 ◦ during protraction, and −78 ◦ for the entire period. In summary, the rowing
motion is characterized by higher beating, accelerations, pitching and stroke plane
deviation. The experiments have shown that the mean flow field is directed parallel
to the stroke plane. In fact, the retraction is the power stroke where the fluid is being
pushed towards the back and the bottom of the insect with an average pitch of −90 ◦,
while the protraction is obviously a recovery stroke where the gill plate is feathering
into the mean flow with a mean pitch of −70 ◦ to -possibly- attempt to minimize the
work done by the mayfly and enhance the efficiency. A plot of the different angles
42
Figure 2.7: Prescribed mayfly flapping kinematics for gill 4. (a) Kinematics projected
onto xy plane. Dashed trails in lighter shades are 0.025 earlier in the beating period
(τ = 1). (b) Retraction (3D perspective). (c) Protraction (3D perspective). Inset
illustrates the 3D perspective for parts (b) and (c). 1©, 3©:Extreme positions of
proximal plate. 2©:Mid-protraction. 4©:Mid-retraction.
43
Figure 2.8: Eulerian angles for mayfly flapping kinematics (◦). Top to bottom: Stroke
or beating angle ψ1−6(t). Stroke plane deviation angle θ1−6(t). Pitch angle φ1−6(t).
Hinge flexion angle γ1−6(t).
44
Figure 2.9: Prescribed mayfly rowing kinematics for gill 4. (a) Kinematics projected
onto xy plane. Dashed trails in lighter shades are 0.025 earlier in the beating period
(τ = 1). (b) Retraction (3D perspective). (c) Protraction (3D perspective). Inset
illustrates the 3D perspective for parts (b) and (c). 1©, 3©:Extreme positions of
proximal plate. 2©:Mid-protraction. 4©:Mid-retraction.
45
Figure 2.10: Eulerian angles for mayfly rowing kinematics (◦). Top to bottom: Stroke
or beating angle ψ1−6(t). Stroke plane deviation angle θ1−6(t). Pitch angle φ1−6(t).
Hinge flexion angle γ1−6(t).
46
describing the motion of all the 6 gills is shown Figure 2.10. Noteworthy is that gill
1 is almost not moving and that the kinematics of gill 2 are not that similar to the
remaining gills 3 to 6.
2.3 Governing Equations
Ideally, the physics of the problem would be investigated by solving for mass
conservation, momentum conservation, and dissolved gases (oxygen) transport. How-
ever, the diffusion properties through the gills surface as well as the exact physiology
inside the tracheal tubes is not well studied. In other words, one would have to solve
for oxygen consumption (by simulating metabolism chemically for example) inside
the mayfly in order to force the correct transport boundary conditions on the sur-
rounding water. In addition, the exact solution of oxygen transport at such low Re
and high P e´ is cost-prohibitive. The estimated computational time is about 3 years if
the same parametric spaces were to be carried out, not including the time needed to
develop and validate the numerical codes. Therefore, including transport was beyond
the scope of the current work, and, as a first approach, only mass and momentum
conservation around the mayfly is considered. In light of that, the PDEs governing
the problem are simply the continuity equation and the Navier-Stokes equations for
incompressible Newtonian fluids, which in Cartesian coordinates read:
∂ui
∂xi
= 0 (2.23)
∂ui
∂t
+
∂(uiuj)
∂xj
= −
∂p
∂xi
+
1
Reω
∂2ui
∂xj∂xj
+ fi i = (1, 2, 3) (2.24)
47
where all the independent (x, t) and dependent (u, p) variables are dimensionless,
xi is the Cartesian coordinate in the ith direction, ui is the velocity component in
the ith direction, p is the pressure, fi represents an external body force field used
to implement the correct boundary conditions (will be introduced in §3). Reω is
the oscillating Reynolds number defined as Reω = fL2/ν where f is the gill beating
frequency, L is the gill length, and ν is the kinematic viscosity. Both equations 2.23
and 2.24 were non-dimensionalized using the gill length, L, as a characteristic length
scale, the beating period, τ (= 1/f), as a characteristic time scale, and U = Lf
as a characteristic velocity scale with inertial scaling for the pressure gradient term
(pref = ρL2f 2).
In the process of non-dimensionalizing the momentum equation, it is customary
to relate the reference time (τ or 1/f), the reference length (L) and the reference
velocity (U) (i.e. U = fL). This results in only one dimensionless group representing
the ratio between the inertial forces and the viscous forces, which is the famously
known Reynolds number. While this is may be a mathematical convenience, it may
truncate some of the physics of the problem. Keeping the reference time, length
and velocity unrelated will result, for example, in St = fL/U sitting next to the
non-dimensionalized local acceleration. Now, it is common to associate St with the
vortex shedding frequency (f) as in the von Ka´rma´n street vortex behind bluff bodies
or the inverted von Ka´rma´n street vortex associated with flapping bodies or even in
flow meters. The question that arises now is what does St mean if there is no vortex
shedding? i.e. the vortices are always attached to the body (which happens at
the lower end of Reω in my parametric space) or there are no vortices at all. In
48
light of that, I found it more insightful to distinguish between Eulerian inertia and
convective inertia rather than just inertia. This may be thought of as an extension of
the well-understood local and convective accelerations by multiplying by the ”mass”
of the fluid particle. Therefore, with the correct normalization and for small St, the
fluid particle changes its inertia as it flows (convects) along its pathline (convective
inertia > local inertia). On the contrary, for large St, the flow moves like a plug
and there is a collective movement of the flow (local inertia > convective inertia).
An example of large St flows is the unsteady Stokes flow. In unsteady flows or
steady oscillatory flows, when St is close to unity (low end of Reω in experimental
studies [58]), there is a competition or balance between the Eulerian (also local) inertia
and the convective inertia. For more insight about dimensional analysis, scaling, and
orders of magnitude, the reader may consult [2]. Noteworthy is that St may also
be regarded as the ratio between the frequency parameter (L2/ντ) and Reynolds
number. In light of that, Reω introduced in (2.24) may be interpreted in many ways:
• ratio of the Eulerian inertia (ρL3U/τ) to the viscous forces (µ(U/L)L2)
• ratio of the convective inertia (ρL3U2/L) to the viscous forces (µ(U/L)L2)
• ratio of a characteristic length scale (L) to a diffusive length scale (Oseen length
= ν/Lf)
All the mathematical tools are now ready for implementation in a discrete envi-
ronment. The details of the numerical scheme along with the immersed boundary
formulation will be introduced in §3. In the next section, two different metrics used
to assess the micro-pumping of mayfly naiads are introduced.
49
2.4 Performance Metrics
In insect locomotion (i.e. flapping flight), it is usually straightforward to define
performance metrics based on the generation of aerodynamic forces such as drag and
lift. Simple models examining the mechanical efficiency for propulsion (work done to
move the animal divided by the work done on the fluid) have shown that a switch
from flapping to rowing kinematics is favored as one decreases the Re below 20 [12].
Mayfly nymphs, however, are not generating forces for the purpose of propulsion, but
instead are flapping their gills to circulate water around them and enhance oxygen
absorption. Therefore, it was necessary to create or look into novel performance
metrics that could be computed based on the current model (mass and momentum
conservation). With the nature of the flow field, a control volume approach such as
the one shown in Figure 2.11 was devised and two performance metrics have been
introduced.
Figure 2.11: Schematic of control volume. Control volume encompasses the entire
array of moving gills.
50
2.4.1 Mechanical Efficiency
One can define an efficiency similar to that of a hydraulic pump, which can be
derived from the mechanical energy conservation equation:
d
dt∗
∫
E∗dV ∗ +
∫
E∗u∗ · dA∗ =
∫
u∗i τ
∗
ijdA
∗
j −
∫
φ∗dV ∗ (2.25)
where E∗ = 0.5ρu∗iu
∗
i is the kinetic energy per unit volume, τ
∗
ij = −p
∗δij + 2µe∗ij is
the stress tensor, φ∗ = 2µe∗ije
∗
ij is the rate of viscous dissipation, e
∗
ij = 0.5(∂u
∗
i /∂x
∗
j +
∂u∗j/∂x
∗
i ) is the strain rate tensor and the star superscript denotes dimensional vari-
ables [41]. Equation (2.25) is non-dimensionalized using the same reference parame-
ters introduced in (2.23) and (2.24) (ρL5f 3) which yields:
d
dt
∫
EdV +
∫
Eu · dA =
∫
uiτijdAj −
∫
φdV , (2.26)
where E = 0.5uiui, τij = −pδij + 2eij/Reω, φ = 2eijeij/Reω, eij = 0.5(∂ui/∂xj +
∂uj/∂xi), where u, p, x and t are dimensionless as previously introduced in (2.23) and
(2.24). Equation (2.26) is then time averaged over an entire beating period. The term
corresponding to adding power to the control volume is identified, which is the time-
averaged rate of work done by the moving gills. Also, the two terms corresponding
to the rate of work done by the control volume are identified, which are the time-
averaged kinetic energy flux through the control volume and the time-averaged rate
of work done by the control volume on the remaining outer domain. The mechanical
efficiency can now be defined as the ratio between the time-averaged output power
and the time-averaged input power, where the difference between these is simply the
time-averaged rate of viscous dissipation or the rate of energy loss throughout the
51
control volume averaged over an entire beating cycle. The mechanical efficiency can
be expressed as follows:
η(%) =
∫
outer Eu · dA−
∫
outer uiτijdAj
∫
gills uiτijdAj
× 100 (2.27)
where ’outer’ refers to the surface defined by the six planes corresponding to the rect-
angular control volume in Figure 2.11, and ’gills’ refer to the surface of the gills. The
overbar denotes time-averaged quantities. Note that the kinetic energy flux through
the control surfaces around the gills is zero owing to the assumed zero thickness of
the gill plates. The mechanical efficiency, η, is analogous to that of a hydraulic pump
where the denominator is the shaft input power while the numerator is the energy
flux in the fluid and stress work done on the outer control volume boundary.
2.4.2 Specific Mass Flow Rate
With reference to the same control volume shown in Figure 2.11, a second
performance parameter is introduced; the specific mass flow rate, ξ. ξ is defined as
follows:
ξ =
m˙
W˙ gills
, (2.28)
where
m˙ =



∫
outer ρu · dA if u · dA > 0
0 otherwise
, (2.29)
and
W˙gills =
∫
gills
uiτijdAj (2.30)
52
where m˙ is the mass flow rate leaving the control volume, and W˙gills is the rate of
work done by the six pairs of gills. ξ is interpreted as the dimensionless time-averaged
mass flow rate towards the mayfly divided by the dimensionless time-averaged rate of
work done by the gills. It can be argued that the mass flow rate towards the mayfly is
associated with the oxygen absorption by the gills. Therefore, the higher ξ, the better
the performance, since in an average sense, higher oxygen absorption is achieved per
unit power. The results of both metrics will be presented in §4 and §5.
53
Chapter 3: Numerical Methodologies and Validation
In this chapter, the numerical methodologies and the computational setup used
to solve for the laminar flow around mayfly nymphs are presented. This includes
the domain spatial discretization, the time advancement scheme and the immersed
boundary formulation used to implement the boundary conditions on the fluid. This
will be followed by comparisons with in-vivo and robotic PIV experiments.
3.1 Spatial and Temporal Discretization
The equations governing the laminar flow of water around the mayfly were
previously introduced in §2.3 and read:
∂ui
∂xi
= 0 (3.1)
∂ui
∂t
+
∂(uiuj)
∂xj
= −
∂p
∂xi
+
1
Reω
∂2ui
∂xj∂xj
+ fi i = (1, 2, 3) (3.2)
where fi is a the forcing function in the ith direction and is non-zero near the solid
boundaries.
The above conservation laws are solved using a finite differencing approach with
a standard second-order central-difference scheme on a Cartesian staggered grid. In
a staggered arrangement, the pressure and other scalar variables are located at the
54
center of the grid cell, while velocity components are located at the cell face centers.
The most important advantage of staggered grids over collocated grids is the stronger
coupling between the pressure and the velocities, which avoids the so-called odd-even
decoupling between the pressure and the velocity. The numerical approximation on
a staggered grid is also conservative of kinetic energy [27]. Staggered grids, however,
can be trickier to implement, since 4 grids are involved in computing the different
gradients.
As for time advancement, a fractional time step method is used to integrate
the unsteady momentum equation. In particular, the low-storage third-order Runge-
Kutta (RK3) scheme is used:
uˆi − u
n−1
i
∆t
= γkA(u
n−1
i ) + ρkA(u
n−2
i )− αk
∂pn−1
∂xi
+ fn−1i (3.3)
∂2φn
∂xi∂xi
=
1
αk∆t
∂uˆi
∂xi
(3.4)
uni = uˆi − αk∆t
φn
xi
(3.5)
pn = pn−1 + φn (3.6)
where n is the time subscript and ranges from 1 to 3 for the RK3 scheme, ∆t is the
time step and fn−1i is the forcing function. The velocity uˆi is the intermediate velocity,
which does not satisfy the continuity equation. Therefore, the Poisson equation (3.4)
is solved and the scalar φ is used to project uˆi into a solenoidal (divergence-free)
field. The spatial operator A contains the convective and viscous fluxes which are
evaluated at the cell faces. The discretization of the above equations on a staggered
grid in the computational space is described in detail in [8, 80]. The coefficients for
55
the RK3 scheme are:
α1 = 8/15, γ1 = 8/15, ρ1 = 0
α2 = 2/15, γ2 = 5/15, ρ2 = −17/60
α3 = 1/15, γ3 = 3/4, ρ3 = −5/12
with:
3∑
k=1
αk =
3∑
k=1
(γk + ρk) (3.7)
In addition, the following stability criterion (or the generalized CFL number
including the time step constraint from the viscous terms) is adopted [8, 80]:
CFL = ∆t
[
|u|
∆x
+
|v|
∆y
+
|w|
∆z
+ 2ν
(
1
∆x2
+
1
∆y2
+
1
∆z2
)]
(3.8)
Theoretically, the stability limit for RK3 is
√
3. However, the actual CFL number
may be lower because the cross-derivatives are not included in the stability analysis.
Consequently, in all simulations, a CFL = 1.0 was found to yield a stable solution.
The Poisson equation shown in vector form in (3.4) can be expressed in discrete
form in Cartesian coordinates as:
(
∂2
∂2x
+
∂2
∂2y
+
∂2
∂2z
)
φi,j,k = gi,j,k =
1
αn∆t
∂uˆi
∂xi
(3.9)
where n is the sub-step index for RK3 and φ = pn−pn−1. With a second order spatial
discretization, the 3-dimensional discretized version of the Poisson equation results in
a hepta-diagonal coefficient matrix, which is extremely expensive to invert directly.
A common approach to alleviate the computational cost is to apply a Discrete Fast
Fourier Transform (DFFT) to the Poisson equation (3.9) in order to reduce the hepta-
diagonal coefficient matrix into a penta-diagonal or even tri-diagonal matrix. For
56
example, applying a DFFT in the y-direction transforms (3.9) into a set of two-
dimensional Helmholtz equations in the uncoupled wave number space:
(
∂2
∂2x
+ k′l +
∂2
∂2z
)
φˆi,j,k = gˆi,j,k (3.10)
where k′l is the modified wave number and is expressed as:
k′l =
2
∆y2
[
1− cos
(
2pil
Ny
)]
(3.11)
where l is the wave number, ∆y is the cell size in the y-direction and Ny is the number
of grid cells in the y-direction not including the ghost cells. Since the above equation
includes derivatives only in the x-direction and z-direction, the resulting coefficient
matrix is penta-diagonal. Hence, a set of Ny uncoupled equations in wave space
may be solved independently. This procedure was implemented using the FFTPACK
library and FISHPACK library, in which a generalized cyclic reduction algorithm
(’BLKTRI’) is used [61, 62]. The solution is transformed back to the physical space
by taking the inverse DFFT of φˆ.
The use of DFFT is some direction comes with a known caveat, in that the flow
must be homogeneous and periodic in that direction. In addition, the numerical grid
must be uniform in that direction. Noteworthy is that DFFT routines are optimized
when the number of points can be factored to the smallest prime numbers.
3.2 Treatment of Immersed Boundaries
Over the past decades, a variety of non-boundary-conforming methods with
various degrees of accuracy and complexity have been proposed, among which are
57
the Cartesian cut-cell formulations and the immersed boundary (IB) methods. The
IB formulation was pioneered by Peskin [51, 52] and has been receiving an increased
attention over the past decade [4, 8, 68, 69,80–82].
IB methods provide a way to impose no-slip boundary conditions on immersed
surfaces that do not conform to the Eulerian grid. The velocity field is reconstructed
around these surfaces by applying a force field (fi) on the discrete momentum equa-
tions in order to satisfy the correct no-slip boundary condition. The most important
feature of the IB methods is that the Eulerian grid (grid where fluid conservation laws
must be satisfied) need not to relate in a special way to the Lagrangian grid (grid
defining the immersed boundary), except that the moving Lagrangian grid must be
sufficiently dense with respect to the fixed Eulerian grid. If the Lagrangian grid
happens to be too coarse with respect to the Eulerian grid, the no-slip boundary
conditions will be incorrectly represented. The simplicity of the IB method facilitates
the numerical integration of the Navier-Stokes equations and is usually applied within
a fractional time step method.
In ’direct forcing’ IB methods, the forcing is done in the discrete equations de-
fined on the Eulerian grid, and it is defined such that the boundary conditions are
satisfied on the boundary itself [4, 26]. These methods are attracting because their
ease of implementation in existing finite differences or finite volumes formulations.
However, moving immersed boundaries introduce additional complications to direct
forcing methods, which leads to hydrodynamic forces that lack smoothness and, there-
fore, stability [64, 81]. A different approach in IB methods is to compute the forcing
function on the Lagrangian (moving) grid and then ’transfer’ this function to the Eu-
58
lerian grid. This approach is introduced next, where the transfer is performed using
the Method of Least Squares (MLS).
3.2.1 MLS Reconstruction
In the formulation proposed by Uhlmann [64], the direct forcing function is com-
puted on each Lagrangian marker, rather than on the Eulerian grid nodes. However,
Uhlmann’s formulation was tailored to a specific problem and provided only integral
quantities of the hydrodynamic forces. Based on the formulation by Uhlmann and
ideas presented in [42,49], Vanella [68,69] proposed a direct forcing scheme that uti-
lizes a versatile MLS approximation in order to build the transfer functions between
the Eulerian and the Lagrangian grids. The MLS approximation can be applied to ar-
bitrary moving or deforming bodies. Vanella also proposed a method to compute the
local traction forces. The overall formulation utilizes very compact stencils and, with-
out compromising accuracy and robustness, gives results that are identical to ‘sharp’
direct forcing methods. A short description of the formulation is given hereunder.
Further details can be found in [68,69] and [49].
First, the velocity component in the ith direction on a Lagrangian marker is
approximated in its support domain by:
U˜i(x) =
m∑
j=1
pj(x)aj(x) = pT (x)a(x) (3.12)
where U˜i(x) is the approximated velocity of the Lagrangian marker in the ith direction,
pT (x) is the basis functions vector of length m, a(x) is the vector of coefficients
and (x) is the position vector of the Lagrangian marker. A linear basis, pT (x) =
59
[1 x y z], was found to be a cost efficient choice that preserves the second order
spatial discretization scheme. Also, a support domain of approximately 2.4 the grid
spacing in each direction was found to be sufficient.
Second, to obtain the coefficient vector a(x), a weighted L2-norm is defined:
J =
ne∑
k=1
W (x− xk)
[
pT (xk)a(x)− uˆki
]2
(3.13)
where xk is the position vector of the Eulerian point k in the Eulerian interpolation
stencil, ne is the total number of Eulerian grid points in the interpolation stencil, uˆki
is the intermediate velocity in the ith direction (calculated from (3.3)) for an Eulerian
grid point k, and W (x− xk) is a known weighting function that uses cubic splines.
The solution of the minimization problem yields the shape function Φ(x) of the
Lagrangian marker l, which is a column vector of length ne. The velocity component
of the Lagrangian marker in the ith can now be calculated from:
U˜i(x) =
ne∑
k=1
φlkuˆ
k
i = Φ
T (x)uˆki (3.14)
where ne is the total number of Eulerian grid points in the interpolation stencil and
φlk is the value of the shape function corresponding to the Lagrangian marker l and
Eulerian grid point k.
Third, the amount of forcing required on the Lagrangian marker in the ith direction
is computed from:
F ni =
Udi − U˜i
∆t
(3.15)
where Udi is the desired velocity of the Lagrangian marker (boundary condition) in
the ith direction.
60
With that, the forcing function is computed on all the Lagrangian markers and the last
step is to transfer the forcing function F ni back to the Eulerian grid. By requiring that
the total force acting on the fluid not to change by the transfer, the forcing function
on the Eulerian grid points can be computed from:
fki =
nl∑
l=1
clφ
l
kF
l
i (3.16)
where cl is a scaling factor related to the average volume associated with the La-
grangian markers and the average volume of the Eulerian grid cells [49, 68,69].
In practice, the MLS formulation works in the following steps: 1) Locate a
Lagrangian marker to enforce its no-slip boundary condition. 2) Locate the nearest
Eulerian neighbor to that Lagrangian marker using bisection routines. 3) Build an
interpolation stencil that consists of a minimum of 5 points for 3D problems (3.12).
4) Find the minimum weighted L2-norm of the difference between the Eulerian inter-
mediate velocity field uˆi and their regressed values U˜i (defined by (3.13)), which yields
the MLS shape functions Φ(x). 5) Use the shape functions to transfer (interpolate)
the Eulerian velocity field uˆi to the Lagrangian marker U˜i (3.14). 6) Compute the
amount of Lagrangian forcing required to satisfy its boundary condition (3.15). 7)
Rescale the shape functions by way of conserving the total force to be transferred.
8) Transfer (interpolate) the forcing function back to the Eulerian stencil (3.16). 9)
Loop on all the Lagrangian markers and sum the forcing field on the neighboring
Eulerian points. 10) Repeat for each direction i.
61
3.2.2 Evaluation of Surface Forces
Non-boundary conforming formulations are faced with complications in com-
puting the hydrodynamic forces generated by the surrounding fluid. This is because
the Eulerian and the Lagrangian grid are almost never aligned. Once the force field is
computed from (3.16), Uhlmann [64] proposed a formulation for computing the hydro-
dynamic forces on rigid bodies provided that all interior points are treated properly.
However, Vanella [68, 69] proposed a different approach for the general case of mov-
ing or deforming bodies, which is based on the MLS formulation discussed previously.
The hydrodynamic force per unit area on a surface element in the ith direction may
be computed directly from the flow field around the body by:
fHi = τijnj =
[
−pδij + 2Re
−1
ω eij
]
nj (3.17)
where
eij =
1
2
(
∂ui
∂xj
+
∂uj
∂xi
)
(3.18)
and nj is the direction cosine of the unit normal vector in the xj direction. Equation
(3.17) has been non-dimensionalized using the same reference variables introduced
in the non-dimensionalization of the continuity (2.23), momentum (2.24) and energy
(2.26) equations, i.e. inertial scaling for the pressure ρL2f 2, fL for the reference
velocity and f for the reference time. The use of (3.17) requires the knowledge of p
and ∂ui/∂xj on the body surface. However, the MLS formulation outlined in §3.2.1
underestimates the actual traction forces due to the smoothness introduced in the
velocity field. Therefore, to find the pressure on a Lagrangian marker l and assuming
62
a linear velocity variation near the immersed bodies, ∂p/∂n is first evaluated from [81]:
∂p
∂n
= −
Du
Dt
· n (3.19)
where n is the unit normal vector passing through the Lagrangian marker l and
Du/Dt is the absolute linear acceleration of the marker l. The value of the pressure
at the surface can then be obtained from:
pl = pe −
∂p
∂n
hn (3.20)
where pe is the pressure at point e in the Eulerian domain, which lies on the unit
normal n at distance hn, where hn is the average grid spacing at point e.
To find the velocity gradients ∂ui/∂xj, the same MLS reconstruction discussed
in §3.2.1 is used, but instead of reconstructing (interpolating) the velocity compo-
nents, the gradients themselves are reconstructed using derivatives of the shape func-
tions. The velocity gradients on the immersed bodies can then be computed from:
∂Ui
∂xj
=
ne∑
k=1
∂φk
∂xj
ui (3.21)
where ne is the number of Eulerian grid points in the interpolation stencil and ∂φk/∂xj
results from the solution of a minimization problem similar to (3.13). Assuming a
linear variation of the velocity near the body, the derivatives, ∂Ui/∂xj, coming from
equation (3.21) are a good approximation of the derivatives ∂ui/∂xj at the surface.
The evaluation of the surface forces is generally more computationally expen-
sive per time step than the evaluation of the forcing functions because it requires the
reconstruction of the 9 components of the velocity gradient tensor while the com-
putation of the forcing functions requires only the reconstruction of the 3 velocity
63
components. Fortunately, the surface forces don’t need to be evaluated at each time
step in the current framework of prescribed kinematics.
3.3 Computational Setup
The numerical tools are almost ready to perform the simulations. To recap,
both the continuity equation and the Navier-Stokes equations (3.2) are solved on a
staggered grid where the advective and diffusive terms are evaluated using second
order differencing and time is advanced using a fractional step method, namely the
low-storage Runge-Kutta scheme (3.3). An IB method is used to implement the
moving boundary conditions, which use a MLS reconstruction technique to trans-
fer the forcing between the Eulerian and Lagrangian meshes. The numerical solver
was parallelized using a domain decomposition strategy [8]. The code was compiled
and executed on the University of Maryland High Performance Computing Cluster
(Deepthought). The code was developed in FORTRAN90, uses the MPIF90 library
and is executed on Intel’s Xeon processors with Linux RedHat. Details about the
computational setup used to simulate the mayfly nymphs follows.
A surface mesh of the mayfly computational model introduced in §2.1 was cre-
ated using Gambit (Gambit 2.4.6). Triangular elements were clustered near the gill
roots, where the surface geometry is most complex. The mesh was checked for quality
(skewness, regularity and inverted elements). In the end, the surface of the mayfly
abdomen, thorax and head required 659216 triangular elements (329610 nodes) to de-
scribe it. Coarser meshes resulted in highly skewed elements near the surface ridges
64
or even elements that don’t describe the surface correctly. On the other hand, only
118580 elements were needed to describe the 14 proximal plates and 79072 elements
to describe the 14 distal plates (planform A) making the total number of elements
required to fully describe the immersed boundaries 856868 elements.
The Lagrangian grid was then immersed into a computational domain, which
is a rectangular prism with dimensions 10L× 30L× 20L in the x, y and z directions
respectively as illustrated in Figure 3.1. The mayfly nymph lies longitudinally along
the y-axis. Plane z = 0 is the median plane (plane of symmetry) of the mayfly with
the left row of gills on the positive side of the z-axis. The x-axis is pointing dorsally
(upwards) from the origin.
As for the domain boundary conditions, the bottom wall is modeled as a no-slip
wall, while the other five walls are modeled as slip walls. It is worth pointing out
that the computational domain is sufficiently large to minimize effects of the domain
boundaries on the inner flow field. In fact, the bounding box encompassing the six
pairs of moving gills represents less than 0.2% of the domain by volume. The domain
is discretized with 200 × 1200 × 400 points in x, y and z directions respectively (96
million points total). Stretching is applied to cluster points in the region around
the mayfly. The cells size were designed so that 40 of them would span the entire
gill length when it is in an orthogonal plane. The code is typically executed using a
minimum of 100 processors. A typical sampling output yields approximately 300GB
of data saved for further post processing. Time discretization resulted in 562 time
steps per cycle for the highest Reynolds number (Reω = 21.6) and 10091 time steps
per cycle for the lowest Reynolds number (Reω = 1.0) making the computational
65
Figure 3.1: Schematic of computational Domain. (a,c) Computational domain lateral
and top view. (b) Closeup of mayfly posterior view (back) showing left gill 4. (d) Gill
immersed in computational domain. Numerical Grid shown is half the actual density
for clarity. (R)Root. (V)Ventral. (D)Dorsal. (M)Midchord. (T1)Tip.
time vary between 5 to 25 days for each simulation.
3.4 Comparison with Experiments
The second order accuracy of the NSE solver used in executing the simulations
was previously validated and led to a series of publications [8–10,80–82]. In addition,
the accuracy of the IB forcing using the MLS reconstruction as well as the calculus of
hydrodynamic forces was established by [68,69], where it was shown that it preserves
the second order accuracy of the numerical schemes. Therefore, in the current work,
the main purpose of the comparison with the experiments is to ascertain that the flow
66
physics observed in-vivo are replicated numerically. In other words, the numerical
grid is fine enough in order to capture the flow physics or more and the prescribed
kinematics are accurate enough to reproduce the main flow features of the rowing
and flapping.
Therefore, in the following section, comparisons with previous PIV experiments
performed on live mayfly nymphs is presented for both flapping at Reω = 21.6 and
rowing at Reω = 2.3. This will be followed by comparisons with robotic experiments
using an idealized rowing type of kinematics at Reω = 17.6. It will be shown that the
physical problem is replicated with sufficient accuracy that enables the execution of
other hypothetical situations without the need for further validation. The main goal
of the comparisons is to establish confidence in the numerical results. The specific
details of the flow fields however, will be deferred to §4.
3.4.1 Comparison with in-vivo Experiments
Previous experimental work by [57–59] investigated the rowing-flapping tran-
sition occurring in C. triangulifer. In order to investigate the beating kinematics,
9 nymphs were studied aging 20 − 42 days with 0.25mm < L < 0.77mm and
2.1 ± 0.18 < Reω < 21.6 ± 1.2. Oscillation frequencies were in the range 20 − 36Hz
with no significant trend associated with Reω. In reporting the results, the nymphs
beating kinematics were classified into three types: 1) rowing for Reω < 5, which was
characterized by a higher pitching range and a negligible hinge flexion angle (of the
same order of the digitization error), 2) flapping for Reω > 20, which was character-
67
ized by a minimal pitch variation and the formation of a distal flap that appeared
to move passively, and 3) intermediate kinematics for 5 < Reω < 20, which marked
abrupt and gradual changes in many parameters as for example the hinge flexion
angle, which increased from 17± 2◦ to 61± 7◦ and decreased stroke plane deviation.
After extracting the beating kinematics, PIV measurements were carried out.
The combination of the lenses and the Argon laser provided a light sheet that was
approximately 20mm wide and 0.15 − 0.2mm thick. The field of view varied from
3mm× 3mm for the smallest nymph and up to 9mm× 9mm for the largest nymph
(animal length varied between 2 − 5.64mm excluding the caudal filaments). The
vertical light sheet was parallel to the medial plane (vertical) and centered between the
lateral and medial extents of the gill plates (±0.2mm). Framing rates of 1kHz allowed
for cyclical resolution of approximately 10◦, with a spatial resolution of approximately
25µm for the 3mm×3mm field of view. Of the approximately 30 cycles recorded in a
typical imaging sequence, an ensemble ofN = 10 frame pairs was found to be sufficient
to converge velocity fields at the desired resolution. The maximum uncertainty for
the velocity was taken to be 2.5mm/s which was typically 5− 10% of the maximum
velocity.
Restricted by the availability of full and reliable experimental sets, comparisons
with two live nymphs were carried out. The first is for a large nymph flapping at
Reω = 21.6, while the second is for a small nymph rowing at Reω = 2.3. The first
comparison was very successful, while the second is less successful, probably due
to the higher uncertainties associated with smaller nymphs. For both comparisons,
the primary source of discrepancy stems from geometrical dissimilarities between
68
Figure 3.2: Comparison of the time-averaged in-plane velocity magnitude of the flap-
ping kinematics at Reω = 21.6. Velocity magnitude is normalized by the local maxi-
mum. Top: Simulation (showing 7 gills). Bottom: In-vivo experiments (showing only
6 gills). Non-moving gill pair 7 was not tracked by experiments. Contour floods below
0.1 are cutoff.
the dissected nymph used to build the computational model and the nymphs used in
extracting the kinematics and performing the PIV measurements. In other words, the
variability and the morphological changes inherent to the in-vivo experiments are the
main source of discrepancy. To give an example of such variability, one notices that
during the digitization of the flapping kinematics of the large nymphs, the beating
frequency was reported at 37Hz, while during the PIV analysis (performed on the
69
same day on the same nymph) showed a frequency of 34Hz. Another prevalent
source of discrepancy is the vertical distance between the nymphs and the bottom
wall. In PIV experiments, this distance was kept constant for all nymphs sizes,
which would not achieve geometrical similarity with the water tank since the nymph
keeps growing. In numerical simulations, the viscosity of the fluid was manipulated
rather than changing the gill size, which maintains geometrical similarity. This is
more desirable since it provides a self-consistent set of numerical results that can be
analyzed later. One more idealization was introduced to the computational model by
maintaining a constant gill length for all 7 pairs. In reality, the mayfly gills get smaller
towards the thorax and the tail making gill 1 and 7 the smallest while gill 3 and 4
the biggest. The actual size of gills 1 and 6 could be smaller than the computational
model by as much as 10%. The uncertainty in Reω (≈ ±1) is another contributing
factor that may affect the size of the different structures in all comparisons.
All PIV results presented here were obtained at the mid-section of the gills, i.e.
at z/L = 1.15. Additionally, to mimic the experimental resolution, the numerical
results were spatially averaged over a thickness equal to that of the laser sheet, which
was approximately 0.2mm or the equivalent of 8 numerical cells. Detailed comparisons
of the time-averaged flow fields of the flapping kinematics at Reω = 21.6 are shown in
Figures 3.2, 3.3 and 3.4. Figure 3.2 shows the time-averaged velocity magnitude nor-
malized by its local maximum. Overall, the qualitative agreement is acceptable and
all the main flow features are reproduced by the computations. As in the experiments,
the flow enters from the bottom and left planes (ventral and anterior), while the top
and right planes (dorsal and posterior) exhibit both inflow and outflow. Also, the
70
Figure 3.3: Comparison of the time-averaged x-component of the velocity u for flap-
ping kinematics at Reω = 21.6. Contour floods below 0.1 are cutoff.
predicted mean outflow direction matches that of the experiments. Quantitatively,
the predicted time-averaged velocity magnitude has the correct range and compara-
ble distribution. In the simulations, however, more localized peaks in the velocity
magnitude can be observed at the gills dorsal and ventral edges. This could be due
to limited experimental resolution or geometrical dissimilarities (inter-gill spacing is
smaller in experiments). It should also be noted that the differences in some of the
mean gill locations (i.e. gills 2 and 3), are due to the fact that the experimental gills
71
shown represent an approximate location of the gills (average of intersections with
laser sheet), while in the simulations, the exact mean location is shown. Noteworthy
is that the simulation shows that the out-of-plane velocity component (also the w
or lateral component) is localized around the distal plates and constitutes no more
than 15% of the velocity magnitude at its highest location. In other words, including
the lateral velocity component in the velocity magnitude calculation does not have a
remarkable qualitative effect on the comparison.
Figure 3.4: Comparison of the time-averaged y-component of the velocity v for flap-
ping kinematics at Reω = 21.6.Contour floods above 0.0 are cutoff.
72
Figures 3.3 and 3.4 show the individual time-averaged velocity components u
and v respectively. Again, the agreement is quiet remarkable given the variability
of the nymphs. The main source of discrepancy is the mean location of the gills.
Specifically, gills towards the back of the mayfly (right) are positioned higher in the
experiments. This could be due to the mayfly nymphs not oriented horizontally
during the experiments and/or a different curvature of the abdomen affecting the
mean location of the gills.
For further validation, instantaneous flow fields were also looked at. Figures 3.5
and 3.6 show comparisons of instantaneous vorticity plots at the end of retraction
of gill 4 and gill 3 respectively. In both figures, contour floods of the z-component
of the vorticity ωz, as well as the instantaneous velocity vectors are shown. The
instantaneous position of the gills is shown in both experiments and simulations.
Additionally, for the simulations, the extreme positions of the flapping gills are shown
in lighter shades.
For both figures, all structures are replicated and the instantaneous velocity
vectors are in good agreement, most importantly on top of most gills and in the inter-
gill spacings. In Figure 3.5, the main difference is the negative vorticity on top of gills
4−6 where in the PIV experiments, the three vortices seems to have started to diffuse
together, whereas in the simulations, the three vortices can still be distinguished. The
same difference appears also in Figure 3.6 with respect to gills 3− 5. Therefore, the
numerical results appears less dissipative, with stronger vorticity in proximity of the
gills and higher vorticity magnitudes at the vortex cores. This could be attributed to
the zero-thickness assumption of the gill plates, to the introduced rigidity of the gills
73
Figure 3.5: Comparison of the instantaneous z-component of the vorticity vector ωz
of the flapping kinematics at Reω = 21.6 and t/τ = 0.102 (end of retraction of gill
number 4).
74
Figure 3.6: Comparison of the instantaneous z-component of the vorticity vector ωz
of the flapping kinematics at Reω = 21.6 and t/τ = 0.326 (end of retraction of gill
number 3).
75
or to the approximation of the hinge, where the in-gill curvature was neglected and
the hinge was approximated as a straight line. As for the location of the gills, gills
2− 4 in the simulation are lagging their counterparts in the experiment compared to
gills 5 and 6, which are well matched. The time lag is, however, less than 5%. This
could be explained by the fact that the gill kinematics in the PIV experiments were
recorded and extracted at a different time than the PIV imaging. In other words, gill
kinematics and PIV were not performed simultaneously. Matching the exact timing
between the simulations and the experiments proved to be challenging due to the
many bodies involved in the flow field. In other words, attempting to match the flow
around gills 4 − 6 for example, will cause a slight time lag around gill 1 − 3 or vice
versa, which again highlights the variability in the experiments.
Overall, the simulation of a flapping nymph at Reω = 21.6 is able to reproduce
the flow physics and the different vortical structures associated with it. The time-
averaged flow fields show excellent agreement, while the instantaneous flow fields show
acceptable agreement. Most of the discrepancies can be explained and are attributed
mainly to insect variability and uncertainties associated with the in-vivo experiments.
This comparison also shows that the grid resolution is sufficient enough to capture
the flow physics, which in turn enables running flapping simulations at lower Reω
provided that the grid resolution is maintained the same.
When it came to comparing the rowing kinematics, there wasn’t a full experi-
mental set to work with because none of the rowing nymphs overlapped between the
kinematics extraction and the PIV imaging. Therefore, a nymph other than the one
used to extract the kinematics had to be selected to perform the flow field compari-
76
son. The matching was done based solely on the Reω and qualitative resemblance in
the rowing kinematics. In other words, 3 nymphs are involved in the rowing compar-
ison; the nymph used in building the computational model and the two nymphs used
by the experiments in extracting the kinematics and performing the PIV imaging.
However, all these 3 nymphs are not geometrically similar in the strict sense. This is
most noticeable in the ratio between the gill length and the inter-gill spacing.
The above dissimilarities are compounded by the higher uncertainty associated
with the smaller nymphs. The largest source of uncertainty may be attributed to the
location and orientation of the laser sheet with respect to the gills. The laser sheet
thickness was half that of the gill length (compared to one fifth of the gill length in
the flapping case) with a conservative uncertainty in its position of 0.1mm (half the
gill length). This can actually result in part of the laser sheet not cutting through
the gills. Clearly a high level of discrepancy is expected.
The comparison of the time-averaged velocity magnitude for rowing at Reω =
2.3 is shown in figure 3.7. The flow field far from the mayfly agrees well with the
experiments. The main flow features are preserved; that is the inflow and outflow
locations and magnitudes and the inflow direction. Furthermore, two recirculation
zones (bottom left and top right) are present in the flow field and they both reside in
regions of low velocity magnitude despite their different location. However, near and
around the gills, discrepancies are observed. The most important source of discrep-
ancy in this comparison is the location of the ventral wall (modeled as a no-slip wall)
with respect to the nymph. In order to match the mean location of the gills and the
inflow direction (top left), the experimental results had to be rotated −15◦ around the
77
z-axis (even though it wouldn’t have a similar effect as modeling the nymphs with a
different orientation). Additionally, the ventral wall in the simulation is closer to the
mayfly when compared to the experiment. This has a combined effect of bending the
flow leaving the mayfly at the bottom right away from the ventral wall. Noteworthy is
that the recirculation zone at the bottom left disappears at higher Reynolds numbers
as observed in simulations (shown later in §4.4), which means that the ventral wall
has a major influence on the results (at least on the ventral side of the flow field below
the abdomen) and therefore providing a reasonable explanation of the discrepancy.
Another feature that is shown by the simulation is the ’wiggly’ shape of the
streamlines in the inter-gill spacing. This could be caused by a larger inter-gill spacing
in the simulations or a low grid resolution on the experiments. Noteworthy is that
the numerical grid resolution in this case is about 4 times that of the experimental
grid. I tried to match the experimental results by changing the location where the
spatial averaging is performed or changing the depth over which the spatial averaging
is performed while remaining parallel to the medial plane (8, 16, 24, 32 cells). The
velocity vectors matched the experiments better when the averaging location moved
towards the gill tips, where the ’wiggly’ effect started to disappear. However the
orientation of the zone of highest velocity magnitude could not be matched.
As with the flapping at Reω = 21.6, I looked into instantaneous flow fields for
the rowing case at Reω = 2.3. A comparison of the instantaneous normalized velocity
magnitude is shown in Figure 3.8. The instant shown is at t/τ ≈ 0.075, which is the
end of retraction of gill 5 (see Figure 2.10 on Page 46 for rowing kinematics), when
the velocity magnitude near the dorsal edge of gill 5 is at its maximum. Given the
78
Figure 3.7: Comparison of the time-averaged in-plane normalized velocity magnitude
for the rowing kinematics at Reω = 2.3. Velocity magnitude is normalized by the
local maximum. Top: Simulation (showing 7 gills). Bottom: In-vivo experiments
(showing only 5 gills). Contour floods below 0.1 are cutoff. The sink-like behavior
of the streamlines in the simulation is because of the forced 2D integration of the
velocity field, which also implies a normal component of the velocity (w).
79
previous discussion, the agreement is still reasonable around most gills.
Overall, the comparison is less remarkable compared to the flapping case at
Reω = 21.6 but still acceptable for the purpose of verifying the correct implementa-
tion of the kinematics. Outflow direction differs by approximately 90◦ compared to
the flapping kinematics, which was shown by experiments and now substantiated by
computations.
To sum up, comparisons with available in-vivo experimental data are suggesting
a correct implementation of both kinematics and that the grid resolution is sufficient
to capture all the flow physics and the general behavior of flapping and rowing mayfly
nymphs. In fact, the numerical grid is finer than required when simulating the lower
Reω cases.
3.4.2 Comparison with Robotic Experiments
To further validate my numerical code, comparisons were made with mayfly
robotic experiments carried out by [43,44]. The robotic mayfly experiment consisted
of an array of 5 gills, which were controlled using programmed micro-servomotors.
The robot used a simplified version of the rowing kinematics previously introduced
in §2.2.3, where the stroke and pitch angles were prescribed, but the stroke plane
deviation angle was neglected. Additionally, for all 5 gills, the stoke plane lateral
inclination angle  was set to 90◦, while the stroke plane inclination angle β was
−60◦. For all 5 gills, stroke and pitch angles were based on the kinematics of rowing
gill number 4 (see Figure 2.10) but varied slightly from each other due to limitations
80
Figure 3.8: Comparison of the instantaneous normalized velocity magnitude for the
rowing kinematics at Reω = 2.3 and tτ ≈ 0.075. Velocity magnitude is normalized by
the local maximum. Top: Simulation (showing 7 gills) spatially averaged over 0.2L.
Bottom: In-vivo experiments (showing only 5 gills). Contour floods below 0.1 are
cutoff.
81
Figure 3.9: Illustration of computational setup and slice at z/L = 1.3
of the servomotors.
The gill plate array had a length of 96mm (from gill 1 to gill 5) and was immersed
in a small aquarium (260× 310× 510mm) partially filled with mineral oil (depth of
fluid = 200mm), and oriented such that the bilateral symmetry plane of the mayfly
body corresponded with the free surface of the liquid. A simplified body shell for
the mayfly was constructed, which only models the abdomen of the mayfly nymph
to which the gills are attached. The model was scaled using the oscillating Reynolds
82
number Reω. The robotic model dimensions are approximately 54 times those of the
original mayfly (L = 40mm). Mineral oil with a viscosity of 175cSt was used as the
working fluid in order to achieve a Reω of 17.4. More details about the setup can be
found in Appendix A.
In order to simulate the robotic experiment, I constructed another Eulerian
grid with dimensions 6.4L× 10L× 12.5L. The same grid density introduced in §3.3
was used, i.e. 40 cells per unit length. Consequently, the rectangular domain was
discretized with 258 × 242 × 194 points in the x-, y- and z-directions, which yielded
a total of 12 million points. Grid stretching was applied in the x- and z-directions
in order to cluster points in the region around the immersed bodies. All domain
boundaries are modeled as non-slip walls. An illustration of the computational setup
is shown in Figure 3.9.
Figure 3.10 shows a comparison of the time-averaged x-component of the vor-
ticity vector ωx. The results are shown for a medial plane z/L = 1.3. The location of
the slice is illustrated in Figure 3.9. The comparison shows excellent agreement of the
vorticity fields, which confirms that a grid density of 40 cells per unit length is suffi-
cient to replicate the flow physics. The comparison is much better than those with the
in-vivo experiments since the environment is more controllable and uncertainties are
much lower than those associated with the in-vivo experiments. More comparisons of
the instantaneous flow fields as well as integral quantities are available in Appendix
A.
This section concludes the comparisons that were performed in order to validate
the implementation of the kinematic chain, the two sets of kinematics that will be
83
Figure 3.10: Comparison of time-averaged flow fields for simplified rowing at Reω =
17.6. Contour floods of ωx at z/L = 1.3. Top: simulation. Bottom: Robotic experi-
ments.
used to carry out the simulations as well as the grid density. With that, the road
is paved for numerical experiments that deals with hypothetical scenarios, such as
rowing at higher Reω and flapping at low Reω. In the other words, the parametric
study of the rowing-flapping transition can be executed. In the next section, I will
introduce the two parametric spaces that were built throughout this study. This will
be followed by the results in §4 and §5.
84
3.5 Parametric Spaces
The primary goal of this study was to investigate the rowing-flapping transition
of the mayfly nymph C. triangulifer, which was first observed by the PIV experiments.
The PIV experiments had reported three types of kinematics throughout the ontogeny
of the nymph; rowing, intermediate and flapping kinematics. Therefore 3 different
Reω were selected so that each lies in one of 3 kinematics regimes revealed by the
in-vivo experiments. These were 2.3, 8.1 and 21.6 belonging to rowing, intermediate
and flapping regimes respectively. For each of these 3 Reω, 2 simulations were carried
out, one using rowing kinematics and the other using the flapping kinematics, making
the total number of simulations 6.
Initially, all simulations were performed using Planform A, which showed some
potential to explain the transition but required running two more simulations at a
lower Reω of 1.0 bringing the total number of simulations to 8. Later, the planform
shape was found to affect some of the results. Therefore, I decided to study the effect
of the gill planform when using rowing kinematics, thereby introducing planforms B
and C. This was probably driven by the higher uncertainties associated with rowing
nymphs which actually use planforms similar to planforms B and C. Therefore, both
planform B and C were used to carry out rowing simulations at the the same 4 Reω,
now 1.0, 2.3, 8.1 and 21.6. This brought up the final number of simulations to 16.
The rowing-flapping parametric space is summarized in Table 3.1 and the results will
be presented in the following chapter; §4. Noteworthy is that optimization techniques
would have been ineffective in this study because of the computational cost associated
85
Reω Ref Flapping Rowing Nature favors Sample time Itr/cycle
1.0 0.26-0.82 A A,B,C Rowing 10.0 10000-10091
2.3 0.60-1.89 A A,B,C Rowing 13.0 4386-4505
8.1 2.11-6.64 A A,B,C Transition 16.0 1333-1429
21.6 5.62-17.71 A A,B,C Flapping 50.0 562-682
Table 3.1: Parametric space used to explore the rowing-flapping transition. 3
independent variables are investigated: 1) Reω {1.0,2.3,8.1,21.6}, 2) kinematics
{flapping,rowing} and 3) Planforms {A,B,C}. Reω = L2f/ν. Ref = fAL/ν where
A is the stroke amplitude (peak-to-peak). The range in Ref for some Reω is caused
by the subtle differences in gills kinematics and different stroke amplitudes between
flapping and rowing.
with the simulations and the relatively small number of simulations that were carried
out. An optimization approach would have required many more simulations in order
to build an initial model. More importantly, the optimized parameter using the
current mathematical model was not known a priori.
Interestingly, a preliminary run using flapping kinematics at Reω = 21.6 without
a hinge was carried out. This was in the testing phase of the kinematics chain, where
only two levels have been introduced. The results did not compare well with in-
vivo experiments and therefore the hinge or the two-plate model was introduced, i.e.
the third level of the kinematic chain previously discussed. This resulted in better
agreement with the experiments as shown in the previous section (§3.4.1). This also
proved the necessity of including the hinge when investigating the rowing-flapping
transition of the nymphs.
In light of the above, I decided to investigate the effect of introducing the hinge
as a first step towards investigating the flexibility of the gills. After all, flexibility
may be achieved in the limit of an infinite number of hinges or plates. Of course,
86
current experimental data is not sufficient to simulate the slight flexibility present in
the flapping gills. On the other hand, the numerical model would require substantial
development in order to simulate deforming Lagrangian meshes. With that, I investi-
gated the role of the hinge by running flapping simulations at Reω = {2.3, 21.6} using
planform A, with and without the flexural hinge (two-plate vs one-plate model). I
also decided to investigate the interaction between the different gills. Therefore, sim-
ulations were carried out using the entire gill array (7 pairs) and using only gill pair
number 4. This led to a parametric space that encompasses 8 simulations, where 2
of this simulations overlaps with rowing-flapping parametric space. In other words,
three independent variable are involved in the second parametric space. These are
the Reω, the flexural hinge (with or without) and the number of gill pairs present in
the array (1 or 7). The results of the second parametric space will be summarized in
§5.
87
Chapter 4: Results I: The Rowing-Flapping Transition
4.1 Introduction
The parametric space introduced in this chapter was designed to illuminate
the rowing-flapping transition by investigating the effect of 3 different parameters,
namely the Reω, the type of kinematics and the gill planform as summarized in
Table 3.1. Originally the goal of this parametric space was to simulate flapping at
low Reω and rowing at high Reω but was subsequently extended to include the effects
of planform morphology. Therefore, this chapter is organized as follows. The first 4
sections deal with mostly qualitative investigation of the different flow fields aiming
to provide better understanding of the kinematics and the morphology effects. In
specific, the time-averaged flow fields of rowing and flapping are shown in §4.2. This
will be followed by a description of the instantaneous fields in §4.3. The effects of
varying Reω is shown in §4.4, while gill planform effects are presented in §4.5. The last
section of this chapter provides a quantitative assessment of the micro-pumping using
the performance metrics introduced earlier, recaps over the qualitative differences
observed and provides a viable explanation of the rowing-flapping transition.
88
4.2 General Description of the Time-averaged Flow Fields
The time-averaged flow field for both flapping and rowing kinematics at Reynolds
numbers corresponding to the in-vivo experiments is shown in Figures 4.1a and 4.1b
respectively. The time-averaged flow field for the flapping kinematics is at Reω = 21.6,
while that for the rowing kinematics is at Reω = 2.3. Only part of the domain is
shown and the left row of gills is omitted for clarity. The sub-domain shown is one unit
length away from the gills and extends from x = −1.11 to x = 1.89, from y = −2.71
to y = 2.69 and from z = −2.73 to z = 2.73, which constitutes only 1.47% by volume
of the entire domain. Also, the right row of gills shown is time-averaged in order to
illustrate the average location of the Lagrangian grid. To highlight the differences
in the mean flow patterns between flapping and rowing, streamlines are shown. The
start and end of many streamlines are marked with a letter corresponding to the plane
it intersects at inlet and exit of the sub-domain. The streamlines are colored with
the velocity magnitude, which is higher near the gills and decays further away. On
the bounding planes, contours (lines or floods) of the velocity component normal to
each plane are shown. Dashed isolines indicate a negative velocity component with
respect to the inertial frame of reference (and not with respect to the plane normals).
The time-averaged velocity has been integrated across the six different planes of the
control surface shown in Figure 4.1. The control surface shown is one unit length (L)
away from the bounding box of the moving gills. This yields the inflow and outflow of
each plane and, in one way, highlights the difference between the two sets of prescribed
kinematics. For the flapping case (Figure 4.1a), the integration reveals that the dorsal
89
(a) Flapping at Reω = 21.6
(b) Rowing at Reω = 2.3
Figure 4.1: Time-averaged flow fields for rowing and flapping. w: velocity component
in the z direction. vm: velocity magnitude. Contour floods and lines: velocity
component normal to the corresponding plane. Streamlines are colored by velocity
magnitude. X,Y,Z denote intersection of the streamline when entering or exiting the
shown sub-domain. x = {−1.11 : 1.89}. y = {−2.71 : 2.69}. z = {−2.73 : 2.73}.
90
plane (x = xmax) contributes 60.27% of the total outflow while the posterior plane
(y = ymin) contributes 28.64%. In other words, 88.91% of the outflow goes through
a plane which is almost perpendicular to the apparent motion of the gills (flapping-
like). The remaining outflow exists through the medial planes (z = zmin&z = zmax)
and through the ventral plane (x = xmax) at 10.94% and 0.15% respectively. As
for the inflow, the flow comes into the control surface from the dorsal, anterior, and
medial planes at 22.26%, 23.65% and 45.49% respectively. For the rowing case, on
the other hand (Figure 4.1b), the primary inflow is through the dorsal plane which
contributes approximately 84.93% of the inflow. The anterior plane contributes only
7.31%, and all other planes have negligible inflow. The outflow contribution through
the medial, posterior and ventral planes is 59.81%, 31.03% and 7.16% respectively.
The mean outflow direction for the rowing kinematics is approximately perpendicular
to the mean outflow for the flapping kinematics. The contribution of each plane to
Kinematics Flapping Rowing
Figure 4.1a 4.1b
Mass flow rate Inflow Outflow Inflow Outflow
Plane % % % %
Ventral (bottom) 4.44 0.15 1.72 7.16
Dorsal (top) 22.26 60.27 84.93 1.99
Posterior (back) 4.19 28.64 2.12 31.03
Anterior (front) 23.65 0.00 7.31 0.00
Medial (left) 22.75 5.47 1.96 29.93
Medial (right) 22.75 5.47 1.96 29.93
Table 4.1: Table showing the dimensionless mass flow rate across the control surfaces
shown in Figure 4.1 for flapping at Reω = 21.6 and rowing at Reω = 2.3. Boldface
notation highlights the major exiting planes and hence the primary difference between
the flapping kinematics and the rowing kinematics.
91
the inflow and outflow is summarized in table 4.1. The above results are in general
agreement with what has been observed in experiments and illustrates the difference
between the flapping kinematics and the rowing kinematics [57–59].
4.3 Dynamics of Rowing and Flapping
Before discussing the parametric spaces and the effect of Reω and gill planform
on the flow field and performance, the micro-pumping dynamics of both flapping
and rowing is introduced in order to elucidate the complicated nature of the flow
field. This is achieved by looking into the instantaneous flow fields, which will also
demonstrate the metachronal wave and consequently highlight the main differences
between both kinematics.
Recall that the metachronal wave of the gills starts posteriorly (tail side), i.e. at
gill 6. Gill 6 starts its retraction, followed by gill 5 then gill 4 and all the way to gill 1.
Also, the phase lag between any two consecutive gills is approximately one quarter of
the period. Since all 6 pairs of moving gills have similar kinematics (beating, pitching
and stroke plane deviation), it is sufficient to follow any of the 6 moving gills in order
to elucidate the transient flow field. In addition, the flow field is symmetric about the
medial plane (z = 0), which means that looking at one side of the mayfly is sufficient
to understand the nature of the flow field. In light of that, the flow field around gill
4 is discussed hereunder.
92
4.3.1 Flapping Dynamics at Reω = 21.6
The pumping sequence is described in Figure 4.2 through examination of four
points in the flapping cycle, covering key features of the gill’s interactions with its
surrounding neighbors. To provide a sense of the beating direction of the gills, their
location at an earlier time step (0.15 t/τ earlier) is shown. As gill 4 protraction
initiates its motion towards the head (t/τ = 0.102), gill 5 immediately posterior
(behind) gill 4 is already near its peak protraction speed, pushing fluid both dorsally
and ventrally. At t/τ = 0.476, gill 4 is nearing the end of its protraction phase and
slowing, while gill 5 initiated retraction and moving toward the tail. This relative
motion causes an increasing interstitial volume (volume between adjacent gills) with
a strong negative suction pressure, sucking fluid into the widening volume. Note
that the biased nature of the hinge (effectively flexing in only one direction), allows
the interstitial volume to be blocked from the space dorsal (above) the gill array,
forcing the fluid to be entrained from the ventral region. Note that this filling process
produces vorticity on the gill surface similar to what would be observed in entrance
region of a planar Poiseuille flow. By t/τ = 0.802, gill 4 is near its peak retraction
speed, while gill 5 has finished retraction is about to reverse direction and start a
protraction stroke. This causes a compression of the interstitial volume, and an
increase in the interstitial pressure (pressure between adjacent gills) . The flexure
on gill 4 allows the distal plate to open, and the fluid is ejected out the dorsal side
of the array, carried forward by the inertial bias from the filling event. Lastly, at
t/τ = 0.976, gill 4 is just past its peak retraction speed, while gill 5 has started its
93
protraction, producing strong relative closing speeds between the gill tips. Counter-
rotating vortices are generated on each adjacent tips (positive vorticity on gill 4,
negative vorticity on gill 5), forming a dipole pair that results in a strong ejection
jet directed upward and behind (dorsal-posterior) the gill array. This event coincides
with peak velocity magnitudes and pressure. The close interaction of the gills and the
formation of this dipole pair was identified as a key feature of the flapping array in
prior work [59]. Due to the relative phasing of the gills, this process is repeated with
each gill, such that the filling and ejection cycle is observed to move in a metachronal
wave proceeding from the tail (posterior) end toward the head (anterior), with a
spatial wavelength corresponding to 5 gills (4 inter-gill spacings).
The above description of the instantaneous flow field was conventionally pre-
sented by moving forward in time and looking at 4 different snapshots of the transient
field. However, another way of presenting the instantaneous flow field my be achieved
by looking at only one snapshot and riding spatially against the metachronal wave.
This is possible because the metachronal wave has a spatial wave length equal to 4
inter-gill spacings (distance between 5 consecutive gills). Had the metachronal wave
length been longer than the entire gill array, the following method of illustration may
have not been possible. Therefore, the pumping sequence may be described if one
follows zones A, B, C then D shown in Figure 4.3. Following these zones is the same
as following one inter-gill spacing through time (illustrated through Figure 4.2). In-
side any of these zones, I refer to the gill at its left as the lagging or anterior gill and
to the gill at its right as the leading or posterior gill. Zone A: the proximal plate of
the lagging gill is starting to protract, while the leading gill is protracting already.
94
Figure 4.2: Instantaneous flapping dynamics at Reω = 21.6. Slice at z/L = 1.15
spatially averaged over 8 cells (z = {1.07 : 1.26}). Focus is on inter-gill spacing 4-5.
(Col a) Vorticity component (ωz). (Col b) Pressure (p). (Col c) Velocity magnitude
(vm). (Row 1) t/τ = 0.102: start of protraction (G4). (Row 2) t/τ = 0.476: mid-
protraction. (Row 3) t = 0.802/τ : mid-retraction. (Row 4) t/τ = 0.976: approaching
end of retraction.
95
Zone A is marked by lower pressure, low velocity magnitude, no vorticity production
at the ventral edge of the lagging gill and positive vorticity generation at the ventral
edge of the leading gill. Zone B: the lagging gill is protracting while the leading gill is
about to retract. This yields to an effective widening of the inter-gill spacing which
causes the fluid particles to accelerate from the ventral side of the mayfly towards the
inter-gill spacing. Zone B is marked by the lowest pressures, higher velocity magni-
tude and the formation of a dipole on the ventral edges of the leading and lagging
gills. Noteworthy is that the position of the distal plate of the lagging gill is oriented
in such way as to prevent a retrograde flow from the dorsal side towards the inter-gill
spacing. Zone C: the lagging gill is about to retract while the leading gill is already
retracting. The inter-gill spacing seems to maintain a constant width, where the
fluid particles keep moving towards the dorsal side and pushing the distal flap of the
lagging gill open. Vorticity is observed along the entire gill surface instead of being
focused on the ventral edges and the pressure gradients are at their lowest. Zone D:
the lagging gill is retracting while the leading gill is finishing its retraction and about
to protract. The effective narrowing of the inter-gill spacing accompanied by the weak
inertia of the fluid particles and the dipole forming on the dorsal edges squeezes the
fluid particles out to the dorsal side. The above sequence can be observed in any
of the 5 inter-gill spacings if one were to follow a specific inter-gill spacing through
time. Each inter-gill spacing produces more, or less current, depending on the subtle
variations in the kinematics, mean position of the gills and mean inter-gill spacing.
One may summarize the inter-gill micro-pumping sequence as follows. The se-
quence starts when the inter-gill spacing is widened which creates a pressure drop
96
Figure 4.3: Instantaneous fields at t/τ = 0.65 for flapping kinematics at Reω = 21.6.
Slice at z/L = 1.15. Gills intersection with slice is shown. Lighter gill shades are
from an earlier time step (0.15 t/τ earlier). (a) z-component of the vorticity vector.
(b) Velocity magnitude. (c) Pressure. (a-c)Velocity vectors shown are in-plane. G1-
G7: Gills 1 through 7. A-D: Pumping sequence for main stream II. E-F: Vortex
detachment from gills 2 and 3 respectively assisting main stream I. Percentages shown
are for mass inflow or outflow.
97
inside the inter-gill spacing, pulling the flow from the ventral side. The flow is en-
hanced by the creation of a dipole at the ventral edge (Figure 4.2, row 2 or Figure 4.3,
zone B). During the “suction” phase, the distal flap reduces the amount of retrograde
flow into the inter-gill spacing from the dorsal side. The “suction” phase is followed
by a “squeeze” or “compression” phase, where the flow is pushed dorsally and poste-
riorly with a slight retrograde ventrally (Figure 4.2, row 4 or Figure 4.3, zone D). The
postero-dorsal flow is enhanced by the creation of a dipole attached to the distal edge
of the gill. Interestingly, Eastham had predicted the periods of suction followed by
periods of compression in many of his publications about mayfly nymph without ac-
tually measuring the pressure and by looking only at velocity fields [15–18]. Eastham
prediction is now confirmed by the full numerical simulation presented in this work.
Also, Sensenig [57] hypothesized that the distal flap could function as the valve of a
pump. Based on the numerical results, the distal flap is acting similar to a non-return
valve by reducing the retrograde flow during the “suction” phase. It is impossible to
prevent all the retrograde flow due to the shape of the distal plate and the fact that
the retrograde flow occurs also laterally.
The periodic widening and narrowing between consecutive gills highlights the
first element of synergy between flapping gills, and appears to play an important role
in fluid circulation. This type of motion, however, primarily accounts for that part of
flow that is directed from the ventral side of the mayfly to the dorsal side (marked by
zone II in Figure 4.3), which is relatively small (5% of inflow) as was demonstrated
previously by table 4.1. While the inter-gill periodic pumping mechanism is inducing
a small percentage of the mass flow rate, its role cannot be overlooked. The periodic
98
suction and compression of the fluid may have its effect on oxygen transport and gas
exchange between the fluid and the tracheal gills (oxygen dissociation).
It turns out that the major contribution of the pumping is induced longitu-
dinally along the dorsal side of the nymph’s body, where fluid enters dorsally and
anteriorly at 22.26% and 23.65% respectively, flows dorsally along the mayfly, then
exits dorsally and posteriorly at 60.27% and 28.64% respectively (see table 4.1). This
flow is primarily induced by the antiplectic metachronal wave that starts form the
back to the animal. To that effect, the flexing of the distal plate during protraction
appears to reduce retrograde flow in the anterior direction. During retraction, the
entire gill is pushing fluid posteriorly, while during protraction, the flexing of the
distal plate reduces the frontal area of the gill thereby reducing retrograde flow in
the anterior direction. In order to illustrate this second element of synergy, namely
the antiplectic metachronal wave, pathlines were integrated over 10 flapping periods
(after reaching the steady oscillatory state) throughout the near field as shown in
Figure 4.4. Seven groups of pathlines are shown, all of which are for particles re-
leased near z/L = 1.15, which is the medial plane cutting through the middle of the
gills. Of particular interest is group ’f’, where particles were released at the dorsal
edge of gill 3. At that instant, gill 3 was retracting and therefore doing work on
those particles and pushing them posteriorly towards gill 4. At Reω = 21.6, it took
only one beating cycle for these particles to reach gill 4, which at this point is also
retracting and therefore pushing the particles further towards gill 5, which in turn
helps ejecting those particles towards the mean flow direction. The above synergy
was found to occur between many consecutive gills. One more example is shown by
99
Figure 4.4: Pathlines integration for flapping kinematics over 10 periods at Reω =
21.6. All particles were released near z/L = 1.15. vm: velocity magnitude. z: depth.
Groups a-d: particles released on the edges of the integration domain, i.e. anterior,
dorsal, posterior and ventral planes respectively. Group e: particles released in inter-
gill spacing 1 − 2. Group f: particles released at the dorsal edge of gill 3. Group g:
particles released in inter-gill spacing 5 − 6. Groups e-g: squares highlight one time
unit. Contour floods: time-averaged velocity magnitude. Gills are shown at their
mean position. Pathlines are colored by depth (z).
following group ’e’ where synergy between gill 1 and 2 is also evident.
So far, the flapping kinematics have demonstrated two elements of synergy.
The first element is the periodic widening and narrowing occurring between any two
consecutive gills which leads to a flow directed from the ventral side of the mayfly to its
dorsal side. This flow is almost perpendicular to the mayfly axis with a slight posterior
tilt depending mostly on the mean position of the corresponding gill. Currently, the
contribution of this flow is at 5% of the inflow, but it is speculated to increase if
the distance between the mayfly nymph and the ventral wall increases. The second
element of synergy is created by the collective effect of the antiplectic metachronal
wave which enables almost the entire gill array to work together by doing work on the
100
same particle as it flows horizontally along the mayfly naiad, which was illustrated in
Figure 4.4. Later, when the effect of Reω is presented, it will be shown that increasing
Reω is actually enhancing both elements of synergies.
Noteworthy is that gills towards the nymph head are tucked closer to the ab-
domen while those towards the back are extended more laterally with the proximal
plates almost parallel to the transverse planes (xz planes). In specific, for gill 1,
 = 97.17◦, while for gill 2,  = 91.37◦. Also, actual nymphs have smaller gills to-
ward the head and the tail which was not incorporated in the numerical model. It
looks that the gradual shift in the mean position of the gills is designed so that gills
toward the head are more streamlined with the oncoming flow while those towards
the back are designed to redirect the outflow in a postero-dorsal direction which is
evident in the time-averaged flow field in Figure 4.1a. In addition, table 4.1 shows
that the dorsal and posterior planes contribute to 88.91% of the outflow while the
anterior, medial and dorsal planes contribute with 91.41% of the inflow with all three
planes showing equal importance for the inflow. These flow ratios complement the
speculation about the function of the shift in the gills mean position.
4.3.2 Rowing Dynamics at Reω = 21.6
To contrast the flapping dynamics, a similar flow sequence is shown for the
rowing conditions (also at Reω = 21.6) in Figure 4.5. While in general there are
similarities with flapping kinematics, the greatest differences to be observed are gen-
erated by the asymmetry in the protraction/retraction duration, the pitch range and
101
Figure 4.5: Instantaneous rowing dynamics at Reω = 21.6. Slice at z/L = 1.15
spatially averaged over 8 cells (z = {1.07 : 1.26}). Focus is on gill 4. (Col a) Vorticity
component (ωz). (Col b) Pressure (p). (Col c) Velocity magnitude (vm). (Row 1)
t/τ = 0.100: end of recovery stroke (left), begin power stroke (right) (ψ4 = 1.43◦, φ4 =
−79.55◦). (Row 2) t/τ = 0.251: power stroke, maximum pitch (ψ4 = 13.06◦, φ4 =
−99.85◦). (Row 3) t = 0.500/τ : ending power stroke (ψ4 = 48.31◦, φ4 = −78.04◦).
(Row 4) t/τ = 0.752: recovery stroke, minimum pitch (ψ4 = 33.01◦, φ4 = −63.16◦).
102
the absence of any gill flexion. This, in effect, isolates the gill during its quicker power
stroke during retraction, while permitting collective synchronized protraction during
the lengthened recovery stroke. In detail, retraction starts for gill 4 at t/τ = 0.100,
positioned at its maximum forward position, while gill 5 has already started its rear-
ward motion. Similar to the flapping kinematics, this creates a low pressure region
within the interstitial space (space between adjacent gills), drawing fluid up from the
lower ventral region. Unlike the flapping case, however, the absence of a flexure in
the gill plate also allows fluid fill from the region above (dorsal) the gills. As gill
4 reaches mid-retraction (t/τ = 0.251), gill 5 has completed its power stroke and is
at its maximum retracted position. Gill 4 is now moving largely in isolation, with a
vortex dipole formed on its dorsal (positive vorticity) and ventral (negative vortic-
ity) edges, which is really the signature of a bound vortex ring generated along the
periphery of the gill. This creates a dominant induced motion in the direction of
the advancing gill, drawing fluid downward and behind the animal. On the anterior
space (between gill 3 and 4), a low pressure is created by the expanding interstitial
volume, inducing flow up from below the animal and down from above. The pressure
in the compressing interstitial space between gill 4 and 5 is now positive, with the flow
moving being pushed both downward and outward (perpendicular to the plane of the
figure), as evidenced by the mean flow shown in Figure 4.1b and table 4.1. As gill 4
ends its retraction stroke (t/τ = 0.500, the compression phase begins on the anterior
side of the gill, reversing the pressure gradient along the anterior surface and also
reversing the sign of the vorticity production along the surface, leading to a reversal
of the vortex dipole bound to the edge. At t/τ = 0.752, gill 4 is in the middle of its
103
lengthy protraction stroke and is surrounded by two neighboring gills, which are also
completing the beginning and ending of the their respective protraction strokes.
During the qualitative analysis of the rowing flow fields, and in contrast with
flapping flow field, no form of synergy could be detected that would bring the entire
gill array to work together. As mentioned earlier, it seems that each gill performs
its quick power stroke in isolation, then slowly retracting while decreasing its pitch
during the recovery stroke in order to decrease the viscous drag.
4.4 Viscous and Inertial Effects
In reporting their results, [58] had classified the nymphs beating kinematics
into three types: 1) rowing for Reω < 5, which was characterized by a higher pitching
range and a negligible hinge flexion angle (of the same order of the digitization error),
2) flapping for Reω > 20, which was characterized by a minimal pitch variation and
the formation of a distal flap that appeared to move passively, and 3) intermediate
kinematics for 5 ≤ Reω ≤ 20, which marked abrupt and gradual changes in many
parameters such as a decreased stroke plane deviation and an increased hinge flexion
angle, which increased from 17± 2◦ to 61± 7◦.
In light of the experimental findings by [58], my initial parametric space was
selected to include 3 different Reynolds numbers that each lie in a different kine-
matical regime. Based on the available experimental data, the choice was limited to
insects corresponding to Reω = {2.3, 8.1, 21.6}. Following my preliminary findings,
the parametric space was subsequently expanded to include a lower Reω = 1.0. In
104
other words, in order to investigate the effect of Reω on the rowing-flapping transition,
4 simulations were carried out at Reω = {1.0, 2.3, 8.1, 21.6} for each set of kinematics
described in §2.2.
It has long been recognized that all sources of vorticity in homogenous fluids
must lie at the boundaries of fluid regions. Studies on vorticity show that its rate of
generation involves the tangential pressure gradient within a fluid and the external
acceleration of the solid boundary. It is also known that vorticity diffuses neither out
of boundaries nor into them, and that the only means of decay is by cross-diffusive
annihilation within the fluid [50], such as between two counter-rotating vortices. In
light of that, it is speculated that the flapping type of kinematics generate more
vorticity compared to the rowing type of kinematics. This is because of the presence
of the distal flap, which, during protraction is almost aligned with the mean flow and
hence providing more tangential acceleration to the surrounding fluid (see Figure 2.7
at t/τ = {0.2 : 0.5} for the kinematics and Figure 4.2a2 for the negative vorticity
generated). This is in contrast with rowing where most of acceleration is normal to
the mean flow and the vorticity generated is used primarily in creating the bound
vortex (Figure 4.5).
In order to illustrate the effects of Reω on the vorticity dynamics of both kine-
matics, contours of the z-component of vorticity (ωz) are shown in Figure 4.6 and
Figure 4.7 for the flapping and rowing kinematics respectively. The focus is on left
gill 4 by the end of its retraction (t/τ = 0.125 for flapping and t/τ = 0.475 for row-
ing). ωz has been normalized by the local maximum that occurs at the core of the
positive vortex near the dorsal edge of the gill. A qualitative criteria has been set to
105
identify vortex detachment. If the local maximum of ωz remains close (attached) to
the dorsal edge of the gill by the end of retraction, then it is concluded that the gen-
erated vorticity is quickly diffused or annihilated by its neighboring counter-rotating
negative vortex. In other words, the momentum diffusion is high enough that no
vortical ’structure’ could detach from the gill. On the contrary, if the local maximum
of ωz slightly convects away from the dorsal edge and lasts for at least 5% of the pe-
riod, then it is concluded that the bound vortex have detached from the gill and have
survived cross-diffusion. The attached or detached distinction can be compared to
flow behind cylindrical structures, where at low Re, two counter-rotating vortices or
recirculation bubbles remain attached to the cylinder, while at higher Re, instability
develops leading to the von Ka´rma´n street vortex and vortex detachment [60].
With reference to Figures 4.6 and 4.7, it is noted that the vortex detachment
occurs at Reω ≥ 8.1 for both types of kinematics. In both Figures parts (a) and (b),
the positive vortex has detached from gill 4, while in parts (c) and (d), the vortex
remains visibly attached and lose coherence quickly under viscous diffusion. Granted
the detached vortex core has smaller magnitude compared to the negative vorticity
being generated by the distal flap, it marks the onset of weak inertia. Also, this
detachment occurs more or less for all six gill pairs when Reω ≥ 8.1 and therefore its
relation to the rowing-flapping transition is unlikely to be coincidental.
The subtle detachment of any of the aforementioned vortices for Reω ≥ 8.1 has
a noticeable effect on the velocity vectors because of its slight interaction with the
negative vorticity being generated by the gill as it ends its retraction and initiates
its protraction. For both types of kinematics, the short-lived instantaneous dipole
106
Figure 4.6: Effect of Reω on ωz for flapping kinematics. Contours floods and lines
of ωz normalized by its local maximum (positive vortex core near the dorsal edge).
t/τ = 0.125: end of left gill 4 retraction. Orthogonal view of a z-slice at z/L = 1.15.
Vectors shown are in-plane velocity vectors. (a) Reω = 21.6. (b) Reω = 8.1. (c)
Reω = 2.3. (d) Reω = 1.0.
is expelling fluid in a horizontal direction (left to right), which acts in favor of the
flapping kinematics, but not so much in favor of the rowing kinematics. Figure 4.6
illustrates the effect on the velocity vectors in the right vicinity of the dipole, where
at low Reω (parts (c) and (d)), the vectors are pointing more upwards, but with
107
Figure 4.7: Effect of Reω on ωz for rowing kinematics. Contours floods and lines
of ωz normalized by its local maximum (positive vortex core near the dorsal edge).
t/τ = 0.475: end of left gill 4 retraction. Orthogonal view of a z-slice at z/L = 1.15.
Vectors shown are in-plane velocity vectors. (a) Reω = 21.6. (b) Reω = 8.1. (c)
Reω = 2.3. (d) Reω = 1.0.
increasing Reω, they start to point more horizontally (parts (a) and (b)). The velocity
vectors at Reω < 8.1 at this instant are also affected by the leading gill (gill 5) because
of the higher diffusion inherent to the lower Reω. Overall, this instant demonstrates a
combined effect of the oncoming fluid inertia (from the left), the instantaneous vortex
108
and the interaction with the leading gill. The short-lived instantaneous dipoles along
with the asymmetry introduced by the distal flap occurring over the entire array
explains the fact that most of the outflow for the flapping kinematics is occurring in
the dorsal direction (see table 4.1). In contrast, for the rowing kinematics (Figure 4.7),
the detached vortex is acting behind the retracting gill so that the instantaneous
dipole is not providing an obvious function so as to help the mean flow.
The previous discussion of vortex detachment (shedding) or attachment was
found to be close to some recent studies. The most similar work is reported by Roper
[54] in §4.2 and precisely Figure 4.4. Using COMSOL, a tethered 2D asymmetric fin,
which is a proxy for Clione antarctica [12], was forced to oscillate (flap) up and down
in a time-reversible fashion. The foil was then able to exert a net pumping force upon
the fluid in a direction orthogonal to the plane of flapping. Two flapping regimes
were identified which were signaled by a reversal in the direction of the time averaged
pumping force at Reω = 20. It was concluded that the pumping mechanism is effective
even at arbitrarily low rates of energy expenditure. At higher frequencies (Reω > 20),
the flow rate is supplied by vortex dipoles shed from the fin edges on alternating halves
of the fin stroke. However, positive thrust is generated only above a critical flapping
Reynolds number (Reω > 20) coinciding with the work in [1, 12, 66, 67]. Roper [54]
pointed out, however, that the mechanism for thrust production is not illuminated by
examination of the instantaneous flow field, since at each half of the fin beat, vortices
build up on the fin edges, but are eliminated by viscous diffusion too quickly to be
shed (Hence, my depiction earlier in 4.6 and 4.7 using normalized vorticity).
Another similar work pertaining to locomotion only (no net pumping investi-
109
gated) was presented in [1]. [1] used in-silico experiments to study a 2D flapping
ellipsoidal body. By introducing a perturbation normal to the flapping direction,
they showed that at sufficiently large Ref , unidirectional locomotion emerges as an
attracting state. Locomotion was shown to be generated in two stages: first, the fluid
field loses symmetry by an instability similar to the classical von Ka´rma´n instability;
and second, precipitous interactions with previously shed vortical structures push the
body into locomotion. Based on my current work and previous studies, it may be
concluded that net pumping occurs even at low Reω provided some sort of asymmetry
is introduced. However, locomotion occurs only above a certain critical Reω.
To further illustrate the effect of Reω on the flapping kinematics, pathlines have
been integrated after reaching the steady oscillatory state for both the highest and
lowest Reω available (Reω = 21.6 and Reω = 1.0 respectively) over 100 periods.
The result of the integration is shown in Figure 4.8. In Figure 4.8, particles were
released all around the gill array at z/L = 1.15, then integrated over 100 periods.
The initial position of the particles is outlined by black circles (◦), and the pathlines
are colored by depth z. Also, the location of the particles is shown every 10 periods
(•). In comparing Reω = 21.6 with Reω = 1.0, one may spot the increasing inertial
effects and decreasing viscous effects by noticing two features of the pathlines. The
first feature is where the particles or pathlines exit the dorsal plane (top or xmax).
At Reω = 21.6, the particles have been displaced further to the right compared to
Reω = 1.0 by about 0.5L, which is highlighted by the red arrows. The second feature
is the oscillating motion of the particles, which is present everywhere for Reω = 1.0,
but only at inflows for Reω = 21.6. Both these features indicate higher inertia present
110
(a) Reω = 21.6
(b) Reω = 1.0
Figure 4.8: Pathlines integration for flapping kinematics over 100 periods. (a) Reω =
21.6. (b) Reω = 1.0. Pathlines are colored by depth z (blue: inwards, red:outwards).
◦: Particles initial seeding position at z/L = 1.15. •: Particles location each 10
periods. Black arrows indicate mean flow direction. Section of gills shown is at
z/L = 1.15 and illustrates their extreme flapping positions.
111
in the flow. In addition, the oscillating motion of the particles with respect to the
actual displacement is an indication of how well or efficient the micro-pumping is. In
specific, it takes on average 4 beating periods for the particles located on the ventral
side (bottom or xmin) to reach the inter-gill spacings when Reω = 1.0, but it takes
fewer than 2 periods when Reω = 21.6. Qualitatively, it may be concluded that the
flow field for Reω = 21.6 provides a higher mass flow rate because of the smaller
oscillations, is doing less work in moving particles because of the fewer number of
periods required to move particles around, and is therefore more efficient than the
flow field for Reω = 1.0. This will be substantiated later when the performance
metrics are presented.
Noteworthy is that the second feature discussed above highlights the second
element of synergy previously discussed in §4.3.1 and how the increased inertial effects
are acting in favor of enhancing the antiplectic metachronal wave. Noteworthy also
is that not all inter-gill spacings are performing the same. Close examination of
Figure 4.8 will reveal that inter-gill spacings 4-5 and 5-6 are pushing more fluid than
their counterparts, followed by inter-gill spacing 3-4, then 2-3 and then 1-2. This is
equally true for both Reω = 1.0 and Reω = 21.6.
To contrast, pathlines integration for the rowing kinematics using the same
planform (planform A) for Reω = 21.6 and Reω = 1.0 is performed and is shown in
Figure 4.9. According to Figure 4.9, changing the Reω has a stronger effect on the
rowing kinematics compared to the flapping kinematics. In Figure 4.9b (Reω = 1.0),
a small recirculation zone is present on the ventral side of the mayfly, below gill pair
2 and 3. This recirculation zone does not extend all the way in the z direction but
112
(a) Reω = 21.6
(b) Reω = 1.0
Figure 4.9: Pathlines integration for rowing kinematics over many periods. (a) Reω =
21.6. (100 periods) (b) Reω = 1.0 (75 periods). Pathlines are colored by depth z (blue:
inwards, red:outwards). ◦: Particles initial seeding position at z/L = 1.15. •: Particles
location each 10 periods. Black arrows indicate mean flow direction. Section of gills
shown is at z/L = 1.15 and illustrates their extreme rowing positions.
113
extends only from the mid-section of the gills (z/L ≈ 1.15) to the their tips (z/L ≈
1.68). The recirculation zone was also observed when Reω = 2.3, but disappeared for
Reω ≥ 8.1. It is interesting that this feature is ’toggled’ during the rowing-flapping
transition, but no apparent connection could be established. The recirculation zone
is probably due to the increasing effect of the no-slip ventral wall located at ≈ 1.5L
below the gill array when decreasing Reω .
Another effect of increasing the Reω is the ratio between the medial (lateral)
flow rate to the posterior (back) flow rate. The higher the Reω, the higher the flux
through the posterior plane as shown in Figure 4.9a. On the contrary, the lower the
Reω, the higher the flux through the medial planes as shown in Figure 4.9b. This
is probably due to the increased interaction between the gills, where a retracting
gill faced by a non-moving or slowly protracting gill will tend to squeeze fluid out
medially rather than pushing fluid posteriorly. In numerical terms, the magnitude of
the y-component of the velocity vector is higher for higher Reω, which aligns with the
direction of the power stroke of the gills while the z-component is higher for lower
Reω. The mass flow rates will be quantified later in §4.6.
The rowing flow fields exhibit one feature that is similar to the flapping kine-
matics previously shown in Figure 4.8, which is the amount of oscillations present
along a pathline with respect to the net displacement of a particle. The higher the
Reω, the higher the net displacement of most particles with respect to the amplitude
of oscillations. This is the telltale sign of increasing fluid inertia. Again, just by
looking at both cases, one may speculate that the higher Reω, the more efficient is
the mayfly in moving fluid around.
114
Finally, at Reω = 21.6 and in comparison with flapping, particles are leaving
the integration domain faster. In Figure 4.9a, particles released on top of the mayfly
naiad left the integration domain in fewer than 10 periods (medially or posteriroly),
while for the flapping case at the same Reω (Figure 4.8a), most of the particles resided
in the integration domain for more than 20 periods. The is primarily due to the higher
accelerations present in retracting half-stroke of the rowing kinematics and the higher
tip velocities associated with the longer strokes.
To summarize this section, one may conclude a few points about the effect of
Reω. First and foremost and for both kinematics, the higher the Reω, the more effi-
cient the fluid is moved around. This was indicated by the smaller oscillations with
respect to the net displacement associated with higher Reω and the fewer number of
periods required to move particles around. In other words, both kinematics are en-
hanced by increasing fluid inertia and decreased viscous effects. Second, the flapping
flow structure is not affected much by increasing the Reω, while the rowing flow struc-
ture shows more sensitivity, though insignificant overall. Finally, the rowing type of
kinematics is able to move fluid around faster and is therefore inducing higher mass
flow rates.
4.5 Planform Morphology
Initially, the parametric space designed to investigate the rowing-flapping tran-
sition had only the Reω as an independent parameter, Reω = {1.0, 2.3, 8.1, 21.6},
yielding only 8 simulations, 4 for flapping and 4 for rowing. After the quantitative
115
assessment using the performance metrics introduced in §2.4, which will be shown
later, the results showed some potential to explain the transition, but did not match
the experimental findings very well. In specific, the quantitative results favored a
transition between 1.0 < Reω < 2.3 when using only planform A (see Figure 2.3 for
planform illustration), which is lower than the regime reported by the PIV experi-
ments (3 < Reω < 8). Knowing that young nymph possess narrower gills than what
was used in the simulations, I decided to investigate the effect of the gill shape or
planforms. Therefore, two more planforms were introduced, which are planforms B
and C (see Figure 2.3 for planform illustration). Planform B was reported by ex-
periments and is used by younger (smaller) nymphs. Planform C is a numerically
distorted version of planform B, obtained by scaling down planform B, but only in its
local x axis. Hence, the area of planform C is 20.72% smaller than B and its aspect
ratio is 26.01% higher than planform B.
The direction of investigating the gill planform was a result of a few key points.
First, there were only two planforms reported by the experiments, namely planforms
A and B for old nymphs and young nymphs respectively. This means that there is a
degree of uncertainty associated with the shape of those gills. In other words, there
is no guarantee that the gill planform is exactly the same across different nymphs or
even within the same nymph. Given that C. triangulifer molts more than a few times
during its aquatic phase, it seems more likely that the gill shape or planform undergo
some sort of ’smooth’ transition rather than jumping directly from planform B to
planform A. This was also speculated because the structural properties of the gills
undergo a transition (bigger gills exhibit more flexibility) and because of the existence
116
of an intermediate kinematics regime. Second, the gill length is not constant within
the same nymph (gill 3 and 4 are the biggest, then gills tend to decrease in size
posteriorly and anteriorly). Finally, during the replication of the rowing kinematics,
I investigated the average location of landmarks ’D’ and ’V’ with respect to the root
’R’ and tried to fit those on planforms A and B. This investigation was completely
independent of the Eulerian angles extraction or the kinematic chain. I was simply
calculating the norms of the different segments (RV ,RD,DV ) and overlaying them
on gill planform B. This in turn suggested that the gills needed to be narrower in
order to better fit the rowing kinematics I was using.
In light of the above, and knowing that the uncertainty in the experimental work
for younger (smaller) nymphs is higher than the uncertainty related to older (larger)
nymphs, I investigated the effect of the gill planform using the rowing kinematics.
In specific, 8 more simulations using the rowing kinematics were carried out, 4 using
planform B and 4 using planform C, both planforms at Reω = {1.0, 2.3, 8.1, 21.6},
with the purpose of comparing those with previous rowing simulations using planform
A. With the introduction of the planform uncertainty, variability or morphology, the
final parametric space encompassed 16 simulations and 3 independent variables, which
are the Reω {1.0, 2.3, 8.1, 21.6}, the gill planform {A,B,C} and the beating kinematics
{flapping, rowing}.
As discussed earlier in §2.1, the area of planforms A, B and C are 0.638, 0.449,
and 0.356 (length squared) respectively and the corresponding aspect ratio is 1.57,
2.23, and 2.81 respectively. Most of the qualitative differences that show up in com-
paring the flow field of different planforms can be attributed to the variation in their
117
(a) Planform A
(b) Planform C
Figure 4.10: Effects of planform morphology on flow fields for Reω = 21.6 at t/τ =
0.476. Iso-surfaces of the second invariant of the velocity gradient tensor (Q2=5),
flooded by the z-component of the vorticity vector ωz. (a) Planform A. (b) Planform
C.
118
(a) Planform A
(b) Planform C
Figure 4.11: Effects of planform morphology on flow fields for Reω = 21.6 at t/τ =
0.476. Iso-surfaces of the viscous dissipation rate (φ = 5 × 10−5), flooded by the
velocity magnitude. (a) Planform A. (b) Planform C.
119
surface area and their aspect ratio. After all, this is the only independent variable or
parameter that is being changed. For example, the (spatial) average velocity magni-
tude behind a retracting gill is higher, when using bigger planforms, i.e. a bigger mass
of fluid is being pushed posteriorly. This will have a direct effect on the mass outflow
exiting the gill array behind the nymphs. The higher the surface of the planform,
the higher the mass flow rates induced for the 3 planforms investigated. A second
qualitative measure is seen when iso-surfaces of flow structures are plotted. In Fig-
ure 4.10, instantaneous iso-surfaces of the second invariant of the velocity gradient
tensor (also known as the Q2 criterion [33]) are shown for the largest planform (A) and
the smallest planform (C) both operating at the highest Reω (Reω=21.6). Also, the
iso-surfaces are colored by the z-component of the vorticity vector, ωz. In addition,
both left and right rows of gills are shown to provided better visualization (remember
that the flow field is symmetric around the medial plane). Generally speaking, the
coherent structures for planform A are bigger, which is most evident by the presence
of the bound vortices around gills 3-6. Not only the vortices are larger, but also
the magnitude of vorticity is higher for bigger planforms as indicated by the colored
contours.
The instantaneous field picked in Figure 4.10 (t/τ = 0.476) actually coincides
with the instantaneous field previously shown in 4.7a when presenting the effect of
Reω on the flow field. At that instant, which is the end of retraction of gill 4, the
onset of vortex detachment was shown, which is illustrated here by the iso-surfaces
of the Q2 criterion. The bound vortex has detached from gill 4 and will quickly lose
coherence into the surrounding field. Again, the detached vortex for planform A is
120
larger that for planform C and exhibits higher vorticity magnitude especially when
moving farther from the gill tip.
With the higher velocities and bigger structures associated with bigger plan-
forms, one may speculate that the rate of viscous dissipation is also higher. For that
purpose, Figure 4.11 shows iso-surfaces of the viscous dissipation rate (φ = 5× 10−5)
flooded by the velocity magnitude at the same instant as Figure 4.10, i.e. at the end
of retraction of gill 4 (t/τ = 0.476). Clearly, the viscous dissipation rate is higher
(evident by the larger iso-surfaces) for planform A at that Reω. This is in fact true
for all Reω considered. The higher velocities near the tip of gill 4 is also evident as
well as gills 5 and 6.
In my efforts to relate the effects of planform morphology to known physical
parameters or dimensionless groups, I found that as many as three Reynolds numbers
could be defined for this particular problem. The first Re is the oscillating Re (previ-
ously defined as Reω = L2f/ν). The second Re could be based on some tip velocity
(maximum or mean), which would take into account the difference between the two
types of kinematics (including the hinge development) and would surely correlate with
Ref = fAL/ν (sometimes reported as the flapping or frequency Re). A third Re may
be defined using some chord length (for example DV ) or the gill aspect ratio, which
would take into account some of the variation in the planform shape. In an overly
simplified situation, each of the 3 independent parameters in this parametric space
(Reω, planform and kinematics) represents some ratio between the viscous forces and
the inertia in one way or another and all 3 parameters are playing a role in explain-
ing the rowing-flapping transition. Noteworthy is that in insect locomotion, where a
121
forward velocity is clearly defined, a manipulation of the first two Reynolds numbers
would yield the so-called advance ratio or reduced frequency, but for a stationary
insect like the mayfly nymph, these parameters do not show up normally.
To sum up this section, the effect of planform morphology was investigated by
simulating 3 different planforms using the rowing kinematics. This was mainly driven
by the higher uncertainty present in smaller nymphs which are known to employ the
rowing type of kinematics and use narrower planforms such as planforms B and C. It
was concluded that larger planforms are inducing higher mass flow rates and bigger
vortices according to the Q2 criteria. The presence of a bound vortices speculated by
previous PIV experiments [57–59] is now corroborated. The effect of increasing the
aspect ratio of the gills may be analogous to decreasing the inertial effects, because
of the smaller fluid mass associated with it (smaller added mass).
4.6 Micro-pumping and Performance Metrics
The results presented in the previous sections attempted to introduce the com-
plicated nature of the flow field while qualitatively comparing the effects of the dif-
ferent independent parameters involved in the rowing-flapping parametric space. In
other words, the effect of varying the Reω, the gill planform and the type of kine-
matics employed was investigated only by looking at the flow fields of the different
scenarios. In fact, the different snapshots of the flow fields were hand-picked from
dozens of animations that were batch-created and inspected in my attempt to ex-
plain the transition. Unfortunately, all of these qualitative results cannot ascertain
122
any possible reason behind the transition reported by the PIV experiments. What
was established qualitatively, with a good degree of certainty, is that higher Reω sim-
ulations are more efficient than lower Reω for both kinematics which doesn’t say much
about the transition. In fact, it is well established that lower Reω are associated with
lower efficiencies [45, 53], including roto-dynamic pumps. It was also illustrated that
the rowing kinematics are able to push fluid faster in comparison with the flapping
kinematics when operating at the same Reω. This was mainly due to the higher
accelerations and longer beating strokes present in the rowing kinematics, where the
effect of latter will appear in Ref (rowing kinematics have higher Ref than flapping
kinematics when both are operating at the same Reω). It can be speculated that the
rate of work done by the rowing kinematics is going to be higher than its flapping
counterpart when both kinematics are employing the same planform (planform A).
This, again doesn’t say much about the transition, since the rowing kinematics in-
duces higher mass flow while and, at the same time, does more work compared to its
flapping counterpart. The effect of increasing the planform aspect ratio on the rowing
kinematics is inducing less mass flow rate and probably doing less work. Narrower
planforms will decrease the interaction between the gills, which might act in favor of
the rowing kinematics since the rowing kinematics does not display much synergy. On
the other hand, the flapping kinematics were shown to have two elements of synergy,
where the first allows two consecutive gills to work together, and the second brings
the entire array to work together by enhancing the antiplectic metachronal wave.
In light of the above, it was necessary to look into quantitative measures that
try to assess the performance of each combination of independent variables. Any
123
quantitative measure ought to be realizable by the current mathematical model, which
solves only for conservation of momentum and mass. This is still acceptable owing to
the exploratory nature of the project and the higher cost associated with simulating
mass transport. Future extensions may be able to tackle the problem with more
complicated models, but only after some understanding of the flow fields have been
reached and the important parameters have been identified.
However, most of the metrics that exist for biological taxa and fluid mechanics
assume the existence of a clear forward velocity, upon which, the metric itself is
based. Obvious examples are the lift and drag coefficients associated with insects
and airplanes. On the other hand, metrics associated with transport are clearly not
an option. This leaves room only for purely innovative, and possibly approximate,
approaches. With the nature of the flow field, which is a periodically moving body
in stagnant fluid, the first engineering approach that was thought of is a CV (control
volume) approach. More importantly, since efficiency is at question, the first metric
aimed to use the conservation of mechanical energy, even though the energy equation
was not part of the mathematical model that solves for the flow around the mayfly
nymph. In the next section, the results of the first metric, η, is introduced along with
the sensitivity of the metric itself to the control volume size.
4.6.1 Mechanical Efficiency
The first performance metric is based on the integral form of the mechanical
energy equation over a fixed volume [41], namely the mechanical efficiency of a control
124
volume. The metric was previously introduced in §2.4.1 and reads:
η(%) =
∫
outer Eu · dA−
∫
outer uiτijdAj
∫
gills uiτijdAj
× 100 (4.1)
where ’outer’ refers to the surfaces defined by the six planes corresponding to the
rectangular control volume (Figure 2.11 on Page 50) and ’gills’ refer to the control
surface in contact with the moving gills. The overbar denotes time-averaged quan-
tities. The effect of the CV size on the metric was investigated by computing η for
3 different CVs. Each CV is larger than the bounding box containing the moving
gill array but does not include them. The 3 CVs are simply rectangular prisms and
are approximately 0.5L, 1.0L and 2.0L away from the bounding box containing the
moving gill array. The 3 CVs are defined by the following coordinates respectively:
• CV1: extends from (-0.59,-1.70,-2.21) to (1.41,1.70,2.21).
• CV2: extends from (-1.09,-2.70,-2.71) to (1.91,2.70,2.71).
• CV3: extends from (-1.40,-3.70,-3.73) to (2.95,3.70,3.71).
where CV1 is the smallest and CV3 and is the largest.
The mechanical efficiency for the rowing-flapping parametric space is shown in
Figure 4.12 for all 3 CVs as a function of Reω. For all 3 planform types and across
the regime considered, the mechanical efficiency of the rowing kinematics is always
higher than that of the flapping kinematics. Also, within the 12 rowing simulations,
larger planforms result in a higher mechanical efficiency. The results of the mechanical
efficiency is not analogous to the kinematic transition in insect locomotion, where a
similar propulsive efficiency favors flapping above some critical Reynolds number [12].
125
Therefore, it is concluded that the transition from rowing to flapping kinematics,
that has been observed in-vivo, is probably not done on the basis of optimizing the
mechanical efficiency.
While the trends of η did not provide a viable explanation to the kinemat-
ics transition, they still provide some insight on the fluid behavior in the regime
investigated. Therefore, the different terms involved in computing the mechanical
efficiency where investigated independently in an attempt to explain the trends ob-
served, especially the increase of η for small CVs at the low Reω side. Recalling that
all the terms are dimensionless and that the findings are for the investigated regime
(1.0 ≤ Reω ≤ 21.6), the key features of the different terms are highlighted hereunder:
For the same CV, an increase in Reω from 1.0 to 21.6 will resulted in:
• a decrease in the average rate of viscous work done by the moving gills by
approximately one order of magnitude.
• a decrease in the average rate of pressure work done by the moving gills by
approximately one order of magnitude.
• a decrease in the average rate of viscous dissipation by approximately one order
of magnitude.
• an increase in the average kinetic energy flux through the control surface by
one to two orders of magnitude depending on the control volume size. In other
words, the fluid gains more kinetic energy as it leaves the control volume.
• an increase in the average rate of work done by the control volume on the
126
Figure 4.12: Mechanical efficiency, η, for 3 control volumes. ∗: Rowing using planform
A. : Rowing using planform B. ◦: Rowing using planform C. +: Flapping using
planform A. Top: CV1 (0.5L away). Middle: CV2 (1.0L away).Bottom: CV3 (2L
away).
127
surrounding fluid by one to two orders of magnitude depending on the control
volume size.
Also, for the same Reω, an increase in the CV size from CV1 to CV3 resulted in:
• a decrease in the average kinetic energy flux through the control surface by one
to two orders of magnitude depending on the Reω. The lower Reω, the faster
the decrease.
• a decrease in the average rate of work done by the control volume on the sur-
rounding fluid by one to two orders of magnitude depending on the Reω. The
lower Reω, the faster the decrease. Only rowing using planform A exhibits a
decrease by three orders of magnitude at the lowest Reω.
• a negligible change in the average rate of viscous dissipation. In other words,
most of the energy loss is occurring in the near field around the moving gills.
The above provides a numerical explanation to the observed trends of η with
respect to Reω. On the high end of Reω, both the average kinetic energy flux and the
average rate of work done by the control volume on the surrounding fluid are of the
same order of magnitude, but the former is increasing its order of magnitude while the
latter is decreasing it. The end result is an increase in the mechanical efficiency η on
the high end of Reω. In other words, for Reω ≥ 21.6, η is driven mainly by the average
kinetic flux, while for lower Reω, η is driven by the average stress work on the outer six
planes. Otherwise put, the first term in the numerator of equation 4.1 is dominating
at high Reω, while the second term is dominating for lower Reω. For bigger CVs,
128
both terms in the numerator are approaching zero and the competition between these
two terms is washed out by their relatively dwindling values. In fact, bigger CVs will
require higher order interpolation and integration schemes in order to compute η.
The trends observed when increasing Reω align with inertial flows, as for example
in turbo-machinery, where the second term in the numerator is neglected along with
the viscous effects resulting in the so-called hydraulic power or the pressure work rate
on the boundaries. The increasing mechanical efficiency for Reω ≥ 8.1 also agrees
with known trends in inertial flows, where the viscous dissipation keeps decreasing
for increasing Reω.
The analysis of the different terms also showed that for the same Reω and the
same CV, all terms that appear in equation (4.1) (kinetic energy flux, stress work rate
and gill work rate) are the highest for rowing using planform A, followed by rowing
using planform B, then rowing using planform C and finally flapping using planform
A. With the order of magnitude discussion above, this would explain why the rowing
kinematics have higher mechanical efficiency than the flapping kinematics.
Noteworthy is that the average total stress work (pressure and viscous) on the
outer six planes was always negative for the entire parametric space. In other words,
the CV was, on average, doing work on the surrounding fluid. However, the instanta-
neous values of the work rate behave differently. To illustrate that, the viscous work
rate and pressure work rate on a plane below the mayfly is shown in Figure 4.13. Fig-
ure 4.13 shows the viscous and pressure work rate of both kinematics using planform
A at the highest and lowest Reω at a plane located at x/L = −0.59. This is one of
the six planes that define the smallest CV of the 3 CVs discussed earlier (CV1). With
129
reference to Figure 4.13, it can be seen that, for that plane and for Reω = 1.0, the
magnitude of the viscous work rate is comparable to the pressure work rate, while for
Reω = 21.6, the magnitude of the pressure work rate is relatively higher. In addition
the fluctuation and the average viscous work rate is decreasing with increasing Reω,
which is just another indication of decreased deformation work or viscous dissipation.
Furthermore, it can be seen that both the pressure and viscous work rate are higher
for the rowing kinematics for both Reω.
Interestingly, the signature of the metachronal wave or the beating gills could
be detected for some of the planes of integration defining the CV. For example, with
reference to Figure 4.13b, the observed 4 local maxima of the pressure and viscous
work rate at t/τ ≈ {0.16, 0.45, 0.7, 0.92} correspond approximately to the end of the
retraction half-stroke of gills {4, 3, 2&6, 1&5} respectively (see Figure 2.8, on Page 44
for flapping kinematics). Likewise, in Figure 4.13d, the observed 4 local maxima of
the viscous work rate at t/τ ≈ {0.25, 0.25, 0.45, 0.75} correspond approximately to the
the end of the retraction half-stroke of gills {6, 5, 4, 3} respectively (see Figure 2.10,
on Page 46 for rowing kinematics). The local maxima of the stress work become
more difficult to spot for higher Reω or bigger CVs due to inertial and viscous effects
respectively.
The fact that the mechanical efficiency failed to explain the rowing-flapping
transition is still interesting on a fundamental level. While it may be completely
intuitive that the mayfly nymphs are not optimizing for mechanical efficiency since it
is already established that mayfly nymphs beat their gills in order to enhance their
metabolic respiration, the question that remains is that why efficient mixing was not
130
(a) Flapping at Reω = 21.6 (b) Flapping at Reω = 1.0
(c) Rowing at Reω = 21.6 (d) Rowing Reω = 1.0
Figure 4.13: Viscous and pressure work rate on the ventral plane (x/L = −0.59) for
flapping and rowing kinematics using planform A.
associated with an optimized mechanical efficiency? A possible answer was discussed
by Purcell [53], where he pointed out that in tiny propelling organisms, the energy
invested in moving around represents a small fraction of their total energy budget. In
other words, propelling insects may not care about mechanical efficiency as much as
engineers do. The most likely outcome of the analysis of the mechanical efficiency of
the mayfly is that enhanced mixing does not correlate with an optimized mechanical
efficiency in the investigated regime. An outcome that is only verifiable using more
131
specific tools such as the Lagrangian Coherent Structures approach [29]. The results
of the mechanical efficiency concluded that the rowing-flapping transition does not
occur due on energetic reasons only, and that another factor must be involved. Con-
sequently, a second performance metric was devised, namely the specific mass flow
rate.
4.6.2 Specific Mass Flow Rate
With the failure of the mechanical efficiency in explaining the rowing-flapping
transition, a metric that combines both energy and mass was thought of. The specific
mass flow rate ξ was presented in §2.4.2 and reads:
ξ =
m˙
W˙ gills
, (4.2)
where m˙ is the mass flow rate leaving (or entering) the control volume, and W˙gills is the
rate of work done by the six pairs of moving gills. ξ is interpreted as the dimensionless
time-averaged mass flow rate leaving the mayfly divided by the dimensionless time-
averaged rate of work done by the gills.
The metric was inspired from a known performance parameter in internal com-
bustion engines, namely the specific fuel consumption. In internal combustion en-
gines, the specific fuel consumption is the flow rate of fuel consumed by the engine
(say in gallons/minute) divided by the output power obtained at the vehicle’s tires
(say in Horsepower), which translates to the vehicle speed depending on the drag
force. The lower the specific fuel consumption, the more efficient is the moving ve-
hicle. The opposite holds of course. It is well known that this metric is lower on a
132
’highway’ and higher ’in city’ and hence a vehicle running at highway speeds consumes
less fuel per unit distance and is, consequently, more desirable.
The metric introduced here carries the same idea, except that a higher specific
mass flow rate is desirable. In other words, the higher the flow rate induced by the
mayfly nymph per unit power expended by the moving gills, the more ’efficient’ is the
mayfly nymphs in circulating water. Again, the metric is limited by the mathematical
model and carries the assumption that a higher flow rate across the control volume
and hence across the mayfly translates into better inter-gill flow, enhanced ventilation
and gas exchange. The main advantage of this metric compared to the mechanical
efficiency is that it incorporates both mass flow and energy, where the former may
be related to gas exchange, while the latter is be related to metabolism and energy
expenditure.
As with the mechanical efficiency, η, the specific mass flow rate, ξ was computed
for 3 different CVs in order to investigate its sensitivity to the CV size. The same
CVs used to compute the mechanical efficiency were used to compute ξ. These are
CV1, CV2 and CV3 which are approximately 0.5L, 1L and 2L away from the moving
gills. Note that for both metrics, only the numerator is function of the CV size, since
all selected CVs are bigger than the bounding box containing the moving gill array.
Therefore, only the numerator is sensitive to the CV size. First, the results for CV2
(1L away) are introduced. This is followed by a discussion of the sensitivity of ξ to
the CV size. Last, a more detailed quantification of the mass flow rates that occur
across the control surface is presented.
Figure 4.14 illustrates the metric ξ computed for CV2, where the abscissa is
133
the average rate of work done by the moving gills (denominator in (4.1)) and the
ordinate is the mass flow rate induced across the control volume (numerator in (4.1)).
Therefore, the slope at any data point in the figure represents ξ. Furthermore, the
figure shows the metric for the entire parametric space (i.e. 16 data points) and hence
introduces the 2 dependent variables (average mass flow rate and average rate of work
done) as function of the 3 independent variables being studied (Reω, kinematics and
planform). Otherwise put, the plot introduces a total of 5 variables.
Figure 4.14: Specific mass flow rate , ξ, for the rowing-flapping parametric space.
∗: Rowing using planform A. : Rowing using planform B. ◦: Rowing using planform
C. +: Flapping using planform A.
134
First, we consider the four data points corresponding to the flapping kinemat-
ics with planform A (+), which will serve as a baseline case for later comparison
with the rowing kinematics. In this case, the average rate of work done increases
(W˙ gills = {4.66, 9.11, 27.30, 60.68}) with decreasing Reω (Reω = {21.6, 8.1, 2.3, 1.0}).
In fact, one order of magnitude drop in Reω causes an order of magnitude in-
crease in the rate of work done. Both pressure and viscous work are contributing
to the rate of work done with the ratio between the former and the latter being
{5.80, 5.47, 5.38, 5.36} respectively. On the other hand, the average mass flow rate
through the control volume is decreasing (m˙ = {1.32, 0.86, 0.60, 0.52}) probably due
to the increase in the average rate of viscous dissipation. Note the four straight lines
in the figure, each starting at the origin and going through one of the four points
of the flapping kinematics. The slope of each line is equal to ξ, which is decreasing
with decreasing Reω (ξ = 1e − 03{281.88, 94.04, 22.15, 8.63}). In particular, ξ drops
two orders of magnitude with one order of magnitude drop in Reω. Overall for the
flapping kinematics, the performance based on ξ is highest at the maximum Reω
considered, which nicely coincides with the qualitative observation of the pathlines
previously presented in §4.4 and with the results of the mechanical efficiency η in
§4.6.1.
Second, we consider the twelve remaining cases with rowing kinematics (three
different planforms at four different Reω). The same symbols are used for each plan-
form (∗,  and ◦ for planforms A, B and C respectively). Note that the symbols
for the same Reω have the same color and are clustered in the neighborhood of the
corresponding lines introduced above, whose slope is the specific mass flow rate, ξ, for
135
the flapping kinematics. Overall, the same trends observed in the flapping kinematics
can also be seen for the rowing kinematics. An increase in Reω, for example, causes
an increase in the average mass flow rate and a decrease in the average rate of work
done. Among the things to note, is that for any planform and while fixing Reω, the
average rate of work done of the rowing case is higher than the flapping case and so
is the average mass flow rate. Also, as previously discussed in §4.6.1, the bigger the
planform area, the higher the rate of work done.
Now, we compare the specific mass flow rate of the rowing cases to the corre-
sponding flapping kinematics. Let us consider clusters of data points for the same
Reω, and use the line going through the flapping kinematics point as a reference.
When the data point for the rowing kinematics is: i) below the line for the particular
Reω then ξ for this case is smaller and it is less efficient than flapping; ii) above the
line for the particular Reω then ξ for this case is larger and it is more efficient that
flapping. Starting at the lowest Reω (Reω = 1.0, black in the figure), we see that all
data points for the rowing kinematics are above the line indicating that ξ is higher
and therefore rowing is more efficient than flapping for all three planforms considered.
As Reω increases to Reω = 2.3, the data point for planform A moves below the cor-
responding line, indicating that rowing using planform A at this Reynolds number is
now less efficient than flapping. Planforms B and C (which are closer in morphology
to the ones found in younger nymphs) are still more efficient than flapping. At the
two highest Reynolds numbers, all data points are on or below their corresponding
lines. For Reω = 21.6, and for all planforms, flapping is always more efficient than
rowing, while for Reω = 8.1 flapping is equally or more efficient than rowing depend-
136
ing on the planform. In other words, if a mayfly were to have gills in the shape of
planform A throughout its entire aquatic phase, the specific mass flow rate metric
would favor a switch from rowing to flapping within 1.0 < Reω < 2.3 as it grows.
Of course, the gill undergoes a change in shape from planform B or C to that of A,
in which case ξ suggests a demarcation from rowing to flapping somewhere within
2.3 < Reω < 8.1. Simulations at Reω = 2.3 clearly show the effect of the planform on
the kinematics; planform A favors flapping, while planforms B and C favor rowing.
Strictly speaking, and based on the numerical results alone, the transition regime
from rowing to flapping kinematics is 2.3 < Reω < 8.1.
Another way of presenting the previous results is shown in Table 4.2, where the
average rate of work done by the gills, the average mass flow rate and the specific mass
flow rate, ξ, are listed for all rowing cases, normalized by their flapping counterparts
at the same Reω. Although it was shown that with increasing Reω, the average mass
flow rate increases and the average rate of work done for either kinematics decreases,
the normalized rate of work done: i) is increasing with increasing Reynolds number;
ii) and the normalized mass flow rate is decreasing, indicating that rowing is falling
behind with increasing Reω no matter which planform is used thereby indicating an
imminent switch to flapping.
To further illustrate the effect of planform morphology on the different depen-
dent parameters like the rate of work done W˙ , the induced mass flow rate m˙ and the
specific mass flow rate ξ, Figure 4.14 was re-plotted while normalizing both axis by
the total surface area of the moving gills (24 × the planform area). In other words,
the x-axis is now the dimensionless time-averaged rate of work done per unit area of
137
Reω
W˙ row,A
W˙ flp,A
W˙ row,B
W˙ flp,A
W˙ row,C
W˙ flp,A
m˙row,A
m˙flp,A
m˙row,B
m˙flp,A
m˙row,C
m˙flp,A
ξrow,A
ξflp,A
ξrow,B
ξflp,A
ξrow,C
ξflp,A
1.0 2.45 1.72 1.30 2.54 2.12 1.73 1.03 1.23 1.33
2.3 2.55 1.81 1.37 2.43 2.06 1.67 0.95 1.13 1.22
8.1 2.91 2.11 1.58 2.27 1.95 1.54 0.78 0.92 0.97
21.6 3.49 2.55 1.88 2.05 1.85 1.51 0.58 0.71 0.79
Table 4.2: Table showing the average rate of work done W˙ gills, average mass flow rate
m˙, and specific mass flow rate ξ for the rowing cases normalized by their flapping
counterparts. row:rowing. flp: flapping. A,B,C: planforms A,B,C. Boldface notation
(ξ ≈ 1) highlights the transition.
the 6 pairs of moving gills while the y-axis is the dimensionless time-averaged induced
flow rate per unit area of the 6 pairs of moving gills. The result of this normalization
is shown in Figure 4.15. Inspecting Figure 4.15 doesn’t change any of the conclusions
deduced from its non-normalized version in terms of the rowing-flapping transition,
i.e. the Figure still suggests a transition in the range 2.3 < Reω < 8.1 if the effects of
planform morphology are taken into account. Figure 4.15 illustrates clearly the effect
of planform morphology on rowing simulations. On a per unit surface area basis, the
rate or work done does not change much by changing the planform area especially
for higher Reω cases. On the contrary the mass flow rate (again per unit gills area) is
decreasing with increasing planform A. This is, in fact, the exact opposite behavior
of the non-normalized version of this figure, i.e the increase in planform area results
in a decrease in the mass flow rate per unit area. Evidently, the increased inertial
effect associated with bigger planforms is not favored by a rowing type of kinematics.
This may also explain why planform A favors a switch that is earlier than the one
observed in-vivo. The normalized results are particularly interesting since the surface
138
area of the gills is definitely tied to the amount of gas exchange that occurs across its
surface.
Figure 4.15: Specific mass flow rate , ξ, for the rowing-flapping parametric space.
Rate of work done and mass flow rate are normalized using the surface are of the
moving gills . ∗: Rowing using planform A. : Rowing using planform B. ◦: Rowing
using planform C. +: Flapping using planform A.
Overall, the results indicate a strong correlation between Reω and the proposed
metric, ξ, which in turn provides a surprisingly simple, yet elegant explanation of the
kinematic switch observed in-vivo. However, before finalizing this result, the effect of
the control volume size on ξ had to be investigated. Therefore, ξ was analyzed for CV1
(0.5L away) and CV3 (2L away) and plotted the same way as for CV2. The plots are
139
shown in Figure 4.16. Analysis of Figure 4.16 yields almost the same result obtained
when analyzing Figure 4.14. This is because ξ stills favors a switch from rowing using
planform B or C to flapping using planform A in the same regime 2.3 < Reω < 8.1
suggested by CV2. In addition, when using planform A for both rowing and flapping,
ξ based on smallest CV (CV1) still favors a switch within 1.0 < Reω < 2.3 which is
the same regime suggested by CV2. The only data point that seems to weaken the
sensitivity analysis belongs to the rowing kinematics using planform A at Reω = 1.0
for the biggest CV (CV3), where this point actually suggests that the transition from
rowing using planform A to flapping using planform A occurs at Reω < 1 as opposed
to the smaller CVs which suggested a transition at 1.0 < Reω < 2.3.
The existence of this point, which represents 1 point out of 48 points, does not
conflict with the in-vivo observations. This is simply because mayflies nymphs don’t
use planform A to row. In fact, they use a planform closer in shape to planform
B or C. Hence, whether the rowing-flapping transition using planform A occurs at
1.0 < Reω < 2.3 (suggested by the two smaller CVs) or at Reω < 1.0 (suggested by the
biggest CV) may not be important to mayflies in specific. However, on a fundamental
level, it would be much more desirable to show that the kinematics switch occurs even
if the planform shape remains the same, i.e. the dependance is only on the Reω. At
this point, there were two choices: 1) run simulations at Reω < 1.0 which is more
expensive and wouldn’t guarantee that bigger control volumes (bigger than CV3)
won’t suggest more conflicting results or 2) keep working on the sensitivity analysis
and look into why ξ behave the way it does, which would bring more understanding to
the flow field and ξ. The second choice was favored because it provides more insight
140
Figure 4.16: Sensitivity of specific mass flow rate ξ to the control volume size. ∗: Row-
ing using planform A. : Rowing using planform B. ◦: Rowing using planform C.
+: Flapping using planform A. Top: CV1 (0.5L away). Bottom: CV3 (2L away).
141
and actually studies the proposed metric.
Recalling that only the mass flow rate through the control volume affects the
value of ξ, it is clear that one needs only to address how the mass flow varies with
the choice of CV, the kinematics and the Reω. It is clear that the mass flow rate
will start to drop as one gets further away from the mayfly and that as the CV size
approaches the entire domain size, the mass flow rate, and therefore ξ, is going to
be effectively zero owing to the no penetration boundary conditions. To explore that
effect, the contribution from each plane of the CV was analyzed for all CVs and all 16
simulations and graphed in manner aimed to quantitatively highlight this difference.
The plots are shown in Figure 4.17 and Figure 4.18.
142
F
ig
ur
e
4.
17
:
Il
lu
st
ra
ti
on
of
th
e
m
as
s
flo
w
ra
te
th
ro
ug
h
th
e
co
nt
ro
l
su
rf
ac
e
fo
r
th
e
ro
w
in
g-
fla
pp
in
g
pa
ra
m
et
ri
c
sp
ac
e.
57
6
da
ta
p
oi
nt
s
w
er
e
us
ed
to
bu
ild
th
is
fig
ur
e
(1
6
si
m
ul
at
io
ns
×
3
co
nt
ro
l
vo
lu
m
es
×
6
pl
an
es
fo
r
ea
ch
co
nt
ro
l
vo
lu
m
e
×
in
flo
w
an
d
ou
tfl
ow
fo
r
ea
ch
pl
an
e
=
57
6
p
oi
nt
s)
.
N
ot
e
th
e
sl
ig
ht
ly
de
cl
in
in
g
m
as
s
ou
tfl
ow
(o
r
in
flo
w
)
fo
r
ro
w
in
g
us
in
g
pl
an
fo
rm
A
fo
r
C
V
3
at
R
e ω
≤
8.
1.
143
F
ig
ur
e
4.
18
:
Il
lu
st
ra
ti
on
of
th
e
m
as
s
flo
w
ra
te
th
ro
ug
h
th
e
co
nt
ro
l
su
rf
ac
e
fo
r
th
e
ro
w
in
g-
fla
pp
in
g
pa
ra
m
et
ri
c
sp
ac
e.
P
er
ce
nt
ag
e
is
of
to
ta
l
in
flo
w
or
ou
tfl
ow
.
57
6
da
ta
p
oi
nt
s
w
er
e
us
ed
to
bu
ild
th
is
fig
ur
e
(1
6
si
m
ul
at
io
ns
×
3
co
nt
ro
l
vo
lu
m
es
×
6
pl
an
es
fo
r
ea
ch
co
nt
ro
l
vo
lu
m
e
×
in
flo
w
an
d
ou
tfl
ow
fo
r
ea
ch
pl
an
e
=
57
6
p
oi
nt
s)
.
144
Figure 4.17 summarizes the effect of the Reω, the planform and the beating
kinematics on the total inflow and outflow to the different control volumes, a quan-
tification of earlier qualitative comparisons of flow fields. In addition, it illustrates
the contribution from each plane of the CV to the total inflow and outflow separately.
In the horizontal direction, the figure contains 4 sets, one set for each Reω. Each set
is divided into 4 subsets corresponding to flapping using planform A, rowing using
planform A, rowing using planform B and rowing using planform C. Finally, each
subset is divided into 3 sub-subsets corresponding to CV1, CV2 and CV3. Vertically,
each sub-subset is a bar plot colored by the plane location and illustrates the contri-
bution of each plane to the inflow (negative) and outflow (positive). In total, the plot
is based on 576 data points. Needless to mention that mass is conserved through the
CV, which is indicated by the top-bottom symmetry of the figure. Also, the left-right
symmetry of the flow is indicated by the equal length of the left and right medial
planes (orange and red).
By looking only on the total inflow and outflow of all simulations, most of the
general trends discussed qualitatively earlier can be deduced from Figure 4.17. For
example, for both kinematics and all planforms, increasing Reω increases the mass flow
rate. Also, for rowing kinematics, bigger planforms result in higher mass flow rates.
Furthermore, it can be seen that flapping kinematics produces lower mass flow rates
compared to rowing and that for lower Reω the difference becomes more pronounced,
i.e. rowing is getting better in comparison with flapping as Reω is decreased.
Now, looking on the contribution of each plane within each sub-subset can
provide us with more details about the changes occurring in the flow field. Most
145
interestingly is the contribution of the medial planes to the total outflow in the rowing
kinematics (orange and red). The mass outflow through this planes does not change
much with decreasing Reω, which means that the medial component of the flow is
increasing with respect to the total outflow as Reω decreases. This was qualitatively
illustrated earlier using pathlines integration for rowing at the highest and lowest
Reω (Figure 4.9). One can also tell the major contributions from the different planes
that highlight the main difference in kinematics. For rowing, the outflow is mainly
occurring anteriorly (cyan) and medially (orange and red), while the inflow is mainly
dorsally (blue) and medially (orange and red). In contrast, for flapping, the major
contributors for outflow are the dorsal plane (blue) and the anterior plane (cyan),
while the inflow is distributed among the medial planes (orange and red) and the
posterior plane (yellow).
By far, the most important observation from this figure is how the mass flow
changes in going from CV2 to CV3 when rowing using planform A. For Reω ≤ 8.1, the
mass flow through the largest control volume (CV3) is smaller than that through the
medium-sized control volume (CV2). The same can’t be said about flapping using
planform A, i.e. CV3 is yielding higher flow than CV2. In short, for Reω ≤ 8.1
and for CVs equal or greater to CV3, the mass flow through the control surface
for rowing using planform A is already declining as one gets further away from the
mayfly, while that for flapping is not. This is caused primarily by the smaller outflow
from the ventral integration plane (dark blue), which is getting closer to the domain
boundaries (ventral integration plane for CV3 is at x = −1.4, while domain boundary
is at x = −1.5).
146
Figure 4.18 is very similar to Figure 4.17, but instead of plotting the absolute
values of the mass flow rate through the different control surfaces, the percent contri-
bution of the control surfaces of the outflow or inflow is shown. In addition, the bar
plots have been reshuffled so that the highest 4 sets denote flapping using planform A,
rowing using planform A, rowing using planform B and rowing using planform C. In
addition, the subsets are arranged by the Reω and then the sub-subsets are arranged
by the control volume. Again, the different trends can be observed by inspecting this
figure such as the effect of decreasing Reω on the rowing type of kinematics using
any of the 3 planforms, which is an decreased medial inflow and an increased medial
outflow.
With the quantitative analysis of the different mass fluxes and the sensitivity
analysis of ξ to 3 different CVs, it is concluded that the observed trends of ξ are a
good representative of the 3 independent parameters and that its sensitivity to the
CV size is not a major factor.
To sum up the results of this chapter, first the main differences in the flow
fields were introduced qualitatively. This included a detailed discussion of the 3
independent parameters being studied, namely the type of kinematics used in §4.3,
which was followed by a discussion of the viscous and inertial effects in §4.4 and finally
the effects of planform morphology in §4.5. Time-averaged flow fields were discussed
and illustrated clearly the primary difference between both types of kinematics in
that the resulting outflows are perpendicular with respect to each other. Pathlines
integration showed that both types of kinematics are enhanced by increasing inertial
effects. Pathlines integration also pinpointed one of the effects of decreasing the Reω
147
on the rowing kinematics, which is an increased medial outflow, which is probably
due to increased gill interaction. Vorticity fields showed that vortex detachment
occurs during the reported transitional regime for both kinematics and is unlikely to
be coincidental. Iso-surfaces of the second invariant of the velocity gradient tensor
showed that bigger planforms were associated with bigger coherent structure and
higher vorticity magnitudes. Instantaneous fields showed that flapping offers a more
synergistic micro-pumping mechanism by bringing the entire array to work together.
Instantaneous fields also showed that rowing gills tend to work in isolation.
Second, quantitative metrics have been devised and computed for the entire
rowing-flapping parametric space. The metrics were restricted by the mathematical
model used, which basically solves for mass and momentum conservation without
incorporating mass transport due to the cost-prohibition of the latter at the time of
the study. The first metric, namely the mechanical efficiency η, failed to present any
explanation behind the transition. This was in contrast to propelling animals where
the mechanical efficiency of flapping animals is higher. The first metric, however,
indicated that energetic reasons may not be the only reason behind the transition and
that transport may not necessarily correlate with the mechanical efficiency within the
investigated regime. Therefore, other factors must be involved and hence the specific
mass flux ξ was introduced.
As simple as it may be, the specific mass flux showed a strong correlation with
the mayfly nymph kinematics and hence favored a kinematic switch within 2.3 <
Reω < 8.1, which matches surprisingly well with the observed in-vivo experiments
(Reω ≈ 5). The regime will be 1.0 < Reω < 2.3 if planform morphology is not taken
148
into account. The sensitivity of the second metric to the control volume size was
meticulously inspected and it was concluded that the control volumes closest to the
gills do provide an adequate interpretation of the results. In addition, the sensitivity
analysis shed more light on the effects of the different independent variables on the
flow field.
Future studies should look into solving for transport while avoiding the higher
computational cost associated with it. This may be achieved by using specialized
tools (e.g. Lagrange Coherent Structures or Monte Carlo simulations), which would
themselves require rigorous validation. Additionally, a study that incorporates trans-
port while studying mayfly nymphs is confronted by another challenge, which is the
unknown mass transport boundary condition on the gills surface and the correct for-
mulation to treat this boundary condition within an immersed boundary formulation.
Only with these types of studies that the actual relation between ξ and mass trans-
port can be revealed. I expect that some relation exists when the mass transport
is convection-dominated but vanishes for diffusion-dominated flows. The P e´ for the
diffusion of dissolved oxygen in water at 25◦C varies from ≈ 103 for small nymphs
to more than 104 for larger nymphs, i.e. more than one order of magnitude change,
which might profoundly alter the mass transport physics throughout ontogeny. Un-
doubtedly, there is room for future improvement.
149
Chapter 5: Results II: Flexibility and the Role of the Hinge
5.1 Introduction
In the first stages of the project, the kinematic chain introduced to simulate the
gills motion contained only two levels, where the second level describes the motion
of the entire gill with respect to the abdomen. However, comparison with in-vivo
experiments wasn’t that successful. Therefore, the third level was introduced, which
describes the orientation of the distal plate with respect to the proximal plate by in-
troducing the hinge flexion angle γ(t), while the second level describes the orientation
of the proximal plate with the respect to the abdomen (as opposed to the entire gill
in the first stage). The rowing-flapping parametric space was then executed using
the 3-level kinematic chain. Recall that the flexural hinge is coupled only will the
flapping kinematics, while the rowing kinematics has negligible hinge flexion.
The investigation of the rowing-flapping transition showed a favorable effect of
the flexural hinge, in that the location of the instantaneous dipoles formed by the
end of a retracting gill enhances the second element of synergy discussed in §4.3.1. In
specific, the flexural hinge enables the instantaneous dipole to eject water in a more
horizontal direction, which aligns the flow partially to the next posterior gill, thus
increasing synergy between neighboring gills. However, I couldn’t ascertain if such
150
short-lived instantaneous dipole is a result of the gill motion only or is a consequence
of the oncoming flow possessing increasingly ’weak’ inertia or a combination of both.
Is it possible that the oncoming flow decreases the vorticity decay near the dipole,
thereby enabling the dipole to live longer? How important is the interaction of the
neighboring gills? What is exactly the role of the hinge? [57] argued that the distal
flap acts as a non-return valve that decreases the retrograde dorsal flow into the inter-
gill spacing. Is there another role associated with the hinge development? With that
in mind, the second parametric space was born. Noteworthy is that the introduction
of the hinge is a step towards a more flexible gill, since it adds one degree of freedom
to the entire gill.
The second parametric space is a set of 8 simulations, 4 running at Reω = 21.6
and 4 running at Reω = 2.3. Each subset of 4 simulations includes 2 simulations
using the full array of 7 gill pairs (6 pairs moving and 1 pair fixed) and 2 simulations
using only gill pair 4 (i.e. gills 1-3 and 5-7 are removed). Finally, each sub-subset of
2 simulations includes 1 simulation that uses the one-plate model (without flexural
hinge) and 1 simulation that uses the two-plate model (flexural hinge). The mayfly
body (abdomen, thorax and head) was included in all simulations. Planform A was
used for all 8 simulations and the proximal plates are prescribed with the same flap-
ping kinematics previously discussed. As with the rowing-flapping parametric space,
three independent parameters are at play. These are the Reω ({2.3, 21.6}), the num-
ber of gills ({1, 7}) and the flexural hinge ({with hinge, without hinge}). In short,
each of the three independent parameters is varied once making the total number of
simulations 8 (23), where 2 of these overlaps the rowing-flapping parametric space.
151
The 8 simulations were adequate to elucidate the role of the hinge. Simulations us-
ing different numbers of gills will be required if one were to determine the optimum
number of gills in an array, but this was not the primary goal of this investigation.
As with the rowing-flapping parametric space, first, the different flow fields are
presented and compared qualitatively. This will be followed by a discussion of the
specific mass flow rate and the different forces acting on the gills and finally concluding
remarks.
5.2 Qualitative Investigation
To elucidate the effect of the spatial asymmetry introduced by the flexural
hinge within the full array of gills, time-average flow fields for both cases are shown
at Reω = 21.6 in Figure 5.1 and Figure 5.2. In both figures the two-plate model (with
flexural hinge) is compared with the one-plate model (without flexural hinge). Also,
in both figures, time-averaged streamlines, gills location and their intersection with
the slice (z/L = 1.15) are shown. Additionally, the dipole centers (vortex cores) are
connected with solid lines in order to illustrate their average location. In both figures,
it is evident that the introduction of the flexural hinge increases the antero-posterior
component of the velocity (v), which is illustrated by the colored arrows on the dorsal
side.
Figure 5.1 shows the time-averaged z-component of the vorticity ωz, while Fig-
ure 5.2 shows the time-averaged normalized velocity magnitude. With reference to
both figures, many features may be distinguished. The first feature is the average
152
Figure 5.1: Time-averaged vorticity component ωz for Reω = 21.6 using the full array
of gills. Slice at z/L = 1.15. Gill edges are shown in brown. Gills intersection
with slice is shown in black. Average dipole centers are connected with black lines.
Streamlines going through each inter-gill spacing are shown. Top: with flexural hinge.
Bottom: without flexural hinge.
153
Figure 5.2: Time-averaged normalized velocity magnitude for Reω = 21.6 using the
full array of gills. Velocity magnitude is normalized by the local maximum. Slice at
z/L = 1.15. Gill edges are shown in brown. Gills intersection with slice is shown
in black. Average dipole centers are connected with black lines. Streamlines going
through each inter-gill spacing are shown. Top: with flexural hinge. Bottom: without
flexural hinge.
154
magnitude of both cores of the dipole. The negative vorticity, which is generated
during the protraction of the gill, is, on average, higher for the two-plate model in
comparison with the one-plate model (global minimum is −6.71 vs. −4.94) and is
stretched along the distal plate. This is most evident for gills 3-6. This is primarily
caused by the protracting distal plate being tangential to the flow as opposed to the
one-plate model. On the other hand, the positive vorticity, which is generated during
the retraction of the gill is, on average, slightly lower for the two-plate model (global
maximum is 4.13 vs. 4.67). The is primarily due to the time-averaging process since
the instantaneous fields show that the vorticity magnitude is higher for the two-plate
model at all times.
The second feature is the orientation of the dipole, where in the two-plate model,
the orientation is more vertical, when compared to the one-plate model. The more
vertical orientation was highlighted in the first parametric space in §4, which enhances
the second element of synergy by aligning the flow from a preceding gill to the next
one. This has the effect of increasing the antero-posterior component of the velocity
(v), which is illustrated by the colored arrows on the dorsal side.
The third feature pertains to the normalized velocity magnitude, where in the
one-plate model, local maxima seem to be isolated near the dorsal edge of each gill,
while in the two plate model, the local maxima near gills 3-5 seem to blend together.
This is primarily due to the longer tip trajectory associated with the two-plate model,
which allows the distal plates to be closer to each other. One also notes that the
local maxima of the normalized velocity magnitude occurs at the center of the time-
averaged dipoles, which indicates that the dorsal dipoles are much stronger than the
155
ventral dipoles for both models.
Based on the above discussion, it is still inconclusive which model is performing
better. We know that the two-plate model shows better agreement with the in-vivo
experiments. We also know that the two-plate model increases the y-component of the
velocity by breaking the spatial symmetry (distal flap). Moreover, we know that both
models break the ’orientational’ symmetry with respect to the ventral wall. Knowing
that the flow induced in the inter-gill spacing is much lower (5%) than the antero-
posterior flow, one may speculate that the two-plate model induces higher mass flow
rates. To further elucidate the effect of the flexural hinge, comparisons using only gill
pair number 4 flapping at Reω = 21.6 are shown next.
Figure 5.3 shows the time-averaged normalized velocity magnitude and stream-
lines for both models (left column). It also shows pathlines integration in proximity
of gill 4 along with time-averaged in-plane velocity vectors (right column). Part (c)
clearly shows that the velocity magnitude for the one-plate model is symmetric about
the mean position of the gill (left-right). At this Reω, the gill is still able to circulate
water, where the outflow is aligned with the gill plane and exits mostly in the dorsal
direction. On the other hand, the two-plate model in part (a) illustrates the broken
spatial symmetry where the outflow is now directed more posteriorly than dorsally.
Both models are inducing flow in the dorsal direction due to the broken symmetry
with respect to the ventral wall.
Pathlines integration shown in parts (b) and (d) also demonstrates the effect
of the hinge on the surrounding field. The specific zone shown was selected because
the time-averaged streamlines (or velocity vectors) approximately coincides for both
156
Figure 5.3: Time-averaged normalized velocity magnitude for Reω = 21.6 using gill
pair 4 and pathlines integration. Top: with flexural hinge. Bottom: without flexural
hinge. Left: normalized velocity magnitude and streamlines. Right: in-plane pathlines
colored by velocity magnitude and time-averaged in-plane velocity vectors for insets
shown in (a) and (c). (N) begin retraction. ( ) begin protraction.
157
models. Also, Pathlines are colored by the velocity magnitude. More importantly,
landmarks highlighting the beginning of retraction (N) and protraction ( ) have been
added. The advantage of the broken spatial symmetry is evident when comparing the
particle displacement during protraction (by following a particle from to N) with
the time-averaged velocity vectors. For the one-plate model the particles retrograde
displacement is 4-5 times the amount of retrograde displacement of the two-plate
model. In addition, the retrograde displacement of the one-plate model is mostly
opposite to the mean flow, while the retrograde displacement of the two-plate model
has a vertical upwards component that helps the mean flow. The latter is probably
due to the upward movement of the distal flap during protraction. For both models,
the magnitude of retrograde displacement was found to be higher for particles closer
to the gill. When comparing the net displacement after a full stroke, one notices
that the net displacement is about 2.5 and 2 numerical cells for the two-plate and
one-plate models respectively. This is about 25% increase in the net displacement of
particles.
Interestingly, the average magnitude of the velocity of the gill tip (T1) in the
two-plate and the one-plate model is about 0.75 and 0.64 respectively, which is about
17% increase. This means that the enhanced flow of the two-plate model in not only
due to the increased tip velocity and that the asymmetry introduced by the distal
flap is playing a major role. Therefore, one may conclude that the two-plate model
is more efficient than the one-plate model due to the smaller retrograde displacement
introduced by the former.
Noteworthy is that the vortex detachment illustrated in the rowing-flapping
158
parametric space (§4.4, Figure 4.6) for Reω ≥ 8.1 is still present in simulations using
only gill pair 4 and in simulations using the one-plate model for Reω = 21.6. The 4
simulations running at Reω = 2.3 exhibit the same behavior discussed in 4.4, most
notably the increased amount of oscillations with respect to the actual displacement
of the particles. The next step step now is to explore the behavior of the specific
mass flow rate ξ within this parametric space.
5.3 Quantitative Assessment
To further evaluate the performance of this parametric space, and based on
findings from §4, the specific mass flow rate is computed for all 8 simulations. Fig-
ure 5.4 shows ξ, where both axis are normalized. Recall that the normalization is
done using the total surface area of the moving gills, which means that m˙ and W˙
are simply divided by 24× the planform area when using the full array or by 4× the
planform area when using gill pair 4. In addition, the results shown are for CV2,
which is 1L away from the bounding box of the moving 6 pairs of gills. Note that the
normalization of m˙ and W˙ does not affect the numerical value of ξ.
One feature in Figure 5.4 actually matches with the rowing-flapping parametric
space, which is the increase of ξ when Reω is increased. This is equally true whether
the full gill array is used or only gill pair 4. This also holds true whether the flexural
hinge is introduced or not. Furthermore, when comparing the two-plate model with
the one-plate model when using the full gill array at Reω = 21.6, one notes that m˙
has increased by 64.31%, whereas W˙ has increased by only 3.14%. Similarly when
159
Figure 5.4: Specific mass flow rate ξ for the parametric space (normalized axis). Rate
of work done and mass flow rate are normalized by the surface area of the moving
gills. Solid lines: full array, two-plate model. Dashed lines: gill pair 4, two-plate
model. Scattered symbols: one-plate model.
160
comparing the two-plate model with the one plate model when using gill pair 4 at
Reω = 21.6, one notes that m˙ has increased by 98.10%, whereas W˙ has decreased
by 3.94%. In short, the introduction of the flexural hinge increases the average mass
flow rate by 64.31% − 98.10% with a negligible change in the average rate of work
done. This results in an increase in ξ by 59.31% and 106.23% for the full array and
gill 4 pair respectively when introducing the flexural hinge at Reω = 21.6.
Likewise, on the lower end of Reω, the introduction of the flexural hinge (two-
plate model) increases the average mass flow rate by 88.10% and 153.88% when using
the entire gill array and gill pair 4 respectively in comparison with the one-plate
model. The corresponding increase in the average rate of work done is 9.34% and
6.53% only. This yields an increase in ξ of 72.03% and 138.32% respectively.
The above quantification shows that ξ favors the two-plate model over the one-
plate model when using gill pair 4. This coincides with known physics in that the
unbroken spatial symmetry of the one-plate model is disadvantageous. However, the
one-plate model is still able to circulate water at low Re, which is probably caused
by the broken ’orientational’ symmetry with respect to the ventral wall. Based on
my work and the work of [46] and [54], I speculate that the one-plate model will still
induce a mass flow at much lower Reω as long as the orientational symmetry with
respect to the ventral wall is broken.
The above results show that when the full array of gills is used, ξ favors the
two-plate model. However, ξ shows that the introduction of the hinge when gill pair
4 is used is more advantageous than when introducing it within the full gill array.
In other words, the removal of the flexural hinge from an entire gill array is not as
161
bad as removing the flexural hinge form one pair of beating gills. This is because the
full array may still benefit from the asymmetry introduced by the metachronal wave.
Moreover, whether using the full array of gills or only gill pair 4, the lower the Reω,
the more beneficial is the introduction of the hinge. In other words lowering the Reω
will worsen the effect of removing the flexural hinge.
Overall, the introduction of the flexural hinge is beneficial for all cases consid-
ered as illustrated in Figure 5.4. In addition, the importance of the flexural hinge
becomes more pronounced at lower Reω and when using a single pair of flapping gills
in comparison with the entire array. Furthermore, ξ favors flapping with gill pair 4
rather than flapping with an entire array whether the flexural hinge exists or not. A
detailed discussion of the hydrodynamic forces acting on the gills follows.
Figure 5.5 shows the different hydrodynamic forces acting on the mayfly when
flapping at Reω = 21.6 using gill pair 4 for the one-plate and two-plate models.
Part (a) shows the tip velocity magnitude along with its u (ventral-dorsal) and v
(posterior-anterior) components. Part (b) shows the rate of viscous and pressure
work done by the gill pair (or mayfly) on the surrounding fluid. Part (c) shows the
total hydrodynamic force acting on the mayfly.
Recall that the proximal plate has the exact same kinematics in the one-plate
and the two-plate model. With reference to Figure 5.5a, one immediately distin-
guishes the different behavior of the gill tip velocity. For the one-plate model, the
symmetry is evident in both velocity components, where their signs are reversed every
≈ 0.5τ . On the contrary, in the two-plate model, the v component of the tip remains
positive for ≈ 0.6τ . While the proximal plate starts to protract at t/τ ≈ 0.575, part
162
Figure 5.5: Flapping using gill pair 4 at Reω = 21.6. (a) tip (T1) velocity (u: ventral-
dorsal, v: posterior-anterior, ‖u‖L2: velocity magnitude). (b) total pressure and
viscous rate of work done by the mayfly (1 pair of gills), gill 4 position is shown at
{0.1,0.5,0.8}. (c) total force acting on the mayfly. Dashed lines: one-plate model.
Solid lines: two-plate model.
163
of the distal plate is still moving anteriorly until t/τ ≈ 0.725, when the distal plate
completely reverses it direction and the tip starts to accelerate posteriorly. This re-
sults in the asymmetry of the velocity magnitude and delayed maxima in comparison
with the one-plate model.
With reference to Figure 5.5b, the effect of the flexural hinge on the rate of
work done is illustrated. The primary effect is observed on the rate of work done by
the mayfly to overcome the pressure forces, where the maximum rate of work done
during protraction is about one third that of the one-plate model. This simply means
a lower pressure drag on the distal plates due to their orientation during protraction.
Furthermore, in both models, negative rate of work done is observed by the end of
retraction owing to the deceleration of the gill. However, only the one-plate model
exhibits a negative work done by the end of protraction. As for the forces acting on
the mayfly (Figure 5.5c), the distinction between both models is not as pronounced
as with the rate of work done. Note the the z-component of the force acting on the
mayfly body is zero owing to the medial symmetry.
Now, a detailed look on the rate of work done by the proximal and distal plates
of left gill 4. Figure 5.6 shows the rate of work done to overcome pressure and viscous
forces for the two-plate model at Reω = 2.3 and Reω = 21.6 when flapping using gill
pair 4 or flapping using the entire array. All curves present the results when using
the two-plate model.
The first feature to notice in Figure 5.6 is how the amount of negative work done
behaves when decreasing Reω. Let us focus on the rate of work done when flapping
using gill pair 4 only (blue curves). In part (b), the fluid is not able to exert any
164
(a) Distal plate, Reω = 21.6 (b) Distal plate, Reω = 2.3
(c) Proximal plate, Reω = 21.6 (d) Proximal plate, Reω = 2.3
Figure 5.6: Rate of work done by gill 4. Left column: Reω = 21.6. Right col-
umn: Reω = 2.3. Top row: distal plate of gill 4. Bottom row: proximal plate of gill
4. Red curves: entire array. Blue curves: gill pair 4. Solid curves: pressure work rate.
Dashed curves: viscous work rate.
work on the distal plate especially by the end of protraction and beginning retraction
(at t/τ ≈ 0.65). The negative pressure work rate occurs only during ≈ 0.07τ when
Reω = 2.3 compared to ≈ 0.19τ when Reω = 21.6 (part (a)). This simply translates
from the weakening inertia of the surrounding fluid. In other words, the surrounding
fluid cannot ’store’ energy in order to send it back to the gill. Coupled with the
negligible mass of the gill, this illustrates why small organisms need to keep moving
around at small Reω in order to propel themselves or circulate water as discussed
165
in [53]. The same is true for the proximal plate shown in parts (c) and (d), where the
amount of negative pressure work rate lasts ≈ 0.15τ and ≈ 0.20τ for Reω = 2.3 and
Reω = 21.6.
Now, let us compare the work rate of the proximal and distal plates of left
gill 4 when it is flapping alone or when it is flapping within the full array of gills
(blue curves vs. red curves). Since the kinematics of gill 4 are exactly the same
in both scenarios, any differences observed between the curves are attributed to the
interaction with the neighboring gills (left gills 3 and 5). For all 4 parts, the rate of
work done is higher when gill pair 4 is present within the entire flapping array. The
use of an entire flapping array has increased the average rate of work done by the
proximal plate by 41.71% and 107.20% for Reω = 21.6 and Reω = 2.3 respectively in
comparison with flapping using gill pair 4 only. Also, the average rate of work done
by the distal plate has increased by 35.27% and 70.15% for Reω = 21.6 and Reω = 2.3
respectively. This shows that the lower the Reω, the more interacting the gills are
and therefore the higher the rate of work done.
Previously, Figure 5.4 showed that ξ favors the scenarios in which gill pair
4 is flapping alone when compared to the flapping using the entire gill array. Now,
Figure 5.6 shows that the interaction between the gills leads to higher work rates. This
simply means that the increase of the rate of work done due to increased interaction
between the gills does not yield a proportional increase in the induced mass flow rate.
Additionally, in §4, it was pointed that the flow due to inter-gill pumping contributes
with only 5% of the inflow (see Figure 4.18). The most likely outcome of all of the
above is that the presence of the entire array of flapping gills becomes increasingly
166
unfavorable at lower Reω. In specific, introducing the full array of flapping gills
decreases ξ by 50.14% and 58.73% for Reω = 21.6 and Reω = 2.3 relative to flapping
using only gill pair 4.
5.4 Last Remarks
The primary goal of this chapter was to investigate the effect of the flexural
hinge on the flow field, on the flapping kinematics and on the specific mass flow rate
ξ. This was achieved by comparing 8 simulations where one half of these simulations
uses the flexural hinge, while the other half doesn’t (Figure 5.4). The results clearly
illustrated the benefit of introducing the flexural hinge because of the broken spatial
symmetry. Furthermore, the results point to the effect of the orientational asymmetry,
where flapping without the flexural hinge is still inducing a mass flow rate at the low
Reω of 2.3. The main benefit of the flexing of the distal plate was a reduced pressure
drag during protraction while increasing the induced mass flow rate.
The investigation of the hydrodynamic forces on the proximal and distal plates
of gill 4 elucidated the interaction between neighboring gills (Figure 5.6). It was
deduced that the increased interaction between neighboring gills increases the rate of
work done by the mayfly while not yielding a proportional increase in the induced
mass flow rate. Therefore, ξ favors cases where gill pair 4 is flapping alone rather than
cases where the entire gill array is flapping metachronally. Based on ξ, the optimum
number of gills may still be more than 1 but seems to be fewer than 7.
To conclude this chapter, two more simulations were performed. These were
167
rowing with gill pair 4 only (i.e. gills 1-3 and 4-7 are removed) while using planform A
at Reω = 21.6 and Reω = 2.3. These two simulations were compared with simulations
using the entire rowing gill array previously shown in the rowing-flapping parametric
space (§4). Interestingly, the use of an entire rowing array seems to decrease ξ in
comparison with rowing using only gill pair 4. However, the lower the Reω, the lesser
is the effect of using the entire rowing array. In specific, introducing the full array of
rowing gills decreases ξ by 33.96% and 24.53% for Reω = 21.6 and Reω = 2.3 relative
to rowing using only gill pair 4.
The above trends are in contrast with the flapping results. Otherwise said,
using the entire array as opposed to using only gill pair 4 has a worsening effect on ξ
for both flapping and rowing in the investigated regime. In comparison with beating
using gill pair 4 only, the use of a flapping array has decreased ξ by 50.14% and
58.73% for Reω = 21.6 and Reω = 2.3 respectively, while the use of a rowing array
has decreased ξ by 33.96% and 24.53% for Reω = 21.6 and Reω = 2.3 respectively.
The observed trends are opposing each other. In short, while the introduction of
the full gill array decreases ξ for both kinematics, a decrease in Reω is beneficial
for the rowing kinematics but is detrimental for the flapping kinematics, all due to
higher interactions between neighboring flapping gills. Therefore, one may speculate
that there exist a certain Rearray (< 2.3) at which using an entire rowing array will
not change ξ in comparison with rowing using only gill pair 4 only. At that same
Rearray, using a flapping array will decrease ξ by more than 58.73%. Furthermore,
if one were to carry the rowing simulations of gill pair 4 using planforms B and C,
it is expected that Rearray is going to be higher than that using planform A. This is
168
because planform B and C will exhibit less interaction with other due to their smaller
areas.
In going from planform A to B to C, the aspect ratio of the gill has increased
from 1.57 to 2.23 to 2.81. Increasing the aspect ratio much further will approach that
of a cilium (≈ 50), which always uses a rowing type of kinematics. The results of ξ
suggest that rigid structures (gills in small nymphs) are unable to work cooperatively
at low Reω due an undesirable increase in their interaction yielding a much higher
rate of work done. Therefore, rowing ’independently’ using narrower gills may be the
solution. Only at a certain Rearray that the interaction of wider gills may be tolerated
and a switch to flapping kinematics is therefore favorable.
Finally, recall that all simulations were executed using prescribed kinematics.
Lowering Reω further will yield higher forces and moments that may not be withstood
by the actual gills. Hence, further evaluation of the above results will definitely require
a strongly coupled Fluid-Structure-Interaction model with properly known structural
properties of the gills.
169
Chapter 6: Conclusions and Future Work Suggestions
Mayfly nymphs are an attractive living testbed for a variety of phenomena re-
lated to the rowing-flapping dichotomy. They operate in the less studied intermediate
Reynolds number regime and abruptly switch kinematics throughout ontogeny. In-
vivo experiments on C. triangulifer have revealed a switch from antiplectic rowing to
flapping, where the latter is coupled with the development of a flexural hinge. In this
study, and building upon the earlier experimental work by [58,59], an accurate numer-
ical model was built. The numerical results for cases matching the in-vivo and robotic
experimental conditions were found to reproduce the basic flow physics within the
uncertainty of the experiments. Afterward, numerical ’experiments’ were executed,
where rowing and flapping kinematics for different gill planforms are used, for a range
of Reynolds numbers not occurring in nature. The result-driven simulations led to
two distinct parametric spaces, where the first investigates the rowing-flapping tran-
sition of the full gill array, while the second looks into the role of the flexural hinge
and the interaction between neighboring gills.
Based on the rowing-flapping parametric space discussed in §4, the following is
concluded+ or corroborated?:
• A newly introduced performance metric, ξ, is being optimized throughout the
170
rowing-flapping transition of C. triangulifer, which shows that the rowing-
flapping transition minimizes the rate of work done by the nymph per unit
mass flow rate induced throughout its ontogeny+. In other words, the switch
from rowing to flapping is most favorable or desirable by ξ.
• To optimize ξ while maintaining planform A, the rowing-flapping transition of
C. triangulifer would be 1.0 < Reω < 2.3+.
• To optimize ξ while morphing from planform A to planform B or C, the rowing-
flapping transition of C. triangulifer is 2.3 < Reω < 8.1+.
• Planforms with higher aspect ratio will delay the rowing-flapping transition, i.e.
the transition will occur at a higher Reynolds number+.
• The rowing-flapping transition is coupled with the onset of vortex detachment
for both rowing+ gills and flapping? gills. The vortex detachment was observed
for Reω ≥ 8.1.
• Within the investigated regime, the rate of work done, rate of viscous dissipation
and induced mass flow rate by metachronally rowing gills is much higher than
the rate of work done by metachronally flapping gills+.
• Contrary to locomotion in animals, the rowing kinematics in ventilating mayfly
nymphs C. triangulifer results in a higher mechanical efficiency than the flap-
ping kinematics+.
• Contrary to locomotion in flapping animals, where locomotion occurs only above
some critical Reynolds number, a net flow rate may be induced by flapping at
171
Reynolds numbers below the transition regime as long as spatial or orientational
symmetry is broken?.
• Low Reynolds number flows have very low mechanical efficiencies due to the
high rate of viscous dissipation?.
Based on the rowing-flapping parametric space discussed in §5, the following is
concluded+, corroborated? or highly speculated!:
• The introduction of the flexural hinge into a flapping gill or an array of flapping
gills increases ξ, i.e increases the induced mass flow rate per unit power+.
• The flexural hinge results in a different distribution of vorticity and instanta-
neous dipoles that increases the synergy between flapping gills+.
• The flexural hinge reduces the maximum pressure drag on the distal plates while
maintaining the same induced mass flow rate in comparison with rigid gills+.
• The lower the Reynolds number, the more beneficial is the flexural hinge?.
• To optimize ξ while morphing from planform A to planform B or C without
developing a flexural hinge, the rowing-flapping transition of C. triangulifer
would occur at Reω > 21.6+.
• In the investigated regime, ξ favors gills rowing or flapping alone rather than
gill arrays rowing or flapping in a metachronal wave+.
172
• For flapping kinematics, the lower the Reω, the more detrimental is the effect
of having an array of gills on ξ in comparison of having only one gill!.
• For rowing kinematics, the lower the Reω, the less detrimental is the effect of
having an array of gills on ξ in comparison of having only one gill!.
• Orientational asymmetry with nearby walls is allowing a symmetrically flapping
gill to induce a net mass flow rate at Reynolds numbers as low as 2.3+.
Throughout this study, it became clear that some parameters are affecting the
flow physics more than others. Future studies may include some of these parameters
as an independent variable in the parametric space. In order of importance, the
following parameters are worth looking at: 1) The inter-gill spacing with respect to
the gill length. 2) The distance from the ventral (bottom) wall with respect to the gill
length. 3) The orientation of the ventral (bottom) wall with respect to the nymph. 4)
The location of the flexural hinge. 5) The orientation and shape of the hinge. 6) The
size of the domain. 7) The number of gills. 8) The varying length of the gills within
the same mayfly. 9) Much higher and much lower Reynolds numbers. 10) Gills with
much higher aspect ratios. 11) An idealized shape of the gill (e.g. ellipse). 12) The
phase lag between neighboring gills (studied by [43]). 13) The function of the fixed
pair of gills (gill pair 7). Most of the preceding points are specific to C. triangulifer,
while only 10-12 would benefit the rowing-flapping transition in general. Point 1
is strongly affecting the interaction between the neighboring gills and is the most
important parameter if one were to explain the rowing-flapping transition within
173
an array (by looking at forces). Point 2 plays a role with orientational symmetry
and might affect the mass flow rate through the inter-gill spacing. The long list of
independent variables coupled with the variability of mayflies might actually provide
an explanation why the observed in-vivo transition occurs in a regime rather than at
a specific critical Reynolds number.
Future work should include a fluid-structure interaction model since the anal-
ysis performed on the work done by the gills in §5 is very sensitive to the extracted
kinematics, the shape of the gill and their mean locations with respect to each other.
However, the structural properties of the gills throughout ontogeny are not known and
would require substantial measurements. Future studies should also aim to incorpo-
rate oxygen transport into the computational model, which is, likewise, confronted by
an uncommon boundary condition on the surface gills, where the physiology governing
the concentration of oxygen inside the tracheal tubes is not well studied. Currently,
a portion of the data generated by the first parametric space is being investigated by
Rodolphe Chabreyrie1 using Lagrangian coherent structures and passive scalar trans-
port. This will result in new quantitative metrics that may substantiate the current
work.
1Rodolphe Chabreyrie, PhD. Department of Aerospace and Mechanical Engineering. The George
Washington University, Washington, DC.
174
RESEARCH ARTICLE
Effect of metachronal phasing on the pumping efficiency
of oscillating plate arrays
Mary Larson • Ken T. Kiger • Khaled Abdelaziz •
Elias Balaras
Received: 16 June 2013 / Revised: 28 April 2014 / Accepted: 28 April 2014 / Published online: 16 May 2014
 Springer-Verlag Berlin Heidelberg 2014
Abstract A programmable oscillating plate array was
constructed in order to study the detailed hydrodynamics of
external pumping by a series of oscillating plates at Rey-
nolds numbers on the order of 10. The array was modeled
after the geometry and kinematics found in the nymphal
mayfly (Ephemeroptera) Centroptilum triangulifer, and
consisted of five plates, each of which could be actuated
independently for stroke and pitch. Scaled tests were per-
formed at a Reynolds number, Re = fLg
2/m = 18, with a
single stroke kinematic pattern modeled after the living
animal. In mayflies, and in many other oscillating plate
systems, an antiplectic metachronal wave is used with a
phase delay of approximately 90, which corresponds to a
travelling wave that moves from posterior to anterior with a
wavelength of approximately four plates. In order to better
understand possible reasons for why the animal system
might favor the observed phase lag, ensemble-correlation
stereo PIV measurements were made to reconstruct the
unsteady three-dimensional phase averaged flow field at a
resolution that allowed a uniform and converged estimate
of the net pumped flux and the total energy dissipation
within and around the vicinity of the gill array. The results
indicate that the baseline case offered an optimal spot in
the mass flux of fluid pumped through the array per unit
energy expended, while also providing a great deal of
flexibility in modifying the stroke amplitude without
interference effects from adjacent gills.
1 Introduction
Small-scale chemical sensors and micro-reactors intrinsi-
cally rely on pumping and scalar mass transport to
accomplish their design goals (Vitko et al. 2004). Tradi-
tionally, the size and flowrate through the mechanical
components used for these functions have been large
enough that inertial convective mechanisms can be effec-
tively used, as exemplified by devices such as centrifugal
fans. One challenge faced by continually shrinking pack-
ages, however, is a concurrent decrease in the operational
Reynolds number and a transit towards viscous pumping
mechanics where such inertial mechanisms are ineffective
(Quin and Grimes 2008). A variety of micropump mech-
anisms have been studied for these applications with the
most prevalent concepts being reciprocating diaphragm
pumps (Thomas and Bessman 1975; van Lintel et al. 1988;
Bourouina et al. 1997), elecrtohydrodynamic pumps (Bart
et al. 1990), electroosmotic pumps (Jiang et al. 2002; Wang
et al. 2009), and magnetohydrodynamic pumps (Jang and
Lee 2000). Many of these pumps exhibit effective fluid
transport for particular conditions, but are limited to fluids
with specific chemical properties, finite volumes of fluid or
small flow rate ranges (Laser and Santiago 2004; Woias
2005). The development of a versatile micropump that can
be effective over a wide range of pumped volumes, flow
rates, and fluid properties in the range of Reynolds number
\100 is desirable but has not yet emerged.
One approach to exploring improved pump performance
across the transition from inertial to viscous dominated
regimes is to observe how animals at this intermediate
scale cope with similar requirements. Animals frequently
use oscillating appendages for propulsion or pumping
across the entire range of Reynolds number, from inertia-
less creeping flow (Re  1) to inviscid potential flow
M. Larson  K. T. Kiger (&)  K. Abdelaziz
Department of Mechanical Engineering, University of Maryland,
College Park, MD 20742, USA
e-mail: kkiger@umd.edu
E. Balaras
Department of Mechanical and Aerospace Engineering, The
George Washington University, Washington, DC 20052, USA
123
Exp Fluids (2014) 55:1741
DOI 10.1007/s00348-014-1741-5
175
(Re  1) (Childress 1981). Distinct patterns of appendic-
ular motion have been observed on either side of this
intermediate transition region, dividing stroke kinematics
into two main categories: flapping and rowing. Flapping is
generally defined as a motion where the net thrust is per-
pendicular to the stroke plane and is predominantly
exhibited at higher Reynolds number (Re [ 100). In con-
trast, rowing is categorized as a drag-based mechanism
with the net thrust being parallel to the stroke plane, and is
found to be prevalent at low Reynolds numbers (Walker
2002). Rowing necessitates use of asymmetry in the motion
or geometry of the appendage, generally with distinct
power and recovery strokes, while flapping can frequently
allow for symmetric motions and still generate a net thrust.
While the simplest geometry calls for a single pair of
oscillating appendages, many animals (such as ciliated
organisms, crustacea and arthropods) use appendage
arrays. This introduces limb spacing and the relative
phase difference between adjacent appendages as addi-
tional parameters for the design. For low Reynolds
numbers, cilia frequently beat out of phase when posi-
tioned in large groups as a result of hydrodynamic cou-
pling of independent adjacent oscillators (Gueron et al.
1997). A phase lag between adjacent cilia creates a
motion referred to as a metachronal wave, and is one
means to introduce a non-time-reversible parameter that is
necessary to generate flow under creeping flow conditions
(Childress 1981). Two-dimensional simulations of discrete
cilia confined within a channel have been conducted to
study the effect of the metachronal phasing and append-
age spacing under creeping flow conditions for both
symmetric (Khaderi et al. 2012) and asymmetric (Khaderi
et al. 2011) stroke motions. The results of the work show
clearly that a net pumping can be induced approximately
perpendicular to the appendages by purely symmetric
individual motions when used in conjunction with a me-
tachronal wave (thus breaking the symmetry), which for
the particular conditions studied yielded an optimal flow
rate for k/L = 4. Allowing for asymmetric stroke motions
provided an increase in the net flowrate over equivalent
symmetric conditions. It also showed at closer appendage
spacings a favorable bias of antiplectic metachrony,
where the wave propagates in the opposite direction of
the effective stroke, over symplectic metachrony, where
the wave propagates in the same direction as the effective
stroke. A computational study by Hussong et al. (2011)
developed a continuum porous media model for sym-
metric beating cilias actuated in a metachronal manner,
extending the previous works to observe the influence of
inertia. The results show that this motion is an effective
transport method for flow in a channel well into the
inertially dominated regime, also demonstrating increas-
ing flow for decreased metachronal wavelengths. At a
much higher Reynolds number (Re * 800), a study on
the locomotion of krill compared how using synchronized
or metachronal pleopod motion affected the efficiency of
propulsion. It was found that the metachronal wave is the
most efficient and also produces the highest body veloc-
ities for the krill (Alben et al. 2010). The metachronal
wave is exhibited on animals that function both in low
and high Reynolds number ranges and may be an
important contributor to effective pumping in both cases.
The current study focuses on a specific animal (the
mayfly nymph Centroptilum triangulifer) that functions
across the intermediate range of Reynolds numbers, tran-
sitioning from a viscous to an inertial dominated regime.
This animal uses seven pairs of external gill plates (see
Fig. 1) to pump water around its body in order to maintain
an acceptable oxygen concentration necessary for respi-
ration. This particular animal is interesting because pre-
vious studies on live mayflies have shown that they
exhibit a distinct shift in their appendage kinematics at a
Reynolds number (Ref = Lg
2f/m, where Lg is the root to tip
gill length, f is the frequency of oscillation and m is the
kinematic viscosity) of order 10 (Sensenig et al. 2009).
Although the flow generated by the animal has been
documented in prior studies, the limited range of behavior
of the animal prevented a detailed understanding of why
and how such a pumping mechanism might be optimized.
One specific question concerning the kinematics of the
array is what is the optimal phase delay of the adjacent
gills for a given set of stroke kinematics? In the mayfly C.
triangulifer, a phase delay of D/ = 60–90 was observed
throughout its life cycle, which corresponds to cycle
wavelength to plate length of k/Lg = 2.3. A similar phase
shift has been noted in other animals using arrays of
appendages operating in a similar Re range, such as
copepods (D/ = 64, k/Lg = 3, Van Duren and Videler
2003), remipedes (D/ = 30, k/Lg = 3.9, Kohlhage and
Yager 1994) and ctenophores (D/ = 40, k/Lg = 4.5,
Barlow et al. 1993), but no detailed explanation has been
given for why this particular phasing may represent a
favorable range of operation at these transitional condi-
tions. The previous detailed studies on cilia array provide
some guidance, but it is not immediately clear that their
two-dimensional results can be generalized to three-
dimensional conditions. Moreover, we wish to ask a
slightly different question than what has been addressed in
previous work: does this phase difference somehow rep-
resent an optimal efficiency point in terms of the cost of
pumping fluid? In order to determine the influence of
individual kinematic and geometric characteristics of a
cost-based pumping efficiency, further experiments that
allow for isolated variation of different properties are
necessary. The current study focuses on the effects of the
phase delay in between the gills.
1741 Page 2 of 13 Exp Fluids (2014) 55:1741
123 176
2 Experimental design and setup
In order to perform proposed studies, it was first necessary
to develop a device that could both replicate the nominal
kinematics of typical mayfly gills as well as extend the
range beyond what was exhibited in nature. The array of
gills used for this experiment was designed using a sim-
plified geometry compared to the live mayfly, which was
controlled using programmed micro-servomotors. The
kinematics and physical design of the gill plate array was
taken from Sensenig et al. (2009). Kinematic definitions
from this previous study can be seen in Fig. 1.
For simplicity, only the most significant kinematic
parameters were considered. Examining the stroke, pitch
and stroke plane deviation angles revealed that the stroke
range was a dominant kinematic parameter throughout the
mayfly’s life cycle (Sensenig et al. 2009). For the lower
Reynolds number conditions, the pitch range appears to
have comparable significance to that of the stroke
angle range, while for the higher Reynolds number case,
the pitch variation is comparably small. The stroke plane
deviation is the least significant factor in both cases. By
retaining stroke and pitch functionality, the device has the
ability to replicate simplified kinematics for both the high
Reynolds number case and the low Reynolds number case.
To simplify the construction, identical gill planform
shape and kinematics were used for all gills. These were
based on the kinematics of gill 4 examined by Sensenig
et al. (2009), as this was located in the mid-point of the
array. Although the living animal consistently used an in-
tragill phase delay of D/ = 60–90, the current study
focuses on effect of varying the phase delays between the
gills, which imposes limitations on the kinematics. It was
necessary to reduce the range of the stroke and the pitch to
64 % of their original values to prevent collision between
adjacent gills and control linkages for the conditions of
DU = 180 (asynchronous motion).
The robotic gill array achieves the desired kinematics
using two micro-servomotors, crank arms and sliding
linkages to control each gill (see Fig. 2). A brief description
of the design is given below, and further details can be
found in Larson (2011). The linkage is designed such that
correlated motion of both servos changes the stroke of the
gill, while a differential movement provides the pitch var-
iation. To control the stroke of the gills, the servomotors are
arranged such that the motor rotation axis is perpendicular
to the stroke plane of the gill. A steel extension arm is
attached to each servomotor and connected to a thin metal
push-rod that is embedded in the gill, passes through the
stroke extension arm, and terminates after an s-bend within
a slot of the pitch extension arm. When the servo arm
moves, the extension arm rotates, which then moves the gill
push-rod the same angular displacement as the servo arm.
The pitch is controlled by varying the difference in the
angle between the two servo arms. When the two servos
move together so that the extension arms connected to them
stay parallel, the pitch remains fixed at a = 90. When there
is a differential movement of the two arms, the metal rod in
the gill is forced to rotate to compensate for the difference,
which causes the pitch of the gill to change. Therefore the
angles of the servo arms used to control the pitch are cal-
culated relative to the positions of the stroke servo.
An important aspect of the pumping done by mayfly
gills is the effect of combining the gills into an array. For
this reason, five gills were used for the experiment to
provide three ‘‘internal’’ gills free of end effects. On a live
mayfly, there is some slight variation for the distance
between the roots of the gills, but for this simplified model
the gills are equally spaced with a root separation to gill
length ratio of 0.6. The gill plates for the robot were fab-
ricated out of transparent, UV cured acrylic (Loctite 3525).
The flow generated by the live mayfly nymph is bilat-
erally symmetric, so it was only necessary to construct a
single side of the mayfly. The gill plate array (length from
Fig. 1 Mayfly (Ephemeroptera)
Centroptilum triangulifer (left,
photo courtesy of David Funk),
as well as kinematic angle and
coordinate system definitions
for stroke motion (right, after
Sensenig et al. 2009)
Exp Fluids (2014) 55:1741 Page 3 of 13 1741
123177
gill 1 to gill 5 = 96 mm) was immersed in a small
aquarium (260 9 310 9 510 mm) partially filled with
mineral oil (depth of fluid = 200 mm), and oriented such
that the bilateral symmetry plane of the mayfly body cor-
responded with the free surface of the liquid. A simplified
body shell for the mayfly was constructed, which only
models the abdomen of the mayfly nymph to which the
gills are attached. This geometry allows for a gimbaled
motion about the gill root, which provides for variation of
the stoke plane inclination angle (b) relative to the fixed
body. In the current work, the stroke plane inclination
angle was fixed at 60.
The model was scaled using the oscillating Reynolds
number Ref = Lg
2 f/m. Considering the velocity ranges of
the servomotors and the spacing necessary to present them
from colliding, the robotic model dimensions are approxi-
mately 54 times those of the original mayfly (Lg = 40 mm).
Mineral oil with a viscosity of 175 cSt was used as the
working fluid for this experiment to achieve a Reynolds of
17.4.
Three-dimensional velocity field data was obtained
using stereo PIV, where two cameras (LaVision Imager Pro
4 M cameras, 2,048 9 2,048 pixels, 50 mm lens) equipped
with Scheimpflug mounts were placed on the same side of
the laser sheet (approximately 2 mm thick) with an
enclosed angle of 66 between them. Because of the shape
of the tank and the index of refraction between mineral oil
and air, it was necessary to construct prisms so that the
camera lenses were normal to the windows (Prasad and
Jenson 1995), giving a resulting field of view of approxi-
mately 154 9 90 mm. Data was taken at 32 planes (nine-
teen planes across a single gill) with a 2 mm spacing
between planes, in order to obtain a three-dimensional
volume of flow data around the array of gill plates. Each
image was interrogated using a multipass algorithm with a
final interrogation window size of 32 9 32 pixels with 0 %
overlap, which gave an effective in-plane spatial resolution
of approximately 2 mm. The largest uncertainty in the
measurements was due to the variability in the position of
individual gills at specific instances in the cycle (typically
\1 mm), resulting from the tolerances in the individual
linkages and the effects of friction. To minimize the impact
of this, 40 images were acquired at the same phase angle
(17 phase angles were acquired over a single cycle), pro-
cessed individually, and then median sorted to remove
outliers from the set. The N = 30 vector fields which
exhibit a minimum deviation from the median vector field
were averaged to produce the ensemble-averaged velocity
field at a given phase angle of the stroke cycle. The
uncertainty of this velocity field is taken by the uncertainty
of the mean of the ensemble, which was found to be less
than or equal to ru/N
0.5 = ±0.05 pixels or ±0.18 mm/s for
all the data, where ru is the square root of the variance of
the ensemble. See Larson (2011) for details of the analysis.
Flux data was directly calculated from this set. For spatial
velocity gradients, the three-dimensional fields were first
Fig. 2 Mayfly (Ephemeroptera)
Centroptilum triangulifer (top,
photo courtesy of David Funk),
programmable oscillating plate
array (left), and details of
linkage assembly at three
different positions in the stroke
cycle (right). The
programmable array mimics the
mayfly geometry, but can be
programmed with arbitrary
kinematics for the stroke and
pitch motion
1741 Page 4 of 13 Exp Fluids (2014) 55:1741
123
178
smoothed using an average of the neighboring cells before
the derivatives were computed using a standard second-
order finite difference. To check the consistency of the
velocity field constructed from the minimum deviation
from the median sample described above, the divergence of
the three-dimensional velocity field was calculated and
found to vary \10 % of the typical velocity gradient
magnitudes. There were isolated locations where spurious
vectors near the gill tips caused errors close to 100 %, but
these were rare and constituted \0.02 % of vectors in a
typical volume.
As a check on the above procedure, a direct numerical
simulation for the case of DU = 90 was conducted,
using the measured kinematics reported above in Fig. 3.
The Navier–Stokes equations for viscous incompressible
flow are solved on a Cartesian grid using an exact pro-
jection method. Both the advective and diffusive terms
are advanced explicitly using a third-order, low-storage,
Runge–Kutta scheme. All spatial derivatives are discret-
ized using second-order central differences on a staggered
grid (Balaras 2004). The body and the oscillating plate
array are immersed into the Cartesian grid, and boundary
conditions are imposed using an immersed boundary
formulation, which utilizes a moving least-squares (MLS)
reconstruction procedure (Vanella and Balaras 2009). The
overall numerical scheme is second-order accurate both in
space and time. The rectangular computational domain
has dimensions of 6.4Lg, 10Lg and 12.5Lg, and is dis-
cretized with 258 9 242 9 194 points in x, y, z directions
respectively (total of 12 million points). Grid stretching is
applied in the x and z directions in order to cluster points
in the region around the immersed bodies so that the
domain density is 40 cubic cells for each gill length (Lg)
around the plate array and the body (which is approxi-
mately twice the spatial resolution of the experiments).
Grid refinement tests (not shown here) showed that the
above resolution was sufficient to capture the dynamics of
the flow. All domain boundaries are modeled as non-slip
walls to mimic experimental conditions.
A comparison of the experimental and numerical results
is shown in Fig. 4. As can be observed, the typical vorticity
distribution is virtually indistinguishable between the
experiment and DNS (Fig. 4a, b). Quantitatively, however,
there are discrepancies, most notably near the core of the
vortices produced near the tips of the gills (Fig. 4c). This is
a result of the slight under-resolution of the interrogation
window size and the smoothing applied prior to calculating
spatial velocity derivatives. Since one of the main goals of
the work is to track the dissipation produced by these flow
structures, Fig. 4d shows a comparison of the temporally
varying dissipation integrated within a control volume
bounding the gills (see Sect. 3.3 for details). By this
measure, both flows track within about 10 % of each other,
with the experimental value tending to be systematically
larger than the DNS results. Thus although the values at the
very core of a vortex are underestimated, the slight over-
estimate in the larger region outside the center compensates
and reduces the net error.
Fig. 3 Stroke kinematics for the robotic mayfly, showing the stroke
and pitch program achieved for each of the 5 gills in the array for each
of the four intergill phase delays (DU) tested. Solid lines indicate the
measured motion extracted from the stereo images during testing,
while the dotted line indicates the intended program
Exp Fluids (2014) 55:1741 Page 5 of 13 1741
123
179
3 Results
The results of this experiment contain flow field data for
four different test conditions. In these four tests, the
amplitude of the stroke and pitch are kept consistent for
each case, but different phase lags are used in between each
gill. The phase lags tested are approximately 0, 90, 180
and 270. The actual stroke and pitch of each gill plate
achieved in the experiment will differ from gill to gill from
the specified motion due to variations in the tolerance of
the actuator construction.
The kinematics realized during the experiment were
documented by tracking 38 points on each gill using the
images from the PIV measurements. Even though the index
of refraction of the gill material (ngill = 1.49) was similar
to that of the oil used in the experiments (noil = 1.47), the
edge of the gill still scattered sufficient light to be visible in
the PIV image. These points were recorded for the five gills
at every phase angle and every plane where the laser sheet
intersected with the gill plates. The results of this kinematic
tracking for each test case (distinguished by their phase
delay, DU) can be seen in Fig. 3. The black circles indicate
the position commands that were sent to the servomotors.
The same pattern was sent for each gill, just at different
times, depending on the phase delay between the gills.
3.1 Phase-resolved and time-averaged velocity
The phase-resolved velocity fields within the X–Y (coronal)
planes for the four different phase delays, as well as the
corresponding mean flow, is shown in Fig. 5. For all cases,
the flow near an individual gill is fairly similar, with
negative vorticity generated at the lateral tip (-y) furthest
from the root during retraction (motion towards ?x), and a
positive vortex during protraction (motion toward -x). A
weaker vortex of opposite sign is generated near the root of
the gill closest to the medial plane (?y). The primary
difference of the flows then stems from the timing of the
interactions of these vortices with their surrounding
neighbors. For the case of synchronous gill motion
(DU = 0), all of the gills peak their circulation at the same
point in time, leading to the formation of an array of vor-
tices of the same sign, which may be viewed as a time-
varying discretized vortex sheet extending along the x-
direction. As it peaks near mid-retraction (t/T = 0.24) the
two oppositely signed sheets (positive at the gill root,
negative at the gill tip) lead to a region of uniform
streaming within the gill region. As the gills decelerate and
reverse direction (t/T = 0.47), the vorticity at the gill tip
diffuses into the outer flow and the external flow is redi-
rected toward the anterior direction. As protraction starts, a
Fig. 4 a Experimental measurement of vorticity, b DNS result for same location and time from the experiment shown in (a), c line profile of
vorticity extracted from y/Lg = -0.79, and d comparison of the net dissipation rate within the control volume for both the experiment and DNS
1741 Page 6 of 13 Exp Fluids (2014) 55:1741
123 180
Fig. 5 Phase-resolved and time-averaged flow through the coronal plane of the gill array. Streamlines are shown with the background colored by
the non-dimensional vorticity, xz/f. For the phase-resolved fields, only 4 of the 17 sampled phases are shown
Exp Fluids (2014) 55:1741 Page 7 of 13 1741
123
181
similar but somewhat weaker flow is repeated in the
opposite direction. The protraction is weaker due to the
asymmetry in the stroke velocity (see Fig. 3).
The DU = 90 is most similar to the actual animal
conditions, and the gill motion propagates as an antiplectic
metachronal wave with a wavelength of k/Lg = 2.4. As
noted by Sensenig et al. (2010), this phasing places the
retracting plate in relative isolation of its neighbors during
the retracting power stroke, creating a dominant dipole
about the root and lateral tip on the same gill, creating an
induced flow in the posterior (?x) direction. The posteri-
orly directed jet follows the metachronal wave anteriorly
during the cycle. During the weaker protraction, the plates
are in their closest proximity and interfere with one another
to prevent significant anteriorly directed backflow.
For DU = 180, the array operates in a nearly anti-
symmetric fashion (the pitch and stroke programs are not
symmetric with respect to each half-cycle, breaking the
symmetry this standing wave), with plates making their
closest approach to neighboring gills at the beginning and
mid-point of the stroke cycle. This leads to dipole pairs on
neighboring gills that produce alternating lateral and
medial directed jets, with a dominant lateral (-y) outward
direction. Finally, DU = 270 is similar to the DU = 90
condition except that the wave exhibits a symplectic
metachronism. This now places the protracting recovery
stroke in isolation, producing a weaker anterior (-x)
directed jet that travels in the posterior direction. The
stronger retraction power stroke occurs in close proximity
to its neighboring gills, preventing effectively directed
pumping.
The time-average of the 17 different equally spaced
samples in the cycle gives a representation of the net flow
through the gill array. In general the flow is seen to enter
the array from the medial plane as well as the anterior
(x \ 0) and posterior (x [ 0) regions to the left and right of
the array. The outflow is generally directed laterally away
from the array (transverse to the attachment line of the gills
along y = 0), but is biased towards anterior or posterior
direction depending on the antiplectic (posterior flow bias)
or symplectic (anterior flow bias) metachronal kinematics.
This is qualitatively in agreement with the findings of
Khaderi et al. (2012), who found that the net flow induced
within their channel was in opposition to the metachornal
wave direction. They argued this resulted from an effective
mean pressure gradient that was created by a spatial
asymmetry between the peak compression (high pressure)
and expansion (low pressure) regions near the tip of the
appendages. While there are undoubtedly similar high and
low pressure regions between the gill plates of the current
system, the comparison can only be speculative due to the
absence of an enclosing channel and the higher Reynolds
number of the current system. Finally, it is observed that
while the symmetric and antisymmetric cases should, in
principle, produce a symmetric flow, there is a slight pos-
terior-directed bias to the outflow due to the fact that the
specified kinematics of the plates are not themselves
symmetric.
3.2 Flux
In order to properly measure the effective efficiency of the
system, one needs to first consider what is being optimized.
For living systems such as a mayfly, the oscillating
appendages are presumably used to pump fresh oxygenated
water close to the skin to enhance the uptake of oxygen into
the animal (Ba¨umer et al. 2000). In this sense, the animal
would seem to benefit most by maximizing the uptake of
oxygen with a minimal expenditure of energy. To calculate
this, one would need to measure (1) the scalar absorption
rate of oxygen onto the surface of the gill plates and the
body of the animal, and (2) the amount work done by the
plates on the working fluid.
The first of these points would require tracking a scalar
with suitable boundary conditions on the plates and the far
field. Such an experiment is beyond the scope of the cur-
rent work. Instead, we assume the absorption rate would be
tied to the local volume flux of fluid through a control
volume surrounding the extents of the plate array. Figure 6
shows three of the six surfaces of the bounding control
volume used for this calculation, which shows similar
trends to the horizontal slices depicted in Fig. 5. Specifi-
cally for all cases, inflow is largely from the dorsal direc-
tion and close to the body from the anterior and posterior
direction. The outflow is directed laterally and dorsally, but
is biased toward the posterior for the antiplectic wave and
biased anteriorly for the symplectic wave. The symmetric
and antisymmetric cases are slightly biased toward the
anterior direction.
The net flux through the array was then calculated by
multiplying the time-averaged out-of-plane velocity com-
ponents on each surface by the area surrounding each
vector, then summing separately all of the flow entering
into the Control Volume (CV, defined as a subset of the
PIV interrogation domain, given by -1.8 \ x/Lg \ ? 1.9,
-1.5 \ y/Lg \ 0.6 and -0.3 \ z/Lg \ ? 1.2) and all of
the flow moving out of the CV. In the equation below, Qin,k
represents the volumetric flowrate going into one side of
the control volume with a unit normal vector in the k
direction:
Qin;k ¼
ZZ
Ukin dxidxj ffi
XNin
lin¼1
Ukin;ldxidxj ð1Þ
where i, j and k represent the directions relevant to the
surface on the volume being analyzed, and l represents the
1741 Page 8 of 13 Exp Fluids (2014) 55:1741
123
182
index over all inward pointing vectors on the surface. i and
j represent the components in this plane, while k represent
the out-of-plane component. The total volume of fluid
moving into the control volume (Qin) is then given by:
Qin ¼
XNsides¼6
l¼1
Qin;l ð2Þ
A similar expression is used for the fluid moving out of
the control volume. The results were made non-dimen-
sional by dividing the volumetric flowrate by Lg
3f, where Lg
is the gill length, 0.04 m, and f is the frequency, 1.85 Hz.
Owing to mass conservation, the sum of the flow in and out
of the CV should equal zero, which is not exactly achieved
due to the accumulated error in the measurement. The total
uncertainty of the measurement was propagated through
the flux calculations using the standard error from the
variation in the flow field calculated from the ensemble of
30 images acquired at a single phase in the cycle. This
results in an uncertainty of about ±1 % of the average of
the total flow in and out of the control volume, consistent
with the observed discrepancy.
The amount of fluid that the array in each test case
pumps is given by the average of the magnitudes of the
total flow in and the total flow out (see Table 1). Com-
paring this value for the different test cases, a phase delay
of DU = 90 produced the greatest magnitude, followed
closely by DU = 180 (Qnet,180 = 0.96Qnet,90) and then
DU = 270 (Qnet,270 = 0.94Qnet,90). The synchronous case
(DU = 0) had a distinctly lower volume flux with
Qnet,0 = 0.60Qnet,90. The marked decrease for the
synchronous case is indicative of the pumping benefit
provided by the relative isolation and close approach that
occurs for the non-synchronous cases. In that the sym-
plectic case is only 5 % less than the antiplectic case is
somewhat surprising, as the pitch program is biased to
favor the antiplectic conditions, due to the fact that the
individual gill kinematics were based on the live animal.
To answer questions of the benefit provided by using a
closely-spaced array, it would be of interest to perform
future studies with a single gill.
3.3 Dissipation
The second parameter calculated is the amount of dissi-
pation in a fluid system, which is the amount of energy lost
due to irreversible viscous shear work done on the fluid.
This is calculated by:
e ¼ l 2 oU
ox
 2
þ2 oV
oy
 2
þ2 oW
oz
 2
þ oV
ox
þ oU
oy
 2"
þ oW
oy
þ oV
oz
 2
þ oU
oz
þ oW
ox
 2#
Due to the relatively small inertia of the fluid in these
flows, the energy imparted to the fluid by the gill plate
motion is dissipated within a short distance of the gill
surface. In a companion study to the work of Sensenig et al.
(2010), direct numerical simulations of a mayfly gill array
confirm that the average work done by the mayfly gills over
a cycle is within 2 % of the total dissipation of energy over
Fig. 6 Time-averaged volume flux across a control volume bounding
the gill array (red = outflow, blue = inflow). Select three-component
velocity vectors are shown on the surface, along with their projection
onto the surface to convey the net three-dimensional flow. The
velocity magnitude is made non-dimensional by v/fLg. a DU = 0,
b DU = 90, c DU = 180, and d DU = 270
Exp Fluids (2014) 55:1741 Page 9 of 13 1741
123183
a complete cycle in a volume within 1 gill length of the
array boundary. Using this information, it is possible to
compare the average amount of work done by the gill
arrays for each test case using the time-average dissipation
within the control volume as an equivalent surrogate for the
work performed by the gills.
The average rate of power dissipated within the control
volume is calculated as:
E ¼ 1
T
ZT
0
Zz2
z1
Zy2
y1
Zx2
x1
e x; y; z; tð Þdxdydzdt
where the subscripts 1 and 2 indicates the boundaries of the
control volume used in the previous flux calculation. The
resulting average dissipations were non-dimensionalized
by qf3Lg
5, where q is the density, f is the frequency of the
cycle, and Lg is the gill length. These results are listed in
Fig. 7. Comparing the time averages of the dissipation
shows that the gills in the DU = 0 case did the least
amount of work. The DU = 180 degree test case required
the most work followed by the DU = 270 and 90 test
case, respectively.
To understand why the average dissipation exhibited the
above trend, it is instructive to observe the peak phase-
resolved dissipation for the each case, as shown in Fig. 8.
The peak dissipation is formed when each gill is per-
forming its power stroke, and the maximum value is largest
when the gill performing its power stroke is relatively close
to adjacent gills. Observing the results of the 180 degree
test case, which generated the most dissipation, the gills
beginning their power strokes are relatively close to an
adjacent gill beginning to move in the opposite direction.
When they end their power strokes, they are moving
towards a gill once again moving in the opposite direction.
These high velocities in opposite directions at close prox-
imity to each other cause high gradients, resulting in the
largest amounts of dissipation.
A similar observation can be made when comparing the
dissipation generated by the DU = 90 and 270 phase
delay cases. For the 90 case, the antiplectic wave allows
the retracting gill to be relatively far away from adjacent
gills when performing its power stroke for the 90 case,
while the symplectic wave for DU = 270 places the plates
comparatively closer during the power stroke. This causes
smaller velocity gradients for the 90 phase lag case, which
Table 1 The non-
dimensionalized volumetric
flow rate going in and out for six
sides of a control volume for the
four different test cases
These were non-
dimensionalized by dividing by
f*Lg, where f is the frequency of
the cycle, 1.85 Hz, and Lg is the
length of the gill 0.04 m
Phase lag (DU) 0 90 180 270
In Out In Out In Out In Out
Surface 1 -0.0926 0.0331 -0.0993 0.0426 -0.0377 0.0915 -0.0988 0.0546
Surface 2 -0.1330 0.1640 -0.2550 0.2770 -0.2900 0.2350 -0.2500 0.2440
Surface 3 -0.0781 0.0503 -0.1580 0.0000 -0.1510 0.0087 -0.1180 0.0648
Surface 4 -0.0364 0.0410 -0.0757 0.1090 -0.0850 0.0279 -0.0935 0.0000
Surface 5 -0.0094 0.0713 -0.0006 0.1780 -0.0046 0.2130 -0.0023 0.1970
Surface 6 -0.0130 0.0086 -0.0146 0.0026 -0.0227 0.0035 -0.0135 0.0052
Total -0.362 0.368 -0.604 0.609 -0.591 0.580 -0.576 0.565
Avg. magnitude of flow
In and Out, Qnet
0.365 0.606 0.585 0.570
Fig. 7 Non-dimensional average dissipation and net pumped flux through the control volume
1741 Page 10 of 13 Exp Fluids (2014) 55:1741
123
184
results in lower peak and total dissipation. The test case
with a phase lag of 0 resulted in the lowest dissipation for
the same reason. The gills are always equidistance from
each other and they are always moving in the same
direction. This synchronized motion results in the lowest
gradients, and therefore produces the least amount of
dissipation.
3.4 Efficiency
To determine which case exhibits the best overall perfor-
mance, a metric to describe the ‘‘pumping efficiency’’ as a
Mass-Specific Volume Flux (MSVF) is introduced. The
MSVF was calculated by dividing the net pumped volume
flux by the average dissipation (representing the work done
by the gills). Uncertainty propagation resulted in uncer-
tainties in the MSVF of ±0.001. The calculated MSVF
values are listed in Fig. 7. This calculation reveals that the
phase lags that yield the highest efficiency are 0 and 90.
With the given uncertainty, there is no statistical difference
between the efficiencies for these two phase delays. The
lowest efficiency was produced by the DU = 180, which
is about 20 % lower than the 90 phase delay efficiency.
Based on observation of natural pumping systems such
as the mayfly, one may hypothesize that the case with a
phase delay of 90 would function with the highest
pumping efficiency. While this is true, a phase delay of 0
resulted in a similar pumping cost per unit of fluid moved
over the array. A mayfly nymph uses its gills to circulate
water around its body in order to maintain the necessary
amount of oxygen in it that allows it to breathe, and hence
the absolute magnitude of the amount of pumped fluid is
also relevant. Since the amount of fluid that the array
pumps with a phase delay of 0 is 40 % lower than the
amount pumped with the array that uses the DU = 90, the
90 phase delay would be more effective for this purpose.
In order to determine if this is a substantial explanation, a
study that specifically focuses on the scalar mass transport
to the surface of the mayfly (and therefore delivers oxygen
to the body) would have to be performed. Finally, one may
also note that the smaller phase delays may be favored to
provide a greater range of flexibility in terms of stroke
amplitude. The current tests were all conducted at a com-
mon stroke amplitude that was reduced from the equivalent
animal case owing to the fact that for the antisymmetric
and symplectic cases, the gills would have collided. Having
the flexibility to increase the stroke amplitude, and pre-
sumably the net volume flux, would be a useful advantage
if additional flux were required by the animal.
4 Conclusions
In order to further explore the pumping mechanisms of a
set of mayfly gills, a two-degree-of-freedom robotic
oscillating plate array was constructed, allowing for an
examination on the effect of gill phasing on the pumping
performance of the array. Stereo PIV was used for four
Fig. 8 Phase-resolved
dissipation rates during power
stroke for a DU = 0,
b DU = 90, c DU = 180 and
d DU = 270 at time t/T = 0.24
Exp Fluids (2014) 55:1741 Page 11 of 13 1741
123185
different test cases to measure all three components of the
unsteady phase averaged velocity field over a three-
dimensional volume surrounding the array. Data was taken
using four different phase delays: the DU = 0, 90, 180
and 270, all with the same programmed stroke and pitch
amplitude. The quantitative measurements acquired for
each case allowed for the net pumping rate, flow induced
dissipation and a ratio of these two, representing a specific
flux efficiency, to be directly computed.
A number of trends emerged when examining which
phase delay resulted in the highest efficiency, highest flux
and highest dissipation. Phase delays that cause the gills to
move asynchronously (and therefore to generate counter-
rotating vortex pairs on adjacent gill tips), were observed to
generate significantly higher amounts of flux. The phase
delay with the gill performing the power stroke the furthest
from adjacent gills, the DU = 90, produced the highest
amount of flux. This makes smaller phase delays, which
allow higher stroke and pitch amplitudes to be reached
without collision, to be advantageous for applications where
high flow rates are required. Synchronized motion
(DU = 0) produces the least amount of dissipation, which
for this low Reynolds number experiment can be equated to
the amount of work done within a region near the gill array.
The total amount of dissipation increases as the phase delay
causes the space between adjacent gills with large velocity
gradients between them to decrease.
These trends show that adjacent gills must have velocity
gradients between them to achieve high flow rates, but that
if the high gradients are accompanied by close proximity to
adjacent gills, the amount of work required will increase,
and the amount of pumping will decrease. This combina-
tion causes a decrease in efficiency. For an optimized
pumping device, it would be necessary to determine the
best balance for spacing and velocity pattern to achieve
beneficial vortex pumping, while minimizing detrimentally
high dissipation.
Acknowledgments This material in this article is based upon work
supported by the National Science Foundation under Grants Nos.
CBET1067066 and CBET0730907.
References
Alben S, Spears K, Garth S, Murphy D, Yen J (2010) Coordination of
multiple appendages in drag-based swimming. J R Soc Interface
7:1545–1557
Balaras E (2004) Modeling complex boundaries using an external
force field on fixed Cartesian grids in large-eddy simulations.
Comput Fluids 33(3):375–404. doi:10.1016/S0045-7930(03)
00058-6
Barlow D, Sleigh MA, White RJ (1993) Water flows around the comb
plates of the ctenophore Pleurobrachia plotted by computer: a
model system for studying propulsion by antiplectic metachro-
nism. J Exp Biol 177:113–128
Bart SF, Tavrow LS, Mehregany M, Lang JH (1990) Microfabricated
electrohydrodynamic pumps. Sens Actuators A 21(1):193–197
Ba¨umer C, Pirow R, Paul RJ (2000) Respiratory adaptations to
running-water microhabitats in mayfly larvae Epeorus sylvicola
and Ecdyonurus torrentis, Ephemeroptera. Phys Biochem Zoo
73(1):77–85
Bourouina T, Bosseboeuf A, Grandchamp JP (1997) Design and
simulation of an electrostatic micropump for drug-delivery
applications. J Micromech Microeng 7:186–188
Childress S (1981) The mechanics of swimming and flying.
Cambridge University Press, Cambridge, MA
Gueron S, Levit-Gurevich K, Liron N, Blum JJ (1997) Cilia internal
mechanism and metachronal coordination as the result of
hydrodynamical coupling. Proc Natl Acad Sci USA 94(12):
6001–6006
Hussong J, Breugem WP, Westerweel J (2011) A continuum model
for flow induced by metachronal coordination between beating
cilia. J Fluid Mech 684:137–162. doi:10.1017/jfm.2011.282
Jang JS, Lee SS (2000) Theoretical and experimental study of MHD
(magnetohydrodynamic) micropump. Sens Actuators A
80:84–89
Jiang L, Mikkelsen J, Koo JM, Huber D, Yao S, Zhang L, Zhou P,
Maveety JG, Prasher R, Santiago JG, Kenny TW, Goodson KE
(2002) Closed-loop electroosmotic microchannel cooling system
for VLSI circuits. IEEE Trans Compon Packag Technol
25(3):347–355
Khaderi SN, den Toonder JMJ, Onck PR (2011) Microfluidic
propulsion by the metachronal beating of magnetic artificial
cilia: a numerical analysis. J Fluid Mech 688:44–65. doi:10.
1017/jfm.2011.355
Khaderi SN, den Toonder JMJ, Onck PR (2012) Fluid flow due to
collective non-reciprocal motion of symmetrically beating
artificial cilia. Biomicrofluidics 6:014106. doi:10.1063/1.
3676068
Kohlhage K, Yager J (1994) An analysis of swimming in remipede
crustaceans. Philos Trans Biol Sci 346(1316):213–221
Larson M (2011) Effect of kinematic variation on viscous pumping by
a robotic gill plate array. MS Thesis Dissertation, University of
Maryland
Laser DJ, Santiago JG (2004) A review of micropumps. J. Micromech
Microeng 14(6):35–64
Quin D, Grimes R (2008) The effect of reynolds number on
microaxial flow fan performance. J Fluids Eng 130:101101.
doi:10.1115/1.2953300
Prasad AK, Jensen K (1995) Scheimpflud stereocamera for particle
image velocimetry in liquid flows. Appl Optics 34(30):
7092–7099
Sensenig A, Shultz J, Kiger KT (2009) The rowing-to-flapping
transition: ontogenetic changes in gill-plate kinematics in the
nymphal mayfly Centroptilum triangulifer (Ephemeroptera,
Baetidae). Zool J Linn Soc 98(3):540–555
Sensenig A, Kiger KT, Shultz J (2010) Transitional ventilation
mechanisms in nymphal mayfly Centroptilum triangulifer. J Exp
Biol 213(19):3319–3331
Thomas LJ, Bessman SP (1975) Prototype for an implantable
micropump powered by piezoelectric disk benders. Trans Am
Soc Artif Organs 21:516–520
Van Duren LA, Videler JJ (2003) Escape from viscosity: the
kinematics and hydrodynamics of copepod foraging and escape
swimming. J Exp Biol 206:269–279
van Lintel HTG, van de Pol FCM, Bouwstra S (1988) A piezoelectric
micropump based on micromachining of silicon. Sens Actuators
15:153–167
Vanella M, Balaras E (2009) A moving-least-squares reconstruction
for embedded-boundary formulations. J Comput Phys 228(18):
6617–6628. doi:10.1016/j.jcp.2009.06.003
1741 Page 12 of 13 Exp Fluids (2014) 55:1741
123
186
Vitko J, Franz DR, Alper M, Biggins PDE, Brady LD, Bruckner-Lea
C, Burge HA, Ediger R, Hollis MA, Laughlin LL, Mariella RP,
McFarland AR, Schaudies RP (2004) Sensor systems for
biological agent attacks: protecting buildings and military bases.
The National Academies Press, Washington, DC
Walker JA (2002) Functional morphology and virtual models:
physical constraints on the design of oscillating wings, fins,
legs, and feet at intermediate Reynolds numbers. J Integr Comp
Biol 42:232–242
Wang X, Wang S, Gendhar B, Cheng C, Byun CK, Li G, Zhao M, Liu
S (2009) Electroosmotic pumps for microflow analysis. Trends
Anal Chem 28(1):64–74
Woias P (2005) Micropumps—past, progress and future prospects.
Sens Actuators B Chem 105(1):28–38
Exp Fluids (2014) 55:1741 Page 13 of 13 1741
123
187
Bibliography
[1] S. Alben and M. Shelley. Coherent locomotion as an attracting state for a free
flapping body. Proceedings of the National Academy of Sciences of the United
States of America, 102(32):11163–11166, Aug. 2005.
[2] G. Astarita. Dimensional analysis, scaling, and orders of magnitude. Chemical
Engineering Science, 52(24):4681–4698, Dec. 1997.
[3] E. Baba´k and O. Foustka. Untersuchungen u¨ber den Auslo¨sungsreiz der Atembe-
wegungen bei Libellulidenlarven (und Arthropoden u¨berhaupt). Pflu¨ger, Archiv
fu¨r die Gesammte Physiologie des Menschen und der Thiere, 119(9-11):530–548,
Sept. 1907.
[4] E. Balaras. Modeling complex boundaries using an external force field on fixed
Cartesian grids in large-eddy simulations. Computers & Fluids, 33(3):375–404,
Mar. 2004.
[5] M. Baltussen, P. D. Anderson, F. Bos, and J. den Toonder. Inertial flow effects
in a micro-mixer based on artificial cilia. Lab on a chip, 9(16):2326–31, Aug.
2009.
[6] D. Barlow, M. A. Sleigh, and R. J. White. Water flows around the comb plates of
the ctenophore Pleurobrachia plotted by computer: a model system for studying
propulsing by antiplectic metachronism. J. Exp. Biol., 177(1):113–128, 1993.
[7] M. F. Beatty. Principles of Engineering Mechanics: Volume 1: Kinematics.
Mathematical Concepts and Methods in Science and Engineering. Springer, 1986.
[8] N. Beratlis. Direct Numerical Simulations of Transitional Pulsatile Flows in
Stenotic Vessels. PhD thesis, University of Maryland, College Park, 2008.
[9] N. Beratlis, E. Balaras, and K. T. Kiger. Direct numerical simulations of transi-
tional pulsatile flow through a constriction. Journal of Fluid Mechanics, 587:425–
451, Aug. 2007.
188
[10] N. Beratlis, K. Squires, and E. Balaras. Numerical investigation of Magnus effect
on dimpled spheres. Journal of Turbulence, 13(July 2014):N15, Jan. 2012.
[11] S. Childress. Mechanics of swimming and flying. Press Syndicate of the Univer-
sity of Cambridge, New York, 1st edition, 1981.
[12] S. Childress and R. Dudley. Transition from ciliary to flapping mode in a swim-
ming mollusc: flapping flight as a bifurcation in Reω. Journal of Fluid Mechanics,
498:257–288, Jan. 2004.
[13] J. den Toonder, F. Bos, D. Broer, L. Filippini, M. Gillies, J. de Goede, T. Mol,
M. Reijme, W. Talen, H. Wilderbeek, V. Khatavkar, and P. D. Anderson. Ar-
tificial cilia for active micro-fluidic mixing. Lab on a chip, 8(4):533–41, Apr.
2008.
[14] M. H. Dickinson. Unsteady mechanisms of force generation in aquatic and aerial
locomotion. American Zoologist, 36(6):537–554, 1996.
[15] L. Eastham. Currents produced by the gills of mayfly nymphs. Nature, page 58,
1932.
[16] L. Eastham. Metachronal Rhythms and Gill Movements of the Nymph of Caenis
horaria (ephemeroptera) in Relation to Water Flow Proceedings of the Royal
Society of London . Series B , Containing Papers of a Biological. Proceedings of
the Royal Society of London., 115(791):30–48, 1934.
[17] L. Eastham. The rhythmical movements of the gills of nymphal leptophlebia
marginata (ephemeroptera) and the currents produced by them in water. J.
Exp. Biol., 13:443–449, 1936.
[18] L. Eastham. The gill movements of nymphal ecdyonurus venosus
(ephemeroptera) and the currents produced by them in water. J. Exp. Biol.,
14:219–228, 1937.
[19] L. Eastham. Gill movements of nymphal ephemera danica (ephemeroptera) and
the water currents caused by them. J. Exp. Biol., 16:18–33, 1939.
[20] L. Eastham. The Abdominal Musculature of Nymphal Chloeon Dipterum L. (In-
secta: Ephemeroptera) in Relation to Gill Movement and Swimming. Proceedings
of the Zoological Society of London, 131(2):279–291, Aug. 1958.
[21] C. P. Ellington. The aerodynamics of hovering insect flight. I. The quasi-steady
analysis. Philosophical Transactions of the Royal Society of London. B, Biological
Sciences, 305(1122):1–15, 1984.
[22] C. P. Ellington. The Aerodynamics of Hovering Insect Flight. II. Morpholog-
ical Parameters. Philosophical Transactions of the Royal Society B: Biological
Sciences, 305(1122):17–40, Feb. 1984.
189
[23] C. P. Ellington. The Aerodynamics of Hovering Insect Flight. III. Kinemat-
ics. Philosophical Transactions of the Royal Society B: Biological Sciences,
305(1122):41–78, Feb. 1984.
[24] C. P. Ellington. The Aerodynamics of Hovering Insect Flight. IV. Aeorodynamic
Mechanisms. Philosophical Transactions of the Royal Society B: Biological Sci-
ences, 305(1122):79–113, Feb. 1984.
[25] C. P. Ellington. The Aerodynamics of Hovering Insect Flight. V. A Vortex
Theory. Philosophical Transactions of the Royal Society B: Biological Sciences,
305(1122):115–144, Feb. 1984.
[26] E. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined Immersed-
Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Sim-
ulations. Journal of Computational Physics, 161(1):35–60, June 2000.
[27] J. H. Ferziger and M. Peric´. Computational Methods for Fluid Dynamics.
Springer London, Limited, 2002.
[28] D. T. Greenwood. Principles of Dynamics. Pearson College Division, 1988.
[29] R. Grigoriev and H. G. Schuster. Transport and Mixing in Laminar Flows: From
Microfluidics to Oceanic Currents. Annual Reviews of Nonlinear Dynamics and
Complexity (VCH). Wiley, 2012.
[30] S. Gueron and K. Levit-Gurevich. Computation of the internal forces in cilia:
application to ciliary motion, the effects of viscosity, and cilia interactions. Bio-
physical journal, 74(4):1658–76, Apr. 1998.
[31] S. Gueron and K. Levit-Gurevich. Energetic considerations of ciliary beating and
the advantage of metachronal coordination. Proceedings of the National Academy
of Sciences, 96(22):12240–12245, Oct. 1999.
[32] S. Gueron, K. Levit-Gurevich, N. Liron, and J. J. Blum. Cilia internal mech-
anism and metachronal coordination as the result of hydrodynamical coupling.
Proceedings of the National Academy of Sciences of the United States of America,
94(12):6001–6, June 1997.
[33] J. Jeong and F. Hussain. On the identification of a vortex. Journal of fluid
mechanics, 285:69–94, 1995.
[34] H. Jiang. Numerical study of the feeding current around a copepod. Journal of
Plankton Research, 21(8):1391–1421, Aug. 1999.
[35] H. Jiang. Chemoreception and the deformationof the active space in freely swim-
ming copepods: a numerical study. Journal of Plankton Research, 24(5):495–510,
May 2002.
190
[36] H. Jiang. Hydrodynamic interaction between two copepods: a numerical study.
Journal of Plankton Research, 24(3):235–253, Mar. 2002.
[37] H. Jiang. The flow field around a freely swimming copepod in steady motion.
Part II: Numerical simulation. Journal of Plankton Research, 24(3):191–213,
Mar. 2002.
[38] S. N. Khaderi, M. H. M. Baltussen, P. D. Anderson, J. M. J. den Toonder, and
P. R. Onck. Breaking of symmetry in microfluidic propulsion driven by artificial
cilia. Physical Review E, 82(2):027302, Aug. 2010.
[39] S. N. Khaderi, J. M. J. den Toonder, and P. R. Onck. Microfluidic propulsion by
the metachronal beating of magnetic artificial cilia: a numerical analysis. Journal
of Fluid Mechanics, 688:44–65, Oct. 2011.
[40] K. Kohlhage and J. Yager. An Analysis of Swimming in Remipede Crus-
taceans. Philosophical Transactions of the Royal Society B: Biological Sciences,
346(1316):213–221, Oct. 1994.
[41] P. K. Kundu and I. M. Cohen. Fluid Mechanics. Elsevier Science, 2001.
[42] P. Lancaster and K. Salkauskas. Surfaces generated by moving least squares
methods. Mathematics of Computation, 37(155):141–141, Sept. 1981.
[43] M. Larson, K. T. Kiger, K. Abdelaziz, and E. Balaras. Effect of metachronal
phasing on the pumping efficiency of oscillating plate arrays. Experiments in
Fluids, 55(5):1741, May 2014.
[44] M. A. Larson. Effect of kinematic variations on viscous pumping by a robotic gill
plate array. PhD thesis, University of Maryland, College Park, 2011.
[45] D. J. Laser and J. G. Santiago. A review of micropumps. Journal of Microme-
chanics and Microengineering, 14(6):R35–R64, June 2004.
[46] E. Lauga. Continuous breakdown of Purcells scallop theorem with inertia.
Physics of Fluids, 19(6):061703, 2007.
[47] J. Lenarcic, T. Bajd, and M. M. Staniˇsic´. Robot Mechanisms. Intelligent Systems,
Control and Automation: Science and Engineering. Springer, 2012.
[48] A. Lilienthal, A. Zell, M. Wandel, and U. Weimar. Sensing odour sources in
indoor environments without a constant airflow by a mobile robot. In Proceedings
2001 ICRA. IEEE International Conference on Robotics and Automation (Cat.
No.01CH37164), volume 4, pages 4005–4010. IEEE, 2001.
[49] G. Liu and Y. T. Gu. An Introduction to Meshfree Methods and Their Program-
ming. Springer-Verlag, Berlin/Heidelberg, 2005.
[50] B. Morton. The generation and decay of vorticity. Geophysical &amp; Astro-
physical Fluid Dynamics, 1984.
191
[51] C. S. Peskin. Flow patterns around heart valves: A numerical method. Journal
of Computational Physics, 10(2):252–271, Oct. 1972.
[52] C. S. Peskin. The immersed boundary method. Acta Numerica, 11, July 2003.
[53] E. M. Purcell. Life at low Reynolds number. American Journal of Physics,
45(1):3–11, 1977.
[54] M. L. Roper. Symmetry breaking and un-breaking in microhydrodynamical sys-
tems: swimming, pumping and bio-ballistics. PhD thesis, Harvard University,
2007.
[55] S. P. Sane. The aerodynamics of insect flight. J. Exp. Biol., 206(23):4191–4208,
Dec. 2003.
[56] K. E. Sapsford, Z. Liron, Y. S. Shubin, and F. S. Ligler. Kinetics of antigen
binding to arrays of antibodies in different sized spots. Analytical chemistry,
73(22):5518–24, Nov. 2001.
[57] A. T. Sensenig. Kinematics of the Mayfly Nymph Gill Array: An Intermediate
Reynolds Number Ventilation Pump. PhD thesis, University of Maryland, College
Park, 2009.
[58] A. T. Sensenig, K. T. Kiger, and J. W. Shultz. The rowing-to-flapping transition :
ontogenetic changes in gill-plate kinematics in the nymphal mayfly Centroptilum
triangulifer ( Ephemeroptera , Baetidae ). Biological Journal of the Linnean
Society, 98:540–555, 2009.
[59] A. T. Sensenig, K. T. Kiger, and J. W. Shultz. Hydrodynamic pumping by
serial gill arrays in the mayfly nymph Centroptilum triangulifer. The Journal of
experimental biology, 213(Pt 19):3319–31, Oct. 2010.
[60] B. M. Sumer and J. Fredsø e. Hydrodynamics Around Cylindrical Structures.
Advanced series on ocean engineering. World Scientific, 1999.
[61] P. N. Swarztrauber. A Direct Method for the Discrete Solution of Separable
Elliptic Equations. SIAM Journal on Numerical Analysis, 11(6):pp. 1136–1150,
1974.
[62] P. N. Swarztrauber. FFT Algorithms for Vector Computers. Parallel Comput.,
1(1):45–63, Aug. 1984.
[63] L. Tang, H. J. Kwon, and K. A. Kang. Studies on effect of sample circulation
on sensing performance of protein C immunosensor. In Proceedings of the Sec-
ond Joint 24th Annual Conference and the Annual Fall Meeting of the Biomedi-
cal Engineering Society] [Engineering in Medicine and Biology, volume 3, pages
1805–1806. IEEE, 2002.
192
[64] M. Uhlmann. An immersed boundary method with direct forcing for the sim-
ulation of particulate flows. Journal of Computational Physics, 209(2):448–476,
Nov. 2005.
[65] L. A. van Duren and J. J. Videler. Escape from viscosity: the kinematics and hy-
drodynamics of copepod foraging and escape swimming. Journal of Experimental
Biology, 206(2):269–279, Jan. 2003.
[66] N. Vandenberghe, S. Childress, and J. Zhang. On unidirectional flight of a free
flapping wing. Physics of Fluids, 18(1):014102, 2006.
[67] N. Vandenberghe, J. Zhang, and S. Childress. Symmetry breaking leads to for-
ward flapping flight. Journal of Fluid Mechanics, 506:147–155, May 2004.
[68] M. Vanella. A Fluid-Structure Interaction Strategy with Application to Low
Reynolds Number Flapping Flight. PhD thesis, University of Maryland, College
Park, 2010.
[69] M. Vanella and E. Balaras. A moving-least-squares reconstruction for embedded-
boundary formulations. Journal of Computational Physics, 228(18):6617–6628,
Oct. 2009.
[70] R. A. Vijayendran, F. S. Ligler, and D. E. Leckband. A computational reaction-
diffusion model for the analysis of transport-limited kinetics. Analytical chem-
istry, 71(23):5405–12, Dec. 1999.
[71] J. J. Vitko, D. R. Franz, M. Alper, P. D. Biggins, L. D. Brandt, C. Bruckner-
lea, H. A. Burge, R. Ediger, M. A. Hollis, L. L. Laughlin, R. P. J. Mariella,
A. R. McFarland, and R. P. Schaudies. Sensor Systems for Biological Agent At-
tacks: Protecting Buildings and Military Bases. The National Academies Press,
Washington, D.C., 2005.
[72] J. A. Walker. Functional morphology and virtual models: physical constraints
on the design of oscillating wings, fins, legs, and feet at intermediate reynolds
numbers. Integrative and comparative biology, 42(2):232–42, Apr. 2002.
[73] J. A. Walker and M. W. Westneat. Mechanical performance of aquatic rowing
and flying. Proceedings. Biological sciences / The Royal Society, 267(1455):1875–
81, Sept. 2000.
[74] J. A. Walker and M. W. Westneat. Kinematics, dynamics, and energetics of
rowing and flapping propulsion in fishes. Integrative and comparative biology,
42(5):1032–43, Nov. 2002.
[75] Z. J. Wang. Dissecting Insect Flight. Annual Review of Fluid Mechanics,
37(1):183–210, Jan. 2005.
[76] T. Weis-fogh. Quick estimates of flight fitness in hovering animals, including
novel mechanisms for lift production. J. Exp. Biol., 59:169–230, 1973.
193
[77] F. M. White. Viscous Fluid Flow. McGraw-Hill, third edition, 2006.
[78] T. A. Williams. A Model of Rowing Propulsion and the Ontogeny of Locomotion
in Artemia Larvae. Biological Bulletin, 187(2):164, Oct. 1994.
[79] T. A. Williams. Locomotion in developing Artemia larvae: mechanical analysis of
antennal propulsors based on large-scale physical models. The Biological Bulletin,
187(2):156–163, 1994.
[80] J. Yang. An Embedded-Boundary Formulation for Large-Eddy Simulation of
Turbulent Flows Interacting with Moving Boundaries. PhD thesis, University of
Maryland, College Park, 2005.
[81] J. Yang and E. Balaras. An embedded-boundary formulation for large-eddy
simulation of turbulent flows interacting with moving boundaries. Journal of
Computational Physics, 215(1):12–40, June 2006.
[82] J. Yang, S. Preidikman, and E. Balaras. A strongly coupled, embedded-boundary
method for fluidstructure interactions of elastically mounted rigid bodies. Journal
of Fluids and Structures, 24(2):167–182, Feb. 2008.
[83] J. Yen and J. R. Strickler. Advertisement and concealment in the plankton:
what makes a copepod hydrodynamically conspicuous? Invertebrate Biology,
115(3):191–205, 1996.
[84] C. Zhou. Accessory gills in mayflies (Ephemeroptera). Stuttgarter Beitra¨ge zur
Naturkunde A, Neue Serie, pages 79–84, 2010.
194