## SUMMARY

Swimming and flying animals generate unsteady locomotive forces by
delivering net momentum into the fluid wake. Hence, swimming and flying forces
can be quantified by measuring the momentum of animal wakes. A recently
developed model provides an approach to empirically deduce swimming and flying
forces based on the measurement of velocity and vortex added-mass in the
animal wake. The model is contingent on the identification of the vortex
boundary in the wake. This paper demonstrates the application of that method
to a case study quantifying the instantaneous locomotive forces generated by
the pectoral fins of the bluegill sunfish (*Lepomis macrochirus*
Rafinesque), measured using digital particle image velocimetry (DPIV). The
finite-time Lyapunov exponent (FTLE) field calculated from the DPIV data was
used to determine the wake vortex boundary, according to recently developed
fluid dynamics theory. Momentum of the vortex wake and its added-mass were
determined and the corresponding instantaneous locomotive forces were
quantified at discrete time points during the fin stroke. The instantaneous
forces estimated in this study agree in magnitude with the time-averaged
forces quantified for the pectoral fin of the same species swimming in similar
conditions and are consistent with the observed global motion of the animals.
A key result of this study is its suggestion that the dynamical effect of the
vortex wake on locomotion is to replace the real animal fin with an `effective
appendage', whose geometry is dictated by the FTLE field and whose interaction
with the surrounding fluid is wholly dictated by inviscid concepts from
potential flow theory. Benefits and limitations of this new framework for
non-invasive instantaneous force measurement are discussed, and its
application to comparative biomechanics and engineering studies is
suggested.

## Introduction

The challenge of quantifying the locomotive forces exerted by fluids on swimming and flying animals has attracted the interest of both biologists and engineers for many years. Early investigations of swimming and flying animals estimated locomotive forces from kinematic data, due to the lack of experimental techniques to quantitatively measure properties of the fluid wake. Quantification and validation of locomotive forces using these approaches relies heavily on semi-empirical models (for reviews, see Webb and Blake, 1985; Vogel, 1994).

Recent studies have taken an alternative approach. Newton's second and third laws dictate that flying and swimming animals generate net locomotive forces by transferring momentum into the wake. The hallmark of the momentum transfer is the creation of vortices. Hence, studies have aimed to quantify the momentum of fluid vortices surrounding animals in order to estimate the associated locomotive forces. Over the past decade, the development of the flow visualization technique digital particle image velocimetry (DPIV) has enabled researchers to realize the direct visualization of flow in the wakes of animals and to measure the corresponding velocity and vorticity fields (e.g. Drucker and Lauder, 1999; Lauder and Drucker, 2002; Nauwelaerts et al., 2005; Spedding et al., 2003; Warrick et al., 2005). These fluid dynamic data have provided information necessary for the estimation of locomotive forces based on kinematics and dynamics of the animal wake.

Several models have been proposed to evaluate locomotive forces using the data from wake measurement data. These studies typically estimate the time-averaged force rather than instantaneous forces. For example, studies have estimated the fluid forces based on measurements of the near-appendage circulation of vortices created by the animals (Dickinson, 1996; Dickinson and Götz, 1996). Locomotive forces experienced by the animal are calculated as the reaction to the momentum of vortex loops shed into the wake (e.g. Drucker and Lauder, 1999; Drucker and Lauder, 2001; Johansson and Lauder, 2004; Stamhuis and Nauwelaerts, 2005). In these cases, the momentum of the vortex is usually measured at the time instant when the vortex ring has just detached from the animal fin/wing. The time-averaged locomotive force over the stroke cycle is then determined by dividing the momentum of the shed vortex by the time duration of the stroke cycle. In other studies, the locomotive forces have been evaluated by examining the wake far downstream, which is equivalent to taking the time-average of what occurred at the site of force generation (e.g. Spedding et al., 2003) (cf. Walker, 2004).

As discussed elsewhere (Peng and
Dabiri, 2007a), viscous dissipation and vorticity cancellation
will reduce the efficiency of the momentum transfer process from 100%,
resulting in an `information loss' in the record of locomotive dynamics
contained in by wake. However, a straightforward viscous scaling argument
shows that these effects are usually negligible on the time scale of
individual stroke cycles. In particular, the distance δ over which
viscosity will act during a single stroke of duration *T*_{S}
goes as δ∼(ν*T*_{S})^{1/2}, where ν is
the kinematic viscosity of the fluid
(Rosenhead, 1963). Regions of
opposite-signed vorticity (e.g. shed from the dorsal and ventral edges of a
pectoral or caudal fin) must be within this distance
(ν*T*_{S})^{1/2} from each other in order to undergo
vorticity cancellation and the associated information loss in the wake. For
repeated swimming or flying motions at frequency *f*_{S}, the
scaling is equivalentlyδ∼
(ν/*f*_{S})^{1/2}. Hence, information
loss in the wake becomes important if the ratioδ
/*L*∼(ν/*f*_{S})^{1/2}/*L*
is of order one or larger, where *L* is the characteristic length scale
of the appendage. A 1 Hz swimming motion in water
(ν∼10^{–2} cm^{2} s^{–1})
corresponds to a characteristic viscous length scale δ∼1 mm, which
is substantially smaller than the length scales of most fish appendages
(although not necessarily small for swimming microorganisms). In air
(ν∼10^{–1} cm^{2} s^{–1}) at 1 Hz,δ∼
3 mm, which is also smaller than the length scales of most bird
appendages. Insects may have appendage length scales on this order, but will
also operate at much higher frequencies. For example, in *Drosophila*
(ν/*f*_{S})^{1/2}/*L*∼10^{–1},
indicating a limited role for vorticity cancellation during the transfer of
momentum from the animal to the wake in a single stroke. Therefore, for the
near-wake (*vis-à-vis* far downstream) studies of concern here,
we will assume no information loss between the dynamics of the animal and the
wake it generates.

Previous studies analyzing wake vorticity have found the measured time-averaged forces to be comparable with the necessary [but not sufficient (Dabiri, 2005)] lift and thrust required to sustain flying and swimming. However, time-averaged forces provide little information about the dynamics of swimming and flying. It is the instantaneous forces that dictate important dynamics of locomotion such as the trajectory, speed and efficiency of swimming and flying. In addition, these previous methods implicitly assume that the flow is steady so that the vortex momentum can be determined from the distribution of vorticity alone. For these reasons, an approach toward the task of estimating instantaneous, unsteady locomotive forces of freely moving animals was recently developed (see Dabiri, 2005; Dabiri, 2006; Dabiri et al., 2006). This model provides a method to empirically deduce unsteady swimming and flying forces based on the measurement of velocity and vortex added-mass in the animal wake. In this method, the momentum of the vortex in animal wake is evaluated as the sum of the linear momentum of fluid inside the vortex and the linear momentum of fluid surrounding the wake vortex, i.e. the added-mass of the wake vortex. Given velocity field measurements in the wake, the vortex boundary in the wake can be determined and the momentum of the wake vortex and its added-mass can be calculated, leading to a quantitative evaluation of instantaneous locomotive forces.

In the present study, this method is applied to quantify the instantaneous
locomotive forces generated by the pectoral fins of the bluegill sunfish
*Lepomis macrochirus* Rafinesque during labriform locomotion. The
overall goal of this study is to develop a case study in which the
aforementioned method can be used to analyze velocity field measurements in
order to estimate instantaneous locomotive forces in freely moving animals. We
use a new data set from pectoral fin locomotion in a well-studied fish, the
bluegill sunfish, to explore the application of this method to fish
locomotion. Due to previous work on this species using the time-averaged
vortex approach (e.g. Drucker and Lauder,
1999; Drucker and Lauder,
2000; Drucker and Lauder,
2001), good estimates of stroke averaged forces are available for
comparison to the instantaneous forces calculated here. In addition, the mean
vertical force is known, as the weight of sunfish underwater is a previously
measured quantity (Drucker and Lauder,
1999), which allows a further check on the instantaneous force
calculations. Furthermore, we aim to explore the utility of this method with a
data set of the kind that is typically available to investigators studying
animal locomotion: a time series of two-dimensional DPIV vector fields.

The outline of our approach is as follows: (1) the time-dependent velocity field of the wake in freely swimming sunfish is measured at high temporal and spatial resolution using DPIV; (2) a dynamical systems theory developed recently (Haller, 2000; Haller, 2001; Shadden et al., 2005; Shadden et al., 2006) is used to identify the boundary of the vortex wake on the two-dimensional (2-D) plane from the velocity field data; (3) a three-dimensional (3-D) approximation based on the 2-D vortex boundary is used determine the momentum of the vortex itself and its added-mass; (4) the corresponding locomotive forces acting on the pectoral fins are calculated. Unlike previous studies that estimated the time-averaged forces over the stroke cycle, this study provides detailed information on how locomotive forces evolve within the fin stroke cycle.

Like most studies involving time-averaged locomotive force estimation of swimming and flying animals, analysis in the present study is based on the flow velocity data from 2-D DPIV measurements. 3-D velocity field data with a large control volume including both the appendage and fish body would give more accurate force estimation [e.g. using the exact equations of motion derived by Noca et al. (Noca et al., 1997; Noca et al., 1999)], but such measurements are not yet currently possible. Limitations associated with the use of flow velocity data from a 2-D DPIV measurement to evaluate instantaneous locomotive forces are investigated in this study, and directions for future studies are suggested to enable more accurate force measurements.

## Materials and methods

### Fishes and experimental conditions

Bluegill sunfish *Lepomis macrochirus* Rafinesque (*N*=7,
total body length *L*=17.3±1.0 cm, mean ± s.d.) were
collected by seine from ponds in Concord, MA, USA and maintained at an average
temperature of 20°C in 40 l aquaria.

In experiments, sunfish swam individually in the center of the working area
(28 cm×28 cm×80 cm) of a variable-speed freshwater flow tank under
conditions similar to those described previously
(Drucker and Lauder, 1999;
Drucker and Lauder, 2000;
Drucker and Lauder, 2001;
Drucker and Lauder, 2003). The
sunfishes were trained to hold station in a current with a velocity of 0.5
*L* s^{–1}. At this relatively low speed, sunfish
swimming usually involves use of the pectoral fins to generate locomotive
forces (labriform locomotion). Only steady rectilinear swimming, during which
the fish maintained a speed within 5% of the flow tank's current speed, was
considered for analysis. To minimize wall effects, the fish were required to
swim near the center of the volume of the working area. Thus the flow
structures are assumed to result directly from movements of the fish pectoral
fins.

### Wake visualization and measurement

DPIV was used to visualize and measure the wake of the sunfish pectoral
fin. General details of the method are provided elsewhere
(Drucker and Lauder, 1999),
although the data generated for this paper resulted from a modified approach
to produce a temporally and spatially much higher resolution data set. The
DPIV technique provides empirical velocity field data for flow in
two-dimensional sections of the swimming fish wake (for details, see
Willert and Gharib, 1991;
Drucker and Lauder, 1999;
Lauder, 2000). A 10 W
continuous-wave argon-ion laser (Coherent Inc., Santa Clara, CA, USA) was
focused into a thin light sheet 1–2 mm thick and 10 cm wide, which
illuminated reflective silver-coated glass spheres (mean diameter 12 μm,
density 1.3 g cm^{–3}) suspended in the water (concentration≈
14 mg l^{–1}). Particle motion induced by pectoral fin activity
was recorded by imaging the laser sheet with a high-speed video camera
(Photron Fastcam-ultima APX, 1024×1024 pixels at 500 frames
s^{–1}); a second camera (Photron Fastcam-X 1280PCI,
1280×1024 pixels, 500 frames s^{–1}) synchronously
recorded a perpendicular reference view showing the position of the fin
relative to the visualized transection of the wake. In the majority of
experiments described herein, the laser was oriented to reveal the flow in the
transverse plane, which crosses the left pectoral fin
(Fig. 1) and intersects the
sunfish body perpendicularly. This provided a full field view of the pectoral
fin wake moving toward the camera, and the high temporal sampling rate
provided a detailed image of pectoral fin wake flow patterns as the wake moved
through the laser plane. Camera images of the transverse light sheet plane
were obtained through a mirror located in the flow tank, well downstream
(5–10 fin chord lengths) of the sunfish pectoral fin. A total of 55
sequences were obtained for analysis from the seven sunfish studied.
Considerable effort was made to obtain DPIV images in which the fish was
swimming as steadily as possible, not maneuvering, and in which the fin was
located in a variety of locations relative to the horizontal light sheet. In
some sequences, we maneuvered the fish into a location where the laser sheet
sliced through the fin itself so that the boundaries of vortex structures
relative to the fin as it moved during the fin beat cycle could be determined.
The sunfish pectoral fin is translucent and the small supporting fin rays do
not significantly obstruct the laser light, so no shadows were formed that
might inhibit vector calculation between the fin and the body. In other
sequences, sunfish were induced to swim upstream of the light sheet so that
the pectoral fin wake alone was imaged.

We focused the analysis in this paper on the vertical (lift or dorsal–ventral) and side (lateral) forces, for several reasons. (1) Time averaged calculations of sunfish pectoral fin forces in both of these directions are available from previous work (Drucker and Lauder, 1999); (2) mean side forces should be near zero if the sunfish is swimming steadily, and the mean lift force should balance the weight of sunfish in water [sunfish are slightly negatively buoyant (Drucker and Lauder, 1999)], providing validation of the calculated instantaneous forces; (3) the DPIV data in the transverse plane were most accurately calculated in the side and vertical directions, while thrust data estimates were not available as flow was moving through the imaging plane.

2-D velocity fields in the wake of the swimming sunfishes were calculated from consecutive digital video images (1024×1024 pixels, 8-bit grayscale) using DaVis software (LaVision GmbH, Göttingen, Germany). Image spatial cross-correlation was performed with convolution filtering iterating from a 64×64 pixel to a 16×16 pixel interrogation area with 50% window overlap, giving about 10 000 vectors per time sample. Erroneous vectors and outliers were removed automatically by removing vectors more than 2 standard deviations away from the mean of the neighbors. Vectors were smoothed with a 3×3 vector averaging filter.

For the purpose of demonstration in this paper, only one representative set of velocity data was chosen for the following analyses. Accordingly, the results shown correspond to the specific set of data analyzed, not a composite analysis of all of the measurement samples.

### Force estimation

Locomotive forces experienced by the fin were calculated as the reaction to the total momentum imparted to the vortex wake. The momentum of the vortex wake consists of two components: the linear momentum of the fluid inside the vortex and the linear momentum of the fluid surrounding the vortex that moves in the same direction as the vortex. The latter component is the added-mass of the vortex in the wake and is identical to the added-mass traditionally associated with fluid surrounding solid bodies in potential flow (Dabiri, 2006). The expression for the wake momentum and its derivation were given previously (Dabiri, 2005) [see also Dabiri et al. (Dabiri et al., 2006) for correction]. Quantification of the momentum of the wake requires not only experimental data with sufficient spatial and temporal resolution, but also an involved mathematical analysis. Thus, approximations and simplifications were made where possible.

If the wake vortex does not deform rapidly, the impulse **I** of the
fluid circulating inside the vortex can be simplified as:
(1)
where ρ is water density (1000 kg m^{–3} at 20°C),
*V*_{V} the volume of the vortex and **U**_{V} the
velocity of the wake vortex center of mass. Furthermore, the impulse of the
wake vortex added-mass can be approximated based on the added-mass tensor
**C** of the vortex, as well as the volume and velocity of each vortex as
it is formed in the wake:
(2)
The contribution of the added-mass of the vortex depends on the vortex shape
and the type of motion the vortex is experiencing (i.e. linear or angular).
Depending on whether or not bulk rotational motion of the vortex (i.e.
rotation of the principal axes of the vortex volume) is accounted for, the
velocity **U**_{V} will be either a 3×1 or 6×1 vector,
and the vortex added-mass tensor **C** is correspondingly a 3×3 or
6×6 matrix, respectively, with elements *c*_{ij}, which
are the dimensionless added-mass coefficients that relate acceleration in the
*i*th direction to the resultant forces in the *j*th direction
(where *i* and *j* can assume translation in *x*-,
*y*- and *z*-axis directions, or rotation in the *xy*-,
*xz*- and *yz*-planes. Repeated subscripts
*c*_{ii} do not indicate summation). For the purpose of
demonstration in this paper, the added-mass effect from bulk rotational motion
of the vortex is not considered here. Also, the added-mass effect is assumed
to be dominated by the *c*_{ii} added-mass coefficients
[suggesting the existence of two planes of vortex spatial symmetry
(Dabiri et al., 2006)];
therefore, the off-diagonal elements of the added-mass tensor **C** are
neglected. The total impulse **I** of the moving vortex can then be
simplified as:
(3)
where the scalar *C* represents the added-mass coefficient for linear
motion along the direction of **U**_{V}. By Newton's second and
third laws, the locomotive force exerted by the fluid on the fin is equal and
opposite to the rate at which the momentum of the vortex **I** changes due
to the interaction between the fluid and the fin:
(4)
To make the expression in Eqn 4 compatible with the format of typical wake
measurements, the time derivative can be written in terms of the finite
difference of data collected at discrete time points *t*_{j},
with *t*_{j+1}=*t*_{j}+Δ*t* and
*j=*0, 1, 2, etc:
(5)
As Δ*t* decreases to zero, Eqn 5 becomes an estimate of the
instantaneous force generated by the swimming or flying animal. AsΔ
*t* increases to *T*, the duration of the propulsive
stroke, Eqn 5 becomes an estimate of the time-averaged locomotive force. In
this study, Δ*t*=10 ms was used. For comparison, the duration of
an entire stroke cycle was approximately 600 ms. In practice, the durationΔ
*t* is only limited by the temporal resolution of the DPIV
data.

### Vortex boundary identification

It can be seen from Eqn 4 that the estimate of locomotive forces requires the determination of the physical boundary separating the vortex from the surrounding flow. The volume, velocity and the added-mass coefficient of the vortex structure all depend on identification of the vortex boundary. For steady flows, the flow structure can be generally identified by examining streamlines derived from velocity fields measured by DPIV, because these streamlines are also fluid particle trajectories. For the unsteady flows of most animal wakes, however, defining the boundary between a vortex and the surrounding flow is not an obvious task. In some simple unsteady flows, streamlines may still reveal the boundary of the vortex if plotted in a reference frame moving with the vortex (Dabiri and Gharib, 2004). While this method can be effective for studying radially symmetric vortex ring wakes such as those generated by some jellyfish [cf. fig. 9 in Dabiri (Dabiri, 2005)], squids and salps, it cannot be used to elucidate the structure of more complex wakes of other swimming and flying animals (Dabiri, 2005).

A more objective, frame-independent technique was recently developed to determine vortex boundaries in unsteady vortex flows measured empirically (Shadden et al., 2006). In this approach, the boundary of the vortex is determined by tracking the relative Lagrangian trajectories of individual fluid particles in the flow, rather than by analyzing the velocity or vorticity fields or streamlines at each time instant as in the Eulerian perspective. The method is briefly summarized here; however, the reader is directed to the papers of Shadden et al. for greater technical detail (Shadden et al., 2005; Shadden et al., 2006).

Given a time-dependent velocity field **u**(**x**,*t*), the
trajectory of a fluid particle **x**(*t*) can be determined by the
ordinary differential equation:
(6)
with given initial conditions. The flow map, which maps fluid particles from
their initial location at time *t* to their location at time
*t*+*T*, can be expressed as:
(7)
where describes
the current location of a fluid particle advected from the location
*x(t)* at time *t* after a time interval *T*. Consider
two adjacent fluid particles **x**(*t*) and
**y**(*t*)=**x**(*t*)+δ**x**(0) in the flow at
time *t*, where δ**x**(0) is infinitesimal. Their locations
after a time interval *T* are
and
. The distance between the two
fluid particles at time *t+T* is therefore:
(8)
where the second equality comes from taking the Taylor series expansion of the
flow about point **x**, and the order O of the last term is indicated in
parentheses. By dropping the small term
O[∥δ**x**(0)∥^{2}] in Eqn 8, a parameter is
introduced that represents the rate of change of the distance between two
initially adjacent fluid particles:
(9)
The parameter , which is
called the finite-time Lyapunov exponent (FTLE), measures the linearized
growth rate of the small perturbation δ**x** over the interval
*T* of trajectories starting near **x**(*t*). In other words,
it characterizes the amount of fluid particle separation, or stretching, about
the trajectory of point **x** over the time interval [*t,t+T*].

The FTLE can be used to identify boundaries in the flow that separate regions with different kinematics, such as the boundaries of vortices. Specifically, curves in the flow corresponding to local maxima of the FTLE field indicate these boundaries, because fluid particles close to the boundary and on either side of it separate much faster than other arbitrary pairs of fluid particles in the flow; therefore, higher FTLE values occur locally near the boundary. For example, consider the two fluid particles in Fig. 2, assuming that they are initially on different sides of a vortex boundary. The two fluid particles separate from each other much faster than other arbitrary fluid particle pairs that both lie on the same side of the boundary. The flow domain boundaries corresponding to ridges in a contour plot of the FTLE field are called Lagrangian Coherent Structures, or LCS (Haller, 2001; Haller, 2002; Shadden et al., 2005).

It is important to note that though the FTLE
is a function of position
variable **x** and time *t,* it is thought of as a Lagrangian
quantity since it is derived from fluid particle trajectories over the time
interval [*t,t+T*]. The absolute value |*T*| is
used instead of *T* in Eqn 10 because FTLE can be computed for
*T*>0 and *T*<0 (Fig.
2). Forward-time integration (*T*>0) reveals repelling
LCS because particles straddling this type of LCS separate faster than other
arbitrary point pairs. Backward-time integration locates attracting LCS
because particles straddling this type of LCS converge faster than other
arbitrary point pairs when advected forward in time, i.e. separate faster when
advected backward in time (Haller,
2001). The integration time |*T*| is chosen
according to the particular flow being analyzed. If a smaller integration time
is used, then less of the boundary is revealed, whereas if a longer
integration time is used, more of the boundary is revealed. Generally, if the
integration time |*T*| is sufficiently long, the
repelling and the attracting LCS usually intersect to give the boundary of the
vortex [in cases where a vortex is known to be present; cf. fig. 6 in Shadden
et al. (Shadden et al.,
2006)]. A larger integration time |*T*| also
gives LCS with higher spatial resolution. However, the choice of|
*T*| is sometimes limited in practice by the
availability of data.

The FTLE fields in the transverse plane were calculated by using the LCS
MATLAB Kit Version 1.0 (freeware download at
http://dabiri.caltech.edu/software.html)
to analyze the 2-D experimental DPIV data. A Cartesian grid was used for the
FTLE computations in the study, with uniform spacing of 0.5 mm. The flow map
at each node was calculated
with a 4th-order Runge–Kutta integration algorithm. Since the velocity
data are discrete, a 3rd-order spatial interpolation was used to provide
necessary spatial resolution. An integration length|
*T*|=200 ms was used due to the nature of the very short
stroke cycle (about 600 ms); this limited availability of velocity data
prevented a calculation with a larger integration time. Longer integration
times result in more clearly defined maxima in the FTLE field, i.e. more
well-defined vortex boundaries. Positive and negative integration time
intervals were used to determine forward- and backward-time FTLE fields,
respectively. Repelling and attracting LCS were determined by locating ridges
(i.e. contours of local maxima) of the forward- and backward-time FTLE fields,
respectively.

### Model approximations

The boundary of the forming wake vortex in the transverse plane was determined by using the method described above. This 2-D boundary of the vortex was used to approximate the following quantities.

First, the 2-D vortex boundary on the transverse plane was used to estimate
the vortex volume *V*_{V}. It has been shown in a previous
study (Drucker and Lauder,
1999) that an isolated vortex ring is generated by the pectoral
fin of the bluegill sunfish during labriform fin kinematics. However, the
vortex boundary may not possess a simple geometry. To approximate the 3-D
shape of the vortex, the calculated vortex boundary on the transverse plane
was mapped into an ellipse, which represents the cross-section of the vortex
in the transverse plane. The calculated vortex boundary and the model ellipse
have the same long axis length and also the same width
(Fig. 3). Assuming the
elliptical vortex cross-section was perpendicular to and symmetrical about the
transverse plane, the volume of the vortex can be approximated by:
(10)
where *w* is the width of the wake vortex on the laser plane view and
*A* is the projected area of the vortex ring on the plane perpendicular
to the laser plane. The added-mass coefficient *C* of the vortex body
can be determined from an equivalent solid body calculation of the
approximated 3-D ellipsoidal vortex shape
(Lamb, 1932;
Dabiri, 2006).

The motion of the LCS was approximated by the motion of the model ellipse.
The long axes of the LCS and the model ellipse are always parallel; hence, the
model is consistent with the known motion of a vortex ring in the direction
normal to its long axis (Lamb,
1932). The projection of the vortex body velocity
*U*_{V} on the transverse plane was approximated by the
velocity of the centroid of the fluid inside the 2-D vortex boundary. Since
the velocity component perpendicular to the transverse plane was not
available, the locomotive forces in that direction were not evaluated.

### Uncertainty analysis

The uncertainty in this study comes from velocity measurements using DPIV,
interpolation of discretized velocity data when calculating FTLE field,
identification of the vortex boundary from the FTLE field, and the 3-D
approximation of the vortex structure based on the boundary information on a
2-D cross-sectional plane. The velocity measurement error using DPIV is on the
order of 1–3% for flows at these velocities and accelerations
(Willert and Gharib, 1991).
This error includes uncertainty due to the fact that the seed particles in the
flow are slightly negatively buoyant. Adrian
(Adrian, 1995) shows that the
characteristic particle lag time is a function of the fluid viscosity ν,
particle diameter *d*_{p}, and particle-fluid density ratio
.
For the present experiments, τ_{lag}=41 μs, or 2% of the
temporal frame spacing (i.e. 1 frame/500 frames s^{–1}=0.002 s)
used to process the DPIV data.

Given that the FTLE is calculated from particle trajectories (integration of velocity), a reasonable concern is that local error in velocity measurement may accumulate upon integration in time. However, it has been rigorously shown (Haller, 2002) that even large local (in space and time) velocity errors do not prevent reliable calculations of the position of LCS, as long as the errors remain small in a special time-weighted norm.

Hence, the major source of uncertainty arises from identification of LCS from the FTLE calculations and from 3-D approximation of the wake geometry based on the 2-D measurements. An uncertainty of 1 mm is approximated for measurements of position and dimensions of the LCS. The uncertainty in parameters derived from these fundamental measurements is determined by the rules of error propagation, i.e. (11) and (12) where σ is the uncertainty (Taylor, 1997).

## Results

### Vortex flow patterns

A stroke cycle of labriform locomotion for the *Lepomis macrochirus*
individuals tested lasts approximately 600 ms at 0.5 *L*
s^{–1}. It involves an oscillatory cycle of (1) anteroventral
fin movement (downstroke, *t*=0–250 ms, with time *t*=0
corresponding to the initialization of fin movement in downstroke), (2)
rotation of the fin around its long axis at the end of downstroke (stroke
reversal, or fin `flip', *t*=250–300 ms), (3) posterodorsal fin
movement (upstroke, *t*=300–550 ms), and (4) a kinematic pause
period during which the fin is held flush against the body
(*t*=550–600 ms).

The velocity field in the transverse plane and the corresponding vorticity component calculated from the velocity field are plotted in Fig. 4 at four instants during the stroke cycle. In the transverse plane, a discrete pair of counter-rotating vortices is seen during the downstroke. One vortex core is located near the dorsal edge of the pectoral fin while the other is near the ventral edges of the fin; the cores possess vorticity of opposite signs (Fig. 4A). A jet-like flow passes between the vortex pair. At the end of the downstroke, the vortex pair appears to contract, giving a vortex structure with a smaller diameter (Fig. 4B). As the upstroke begins, the vortex generated in the downstroke sheds from the fin while a new vortex becomes visible at the dorsal tip of the fin (Fig. 4C). Later in the upstroke, the vorticity map becomes complicated and it is difficult to identify a coherent vortex structure (Fig. 4D).

### Lagrangian Coherent Structures (LCS) and vortex boundary

In Fig. 5 color contour
plots of the finite-time Lyapunov exponent (FTLE) fields computed with
integration times of *T*=–200 ms and *T*=200 ms at the
arbitrary time *t*=250 ms are superimposed on the velocity map at the
same instant in time. The ridges of high FTLE values in each plot indicate
LCS. For Fig. 5A, the FTLE is a
backward-time FTLE since *T*<0, showing an attracting LCS; and for
Fig. 5B the FTLE is a
forward-time FTLE since *T*>0, showing a repelling LCS.

The FTLE fields shown in Fig.
5 do not give a vortex boundary that is as sharply defined as in
previous studies of isolated vortex rings [cf. fig. 6 in Shadden et al.
(Shadden et al., 2006)]. This
is due to the aforementioned limitation in integration time|
*T*|. Nonetheless, the location of the LCS boundary can
still be approximated to lie at the centerline of the FTLE contours in each
frame. The backward- and forward-time LCS derived from
Fig. 5A,B are plotted together
to give the vortex boundary in Fig.
6. A best-fit spline connection was used if the backward- and
forward-time LCS did not intersect in a given frame. Ideally the attracting
and repelling LCS should always completely intersect to give a well-defined
wake vortex boundary when a wake vortex is present. However, due to the
aforementioned short integration time, in some frames less of the boundary is
revealed by the LCS. In Fig. 6,
the boundary of the vortex at the transverse plane is superimposed on the DPIV
velocity field data at *t*=250 ms. Notice that it is impossible to
define a vortex boundary from inspection of the velocity field alone, whereas
the theory governing the LCS ensures that the boundary is captured by the
present FTLE measurements (Shadden et al.,
2005; Shadden et al.,
2006).

The boundary of the vortex on the transverse plane does not have a regular elliptic shape, a possibility anticipated in the previous section. The fin (the portion with high brightness alongside the sunfish body in Fig. 6) can be seen embedded inside the vortex, i.e. the shape of the vortex corresponds to the shape and location of the fin. This result that the fin is enclosed in the vortex is in accordance with the fact that the pair vortices of the vortex ring are generated at the dorsal and the ventral edges of the fin. Given the fact that the fin density is only 10% higher than water and its thickness is on the order of 100 μm (Alben et al., 2007; Lauder et al., 2007), its contribution to the total mass of the vortex can be neglected in the subsequent force calculations.

### Time evolution of the vortex

The FTLE fields were calculated at a series of discrete time instants
during the entire fin stroke cycle (*t*=0–600 ms), with time
interval between consecutive time instants Δ*t*=10 ms. The LCS
could not be determined for the very early part of the downstroke
(*t*=0–100 ms), when insufficient DPIV data is available to
determine the backward time structure (i.e. no data is available before
*t*=0, preventing knowledge of fluid particle behavior for backward
integration times *T* such that
*t*–|*T*|<0). Neither was the LCS
determined for the upstroke (*t*=400–600 ms), during which
increasing three-dimensionality of the wake limited the fluid particle
behavior that could be deduced from DPIV data in the 2-D transverse plane. The
LCS was calculated on the time span from *t*=100–400 ms, during
which an isolated vortex pair can be clearly seen on the vorticity map. This
time span covers most of the downstroke and the stroke reversal. The time
evolution of the vortex boundary is plotted in
Fig. 7. It can be seen that the
shape of the vortex boundary changes with time. Compared with the early
portion of the downstroke, the projection of the vortex on the transverse
plane becomes wider, suggesting that the pair of vortices moves closer to each
other during the late downstroke, or that the vortex pair is advected out of
the plane of the laser sheet, or a combination of these two effects.

The position of the projection of the vortex centroid on the transverse
plane was determined by calculating the centroid of the area enclosed in the
vortex boundary. The trajectory of the projected vortex centroid on the
transverse plane is plotted in Fig.
8 and the change of the position with time is plotted in
Fig. 9. The movement of the
vortex is consistent with the fin kinematics. The velocity of the projected
vortex centroid on the transverse plane was calculated from the trajectory by
the change of position in unit time (Fig.
10). It can be seen from the velocity profile that during the
downstroke, the vortex accelerates vertically most significantly during
*t*=100 to 150 ms and accelerates horizontally most significantly
during *t*=160–210 ms. It will be shown later that this pattern
of movement plays a major role in the force generation.

### Vortex geometry

The LCS shows the boundary of the vortex on the 2-D transverse plane where PIV data were taken. As mentioned previously, the 3-D shape of the vortex was derived from the 2-D boundary. Since the boundary has an irregular shape, it was approximated by an ellipse to construct the 3-D vortex shape. The time evolutions of the volume, the width and the area of the vortex (all defined previously and shown in Fig. 3) are plotted in Fig. 11A–C. The volume of the vortex remains relatively constant during the early part of the downstroke and begins to decrease significantly in the latter part of the downstroke and the stroke reversal. This effect, which may be an artifact of planar measurement of a 3-D flow, is examined more closely in the Discussion.

### Vortex added-mass

The added-mass coefficient of the vortex was determined by the approximated 3-D shape of the vortex ring (Lamb, 1932). The relationship between added-mass coefficient and time is plotted in Fig. 11D. In the early part of the downstroke, the vortex resembles a thin disk, giving a larger added-mass coefficient, while in the latter part of the downstroke, the vortex is more sphere-like, giving an added-mass coefficient with a value approaching 0.5, which is the added-mass coefficient of a sphere. Since the added-mass coefficient is determined by the shape of the vortex, it follows the change in vortex shape closely.

### Locomotive force

The forces exerted by the fluid on the pectoral fin over the time span of interest are plotted in Fig. 12, with horizontal component (lateral force) and vertical component (lift force) separated. Also plotted are the time-averaged forces calculated using the vorticity method (Drucker and Lauder, 1999). The positive horizontal force is directed toward the body while the positive vertical force is directed upward. The data indicate that at the early downstroke there is a relatively large lateral force in the direction opposite to the fin movement. The fin also generates a significant lift at the early downstroke. In the late part of the downstroke, the horizontal force is close to zero and the magnitude of the vertical force also becomes much smaller than it was in the early part of the downstroke.

There is a phase difference between when lift force and lateral force are
generated during the downstroke. Comparing
Fig. 12A with
Fig. 12B, it is apparent that
the fin generates significant lift (at *t*=100–160 ms) before it
generates lateral force (at *t*=160–210 ms). This effect arises
because the acceleration of the vortex in the vertical direction occurs prior
to its acceleration in the horizontal direction
(Fig. 10).

## Discussion

In this paper we present a case study to illustrate a framework for combining traditional DPIV measurements with a new class of Lagrangian analysis tools that analyze the wake as a distinct fluid structure with its own added-mass dynamics. The fin of the bluegill sunfish was observed to be embedded within this vortex wake structure. It is expected that a similar phenomenon will observed for other animal appendages moving in air or water. This result suggests that the dynamical effect of the vortex wake on locomotion is to replace the real animal fin with an `effective appendage,' whose geometry is dictated by the forward- and backward-time Lagrangian Coherent Structures (LCS) and whose interaction with the surrounding fluid is wholly governed by inviscid concepts from potential flow theory (Peng and Dabiri, 2007a; Peng and Dabiri, 2007b). A similar concept, the `displacement thickness', is well known in the theory of steady aero- and hydrodynamic flows (e.g. Rosenhead, 1963).

The benefit of this perspective based on an `effective appendage' is that
it facilitates a direct correlation between the morphology and kinematics of
the real appendage and the dynamics of locomotion. The expression for the
locomotive force (Eqn 4) indicates that there are three parameters
contributing to the force generation: the vortex added-mass coefficient
*C*, volume *V*_{V}, and velocity **U**_{V}
of the wake vortex body. Of these three parameters, *C* and
*V*_{V} are more directly related to the geometry of the vortex
while **U**_{V} is more closely related to the fin motion. Since
locomotive forces are generated proportionally to the change of momentumΔ
[(1+*C*)*V*_{V}**U**_{V}] in unit
time, the contribution of each parameter can be evaluated by comparing the
three logarithmic terms on the right-hand side of the following equation:
(9)
Fig. 13 shows the change of
the logarithm of each individual parameter and the total change in wake
momentum. It can be seen that early in the stroke cycle the changes in the
shape of the vortex are small and the total change in momentum follows closely
the change in vortex velocity. As mentioned previously, the vortex velocity
closely follows the fin motion during the stage of the stroke cycle.
Therefore, force production during the early fin motion is primarily dependent
on the fin kinematics. After the initial fin motion, force generation is
primarily dictated by changes in the vortex shape; hence, fin morphology
(which dictates vortex shape) governs force production later in the fin stroke
cycle. Analyses of the vertical and horizontal forces are consistent, showing
that the shape of the vortex remains relatively constant early in the stroke
cycle but changes in the latter stages. Hence we conclude that locomotive
forces can be generated not only by accelerating or decelerating the fin, but
also by changing the shape of the vortex *via* fin morphology. The
clarity of the parameter relationships in
Fig. 13 is a benefit of the
`effective appendage' approach taken presently. Since the timing of forces on
the `effective appendage' is coincident with those on the real appendage (i.e.
instantaneous forces are being evaluated), this perspective is compatible with
investigation of muscle dynamics and activity patterns
(Lauder et al., 2007;
Mittal et al., 2007).

The ability of Lagrangian methods as implemented here to identify the boundary of this `effective appendage' in a fully unsteady flow without direct appeal to the vorticity field makes it possible to develop new models of fluid dynamic locomotion that are simultaneously more accurate and less complex than existing ones. Furthermore, the current heavy reliance on vorticity and other Eulerian field concepts, which are known to be inadequate for vortex identification in unsteady flows (Haller, 2005), is made unnecessary. The present results clearly illustrate that, as suggested by Haller (Haller, 2005), the correlation between the spatial distribution of vorticity and the actual vortex boundary can be quite poor in an unsteady flow, especially where global rotational motion is present. The downstroke of the bluegill sunfish pectoral fin is dominated by such rotational kinematics.

In this study, locomotive force is quantified only on the downstroke and early portions of the stroke reversal, when the vortex structure is primarily attached to the fin. The forces estimated are consistent with the magnitude of the time-averaged forces calculated using the same vorticity method as in Drucker and Lauder (Drucker and Lauder, 1999). The time-averaged vertical and lateral forces are 2.23 mN and 4.55 mN for this study (Fig. 12), which are consistent with the time-averaged forces quantified in Drucker and Lauder (Drucker and Lauder, 1999) for the same species (3.24 mN vertical and 6.96 mN lateral), noting that the data used in that study were taken from a different animal sample and using a different experimental arrangement. The evolution of forces is consistent with the motion of the animal. The lateral forces from the paired pectoral fins are expected to approximately cancel, resulting in the absence of substantial lateral motion as was observed empirically. The resultant vertical force, i.e. lift minus weight, is positive during the early phase but negative during the later phase, consistent with the kinematic measurements of Gibb et al. (Gibb et al., 1994). Since Gibb et al. only present kinematic descriptions of the locomotion (i.e. no force measurements) (Gibb et al., 1994), the comparison that can be made with the present work is only qualitative. Nonetheless, the force peaks in Fig. 12 are suggestive of a rise during the early phase of the fin beat and sinking during the late phase, as seen by Gibb et al. (Gibb et al., 1994). Perhaps coincidentally, the evolution of the locomotive forces resembles the pattern of forces generated by an insect wing, which has an early peak when the wing starts from rest followed by decay to a stable level (Birch and Dickinson, 2003).

The greatest challenge in this study was to identify the 3-D boundary of the vortex and its motion given 2-D planar DPIV data. As seen from Eqn 4, only by identifying the 3-D boundary of the vortex can the volume, added-mass and velocity of the vortex body be determined to evaluate locomotive forces. Since only 2-D velocity data are available in the present study, approximations were required.

A concern in the calculations is the assumption that the vortex ring is symmetrical about the transverse laser plane. This may not be true during the entire stroke cycle, and an asymmetrical distribution of the vortex ring on either side of the laser plane would cause underestimation of the vortex volume. Since the vortex cannot move far from the laser plane as long as it is attached to the fin, which is oscillating near the laser plane, the approximation used in this study of the attached vortex wake is reasonable. However, after the vortex is shed from the fin, it is advected out of the laser plane by the ambient current, making the approximation invalid. This is likely the effect observed toward the end of the present measurements (e.g. Fig. 11A,C).

Attempts were made to increase the accuracy of the 3-D vortex boundary approximation. FTLE fields were calculated on horizontal planes from DPIV data, with the aim of identifying the vortex boundary corresponding to the transverse plane measurement and thereby constructing the 3-D vortex boundary. However, boundaries were not clearly revealed in the corresponding FTLE fields due to the effect of net ambient flow on these two planes, which dominated the vortical motions. Even if it is possible to identify the boundary of the vortex on three perpendicular planes, approximation is still required to construct the 3-D boundary of the vortex (though the result may exhibit higher accuracy than the present methods when measuring more complex wake geometries). The identification of 3-D vortex boundaries would also enable evaluation of the effect of rotational added-mass, which is neglected in the present study. The determination of the 3-D vortex boundary requires the development of flow visualization techniques or numerical methods that can provide 3-D velocity field information, which can then be analyzed using the LCS method implemented in this study.

Given volumetric DPIV data for swimming or flying animals (or corresponding data computed in numerical simulations), a true validation of non-invasive force measurements demands a comparison of the animal body trajectory predicted by the force measurements with the body trajectory measured empirically. To our knowledge, such a comparison of measured and predicted kinematic data has not yet been performed by any study. Ideally, wake vortex dynamics and body kinematics should be measured simultaneously and in three dimensions. A demonstration that the estimated forces agree with the time-averaged force required to sustain lift of neutral buoyancy is necessary, but not sufficient by itself to validate instantaneous force estimates (cf. Dabiri, 2005). An additional avenue for validation would be comparison of calculated force profiles with those resulting from computational fluid dynamic analysis of the same fin beats using measured 3-D kinematics from that beat (with coupled fluid-structure interactions included in the computation for the swimming case, due to the comparable density of the fluid medium and the appendage). Once the force model has been validated, other aspects of animal behavior such as energetics can be examined according to established models (Schultz and Webb, 2002).

Finally, the relative importance of unsteady effects in this study was
deduced *a posteriori via* direct examination of the contribution from
wake vortex added-mass. Alternatively, the wake vortex ratio *Wa*
introduced in Dabiri (Dabiri,
2005) and refined in Dabiri et al.
(Dabiri et al., 2006) provides
an *a priori* indication of the need to examine a particular flow with
the level of scrutiny applied here.

## List of symbols

- A
- area
- C
- added-mass coefficient for linear motion along the direction of
**U**_{V} - C
- added-mass tensor of the vortex
*d*_{p}- particle diameter
- DPIV
- digital particle image velocimetry
*F*_{L}- locomotive force
*f*_{S}- stroke frequency
- FTLE
- finite-time Lyapunov exponent
- I
- impulse
- L
- length
- LCS
- Lagrangian Coherent Structures
- t
- time
*T*_{S}- stroke duration
**U**_{V}- velocity of the wake vortex center of mass
*V*_{V}- volume of the vortex
- w
- width
- δ
- distance
- ϕ
- Lagrangian flow map
- ρ
- water density
- σ
- uncertainty
- ν
- fluid viscosity

## ACKNOWLEDGEMENTS

The authors thank Shawn Shadden, Jerry Marsden and Parviz Moin for discussions of the LCS theory. This work was supported by NSF grant OCE0623475 to J.O.D.; an ONR-MURI Grant N00014-03-1-0897 on fish pectoral fin function, monitored by Dr Thomas McKenna and initiated by Dr Promode Bandyopadhyay; and NSF grant IBN0316675 to G.V.L.

- © The Company of Biologists Limited 2007