PHYSICS OF FLUIDS 17, 113601 !2005"
Numerical experiments on flapping foils mimicking fish-like locomotion
Paolo Blondeaux, Francesco Fornarelli, and Laura Guglielminia!
Department of Environmental Engineering, University of Genova, Via Montallegro 1, 16145 Genova, Italy
Michael S. Triantafyllou
Department of Ocean Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue,
Cambridge, Massachusetts 02139
Roberto Verzicco
DIMeG and CEMeC, Polytechnic of Bari, Via Re David 200, 70123 Bari, Italy
!Received 9 February 2005; accepted 26 September 2005; published online 16 November 2005"
The results of numerical experiments aimed at investigating the topology of the vortex structures
shed by an oscillating foil of finite span are described. The motion of the foil and its geometry are
chosen to mimic the tail of a fish using the carangiform swimming. The numerical results have been
compared with the flow visualizations of Freymuth #J. Fluids Eng. 111, 217 !1989"$ and those of
von Ellenrieder et al. #J. Fluid Mech. 490, 129 !2003"$. The results show that a vortex ring is shed
by the oscillating foil every half a cycle. The dynamics of the vortex rings depends on the Strouhal
number St. For relatively small values of St, the interaction between adjacent rings is weak and they
are mainly convected downstream by the free stream. On the other hand, for relatively large values
of St, a strong interaction among adjacent rings takes place and the present results suggest the
existence of reconnection phenomena, which create pairs of longitudinal counter-rotating
vortices. © 2005 American Institute of Physics. #DOI: 10.1063/1.2131923$
I. INTRODUCTION
Fish usually propel themselves in water by means of a
transverse wave that moves backwards along the body from
the head to the tail. This propulsive mode can be broadly
divided into two classes: the anguilliform mode and the carangiform mode. Fish using the anguilliform mode are flexible over the whole body and the propulsive wave has an
amplitude that is almost constant along the body. In the carangiform mode, the amplitude of the wave becomes significant only in the posterior part of the fish, the remainder of
the body being relatively rigid. Hence, thrust generation is
confined almost exclusively to the caudal fin, which can be
considered to oscillate in an otherwise undisturbed free
stream.
An analysis of the tail kinematics shows that, in order to
investigate the propulsion mechanism, the tail can be approximated by a foil that moves with a constant forward
speed, oscillates in the transverse direction, and simultaneously rotates, the two periodic motions being characterized
by a 90° phase shift.
To understand fish propulsion, a great deal of work has
been devoted to the study of the dynamics of oscillating
foils. The papers of Triantafyllou et al.1 and Mittal2 provide
an excellent review of the experimental and numerical investigations, respectively. On the basis of two-dimensional experiments, Anderson et al.3 noted that the production of a
significant leading-edge vortex, which travels downstream
and interacts with the developing trailing edge vortex, is important for high efficiency. This conclusion has been supported by the recent two-dimensional numerical simulations
a"
Electronic mail: [email protected]
1070-6631/2005/17"11!/113601/12/$22.50
of Wang,4 Lewin and Haj-Hariri,5 Pedro et al.,6 and
Guglielmini and Blondeaux.7 In particular, Guglielmini and
Blondeaux7 obtained results for a wide range of parameters
#in particular, the Strouhal number St ranged between 0.1
and 0.6, the ratio between the amplitude of the heaving oscillations and the foil chord ranged between 0.1 and 1.1, the
phase between the pitching and heaving oscillations was varied between 65° and 105°, and the distance of the pitching
axis from the leading edge of the foil made dimensionless
with the foil chord fell in the range !−0.25, 1"$ and identified
those giving rise to optimal propulsive performances, even
though only moderate values of the Reynolds number were
considered #in particular, the Reynolds number fell in the
range !1000, 3000"$.
In the two-dimensional approximation, only the crossstream component of the vorticity is taken into account.
Therefore, the propulsive efficiency is overestimated, due to
the underestimation of the wake energy loss. A threedimensional analysis is clearly required for more accurate
and realistic estimates of the propulsive efficiency. Moreover, a three-dimensional approach able to handle different
foil geometries should clarify the apparent evolutionary convergence on the fish tail toward a lunate form.
A recent contribution to the understanding of the flow
structure behind an oscillating foil of finite span is the work
of von Ellenrieder et al.,8 where the overall structure of the
flow is examined in detail. von Ellenrieder et al.8 used flow
dye visualizations and tried to identify the vortex structures
shed in the flow by looking at the time development of streak
lines. Even though the approach is certainly valuable to determine the gross features of the flow, it has a few limitations. As discussed by von Ellenrieder et al.,8 dye filaments
17, 113601-1
© 2005 American Institute of Physics
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-2
Phys. Fluids 17, 113601 "2005!
Blondeaux et al.
contain information that starts at the point of introduction
and is integrated along the path of the dye lines. Hence, it is
difficult to ascertain whether the shape of a streak line is the
result of fluid distortion or flow memory, and a great deal of
caution must be exercised in interpreting its behavior.9 For
example, Hama10 analyzed the case of a shear flow perturbed
by an unamplified traveling sinusoidal wave and showed that
the streak line that emerges from the critical layer always
appears to roll up regardless of the true nature of the superposed perturbation. Moreover, von Ellenrieder et al.8 pointed
out that the Schmidt number of the dye was O!103" and,
hence, it is likely that the dye followed the vorticity only
during the initial stage of the flow development.
Quantitative information on the wake structure behind a
thrust-producing flapping foil of finite span began to be collected with use of the particle image velocimetry !PIV"
technique,11,12 which allows the measurement of the velocity
field to be made along with the evaluation of its phaseaveraged and random components. Moreover, the strength of
the vortices shed by the foil can be estimated too.
Presently, we investigate the flow that is generated by an
oscillating finite-span foil by means of a numerical approach.
One popular method that is available to compute the
unsteady flow generated by arbitrary bodies moving at high
Reynolds numbers is the potential flow panel method. A review of the various panel methods is given by Katz and
Plotkin.13 Recently, Vest and Katz14 determined the shape of
the trailing wake, by assuming that all vorticity is confined
into the boundary layer on the surface of the body and in a
thin vortex wake behind the body and by using a timestepping procedure. Smith et al.15 also developed an unsteady aerodynamic three-dimensional panel method that can
realistically describe the kinematics of a flapping wing and
can also account for wing flexibility. However, these panel
methods assume an infinite value of the Reynolds number
and, by ignoring the viscosity of the fluid, they fail to predict
the details of the near-field flow, e.g., the large attached
leading-edge vortex that has been recognized to play a key
role in the propulsive efficiency3 and in the lift generation.16
In the present work, the solution of the full NavierStokes equations for the three-dimensional, unsteady flow
generated by a flapping foil of finite span is determined.
Numerical simulations have some advantages with respect to
laboratory experiments. First, via numerical simulations one
has access to velocity and pressure fields in the threedimensional space and in time. Hence, it is possible to make
simultaneous velocity and pressure measurements in different points of the domain, while this would be prohibitively
expensive and would be subject to unacceptable errors due to
probe interference in laboratory experiments. Moreover, using a numerical approach, it is possible to “measure” quantities such as vorticity, strain rate, dissipation, etc., which are
impossible to measure accurately in the laboratory. Even
though computational approaches are limited in their reproduction of the actual flow by the relatively small values of
the Reynolds number that can be efficiently computed, by
the accuracy of the numerical methods !which rapidly deteriorates as the grid size is increased", and by the realism of
the initial and boundary conditions, the use of numerical ap-
˜ , ˜Y , ˜Z" moving with the forward velocity −U
FIG. 1. The reference frame !X
0
of the foil and the reference frame !X , Y , Z" moving with the foil.
proaches is certainly valuable to investigate the details of the
flow structure.
In this paper, our attention is mainly focused on the vortices that are generated by flow separation both at the leading
and trailing edges of the foil and on the vorticity shed at its
tips. As in the paper by Dong et al.,17 a study of the topology
and dynamics of the vortex structures is carried out to provide some physical insight into the mechanisms that control
momentum transfer and the generation of lift and propulsive
thrust. The procedure used in the rest of the paper is as
follows. In the next section, the problem is formulated and
the dimensionless parameters that control the results are
listed. In Sec. III, the numerical code used to perform the
numerical simulations is briefly described. In Sec. IV, the
code is validated by comparing the results with the experimental visualizations of Freymuth.18 Then, the phenomenon
is investigated for values of the parameters close to those
used in the laboratory experiments of von Ellenrieder et al..8
The numerical results are compared against dye visualizations and discussed to point out the main features of the
vortical flow. Some conclusions are drawn in Sec. V, where
both ongoing and future research are also described.
II. THE PROBLEM
A foil, characterized by a finite aspect ratio, moves with
a constant forward speed −U0 and oscillates with a combination of harmonic vertical and angular oscillations !see Fig.
1". The foil has a symmetrical cross section and a rectangular
planform. Its chord and thickness are denoted with C and S,
respectively, while its span is B. If a right-handed reference
˜ , ˜Y , ˜Z" moving with velocity −U is introduced, such
frame !X
0
that the ˜X coordinate is aligned with U0 and the ˜Z coordinate
is in the spanwise direction, the motion of the foil is described as function of the time T by providing the ˜Y coordinate A!T" of the point O around which the foil pitches and
the angle !!T" that the foil chord forms with the ˜X axis,
˜Y = A!T" = Aˆ cos!2"FT",
O
!1"
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-3
Phys. Fluids 17, 113601 "2005!
Numerical experiments on flapping foils
!!T" = ¯! + !ˆ cos!2"FT + #".
!2"
Here Aˆ and !ˆ denote the amplitudes of the heaving and
pitching oscillations, respectively, F is the frequency of the
oscillations, and # is the phase lead of the angular motion
with respect to the lateral motion. Moreover, ¯! denotes the
mean angle of attack.
The flow generated by the fluid flowing around the oscillating foil can be determined by solving momentum
!Navier-Stokes" and continuity equations. By choosing the
foil chord C as length scale, the period 1 / F of the foil oscillations as time scale, and U0 as velocity scale, it can be easily
verified that the dimensionless problem is characterized by
eight dimensionless parameters beside the geometry !thickness and shape" of the foil section: !i" the ratio b = B / C between the span B of the foil and its chord C, !ii" the ratio
a = Aˆ / C between the amplitude of the heaving oscillations
and C, !iii" the amplitude !ˆ of the pitching oscillation, !iv"
the phase # between the pitching and heaving oscillations,
!v" the mean angle of attack ¯!, !vi" the ratio d = D / C, where
D is the distance of the pitching axis from the leading
edge of the foil, !vii" the Reynolds number of the flow
Re= U0C / $, defined with the forward velocity of the foil, the
chord of the foil, and the kinematic viscosity $ of the fluid,
and !viii" the Strouhal number St= 2FAˆ / U0, defined using
the frequency F of the foil oscillations.
The numerical approach solves the problem in the
Cartesian coordinate system !X , Y , Z", which moves with the
foil !see Fig. 1", using primitive variables. The details of the
problem are given in the Appendix.
III. THE NUMERICAL APPROACH
In the present work, the problem has been solved numerically looking for the values of the velocity and pressure,
at a finite number of grid points in the fluid domain, which
approximate the continuous distribution of the actual velocity and pressure fields.
The numerical method uses a finite-difference approach.
Standard centered second-order finite-difference approximations of the spatial derivatives are used. The time advancement of the Navier-Stokes equations, written in the noninertial reference system !X , Y , Z", employs the fractional-step
method extensively described in Refs. 19 and 20. The nonsolenoidal intermediate velocity field is evaluated by means
of a third-order Runge-Kutta scheme to discretize convective
terms, together with a Crank-Nicholson scheme for the diffusive terms. The implicit treatment of the viscous terms
would require the inversion of large sparse matrices, which
are reduced to three tridiagonal matrices by a factorization
procedure with an error of order !%T"3.21 Then, by forcing a
continuity equation, a Poisson equation for the pressure field
is obtained, which is solved by standard methods. The complex geometry that characterizes the foil section is handled
by means of the immersed boundary technique,22,23 which is
based on the introduction of fictitious body forces in the fluid
domain along the contour of the foil. The strength of these
fictitious forces is chosen in such a way that the no-slip con-
dition at the foil surface is forced without the need to introduce any special treatment of the problem at the surface of
the foil.
A nonreflective boundary condition at the outflow
boundary allows the use of a short computational domain in
the X direction. Different computational boxes are used, the
size of which is given in the following. The mesh is uniform
in the spanwise direction, while in the streamwise and vertical ones a nonuniform mesh has been used to cluster the grid
points in the vicinity of the foil where gradients are expected
to be stronger.
Visualizations of the trajectories of passive tracers have
been used to obtain preliminary information on the velocity
field. However, the results are not discussed hereafter because the analysis of the streak lines cannot be used to surely
identify vortex structures. After having gained a qualitative
picture of the flow structure, a quantitative investigation is
made of the coherent vortices, having access to velocity and
vorticity fields in the three-dimensional space and time. In
addition, since the analysis of the isovorticity surfaces may
be considered not fully adequate to detect vortices, we also
use the definitions of vortex structures proposed by Jeong
and Hussain24 and Zhou et al.25 If the unsteady irrotational
straining and the viscous effects are neglected, the pressure
Hessian is proportional to D2 + R2, where D and R are the
symmetric and antisymmetric part of the velocity gradient
tensor, respectively. The above fact can be used to determine
the existence of a local pressure minimum due to vortical
motion. The vortex core is so defined as a connected region
with two negative eigenvalues of D2 + R2.24 Since D2 + R2 is
a symmetric tensor, all eigenvalues are real. If the eigenvalues are ordered in accordance with their decreasing values
!&1 ' &2 ' &3", the new definition is equivalent to the requirement &2 ( 0 within the vortex core. As discussed by
Jeong and Hussain,24 this definition captures the pressure
minimum in a plane perpendicular to the vortex axis at high
Reynolds numbers and also accurately defines vortex cores
at low Reynolds numbers, unlike a pressure minimum criterion, which sometimes is used as an intuitive indicator of a
vortex. Moreover, as suggested by Zhou et al.,25 we visualize
also the imaginary part &ci of the complex eigenvalues of the
velocity gradient tensor. Indeed, if &r is the real eigenvalue
of the velocity gradient tensor, with a corresponding eigenvector Vr, and &cr ± ı&ci are the conjugate pair of the complex eigenvalues, with complex eigenvector Vcr ± ıVci, the
local flow is either stretched or compressed along the axis
Vr, while on the plane spanned by the vectors Vcr and Vci the
flow is swirling. Moreover, the strength of the local swirling
motion is quantified by &ci and hence the imaginary part of
the complex eigenvalue pair is referred to as the local “swirling strength” of a vortex. Numerical experiments suggest the
values of the modulus %!%, the values of the second eigenvalue &2 of the matrix D2 + R2 and of &ci that should be used
to visualize the vortex structures generated by the foil oscillations. While %!% and &ci are of order U0 / C, the second
eigenvalue of D2 + R2 is of order !U0 / C"2. Of course, the
larger the values are, the smaller are the regions visualized
by the approaches. Moreover, it is worth pointing out that the
approaches do not visualize the vortex tubes. Hence, even
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-4
Blondeaux et al.
Phys. Fluids 17, 113601 "2005!
FIG. 3. Square wing in rapid pitch from 0° to 30° in 0.5 s. On the left, the
vortical structures observed by Freymuth !Ref. 18" are sketched on the basis
of his Fig. 3. On the right, the isosurface of %!% = 10U0 / C obtained by means
of the present numerical simulation is plotted. %T = 1 / 8 s.
results it provides with the results obtained with the code
N)*+!r,26,27 developed at Brown University and based on a
finite-element approach. The interested reader is referred to
Blondeaux et al.,28 where a detailed quantitative comparison
is described.
A. Simulation of Freymuth’s experiments
FIG. 2. Rectangular wing in pitch up-down maneuver from 0° to 30° to 0°
within 0.5 s. On the left, the vortical structures observed by Freymuth !Ref.
18" are sketched on the basis of his Fig. 4. On the right, the isosurface
%!% = 12U0 / C obtained by means of the present numerical simulation is plotted. %T = 1 / 8 s.
though sometimes the vortex structures described in the following seem to begin and end in the flow, the Helmholtz
vortex theorems do not break down.
IV. THE RESULTS
The reliability and the accuracy of the numerical code
presently used is supported by a comparison between the
To further support the numerical results, some of the
flow visualizations of Freymuth18 are reproduced. In a first
run, b is set equal to 2 and the foil profile is chosen similar to
a NACA0012 profile. The computational box in the !X , Y , Z"
coordinate system is 5.5 foil chords long, 3 chords wide, and
4 chords high. Moreover, 240, 151, 220 grid points have
been used in the streamwise, spanwise, and cross-stream directions, respectively. The foil is described using 60, 50
, 13 points, respectively. In order to limit the computational
costs, the Reynolds number is fixed equal to 1000, rather
than 5100 as in Ref. 18. For both Re= 1000 and Re= 5100,
the size of the vortex structures shed by the foil is much
larger than the thickness of the boundary layer that develops
along the foil surface and the thickness of the free shear
layers shed by the foil. Hence, no relevant difference is expected between the dynamics of the vortex structures presently computed and that observed by Freymuth.18 Only, the
viscous decay of the vortices, taking place far from the foil,
is envisaged to be different. Since Freymuth18 visualized the
flow generated by an aperiodic motion of the foil, relationships !1" and !2" are modified in accordance with the experi-
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-5
Numerical experiments on flapping foils
Phys. Fluids 17, 113601 "2005!
FIG. 4. Isocontours of the streamwise vorticity -x at . = −" / 2, when the foil is at its middle heave position, for a = 0.5, !ˆ = 5 ° , St= 0.35, # = 90° , Re
˜ = 0. Second row: portion of the plane ˜Z = 1.55C, located at distance equal to C / 20 from the foil tip. The results on the
= 164. First row: portion of the plane X
left have been obtained with 311, 188, 325 grid points, the results on the right have been obtained with 311, 188, 257 grid points. In the figure, the spatial
coordinates are made dimensionless using the foil chord C.
mental data. In particular, the transverse position of the foil
is kept fixed and equal to zero. Starting from ! = 0 at T = 0,
the foil is rotated up to a value of ! = 30° at T = 0.25 s and
then rotated back to 0° at T = 0.5 s, where it is again kept
fixed. Figure 2 shows the numerical predictions of the vorticity field at different times, along with the streak lines visualized by Freymuth.18 A vortex ring is formed by the two
counter-rotating trailing edge vortices that develop during
the pitch-up and pitch-down phases and by the vorticity shed
by the tips of the foil. The remaining parts of the edge vortex
loops, which were previously connected to the leading-edge
corners, annihilate. Also, the leading-edge dynamic stall vortex filaments seem to form a ring. However, in the experiments the dye pattern becomes so entangled that the leadingedge ring has the appearance of a puff. In the numerical
simulation, the leading-edge vortices give rise to a vortex
structure that moves downstream along the lee side of the
foil and then is shed in the free stream. The qualitative agreement between the numerical simulations and the flow visualizations of Freymuth18 is good and supports the reliability
of the numerical results, even though the vortex structures in
the laboratory experiments are slightly farther away from the
foil than the numerically simulated vortices. Some of the
discrepancy might be due to the deficiencies of the numerical
approach and in particular to the difference between the experimental and numerical Reynolds numbers. However, part
of the quantitative discrepancy can be certainly ascribed to
the possible differences between the motion of the foil in the
laboratory experiment and that in the numerical simulation.
Indeed, the details of the foil kinematics are not given in Ref.
18 so that a sinusoidal time development has been presently
chosen for the foil angular velocity.
Figure 3 shows a comparison between the computed and
the visualized vorticity fields for a foil that, starting from a
vanishing value of the angle of attack, is rotated up to 30° in
0.5 s and then kept fixed. In this case, the dimensionless span
of the foil is equal to 1 !b = B / C = 1". The large vortex loop
that develops on the right of the foil is formed by the vorticity generated along the foil and separated from the trailing
edge and from the wing tips during the pitch-up. The loop
connects to the front corners of the wing. Simultaneously,
because of the high angle of attack, vorticity also separates
from the leading edge and a dynamic stall vortex appears.
This connects to the front corners of the wing, thus completing a closed vortex system in accordance with Helmholtz’s
law. Subsequently, the anchoring of the combined leadingand trailing-edge vortex system to the front corners gets
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-6
Blondeaux et al.
Phys. Fluids 17, 113601 "2005!
FIG. 5. Isosurface of %!% = 1.5U0 / C for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35. Side and plan views at !a" . = −" / 2, !b" . = −" / 4, !c" . = 0, !d"
. = " / 4, after about four cycles. In the figure the spatial coordinates are made dimensionless using the foil chord C.
somewhat lost, presumably because of viscous and turbulent
annihilation of the nearly counter-rotating leading- and
trailing-edge tip vortices. It should be mentioned that in the
experimental procedure of Freymuth,18 some smoke was inadvertently introduced in regions where no vorticity exists
and contributed to the “optical” noise of the photograph.
Also in this case the qualitative agreement between the flow
visualizations and the numerical results is good, even though
quantitative differences are present. As previously pointed
out, these differences might be due to deficiencies of the
numerical simulation and to the possible differences between
the kinematics of the foil in the laboratory and numerical
experiments. Indeed, after 1 / 8 s, the trailing-edge vortex is
well developed and it has moved far from the foil in the
laboratory experiment, while it is weaker and still close to
the foil in the numerical simulation, thus indicating that the
actual rotation of the foil is different.
B. Simulation of von Ellenrieder et al.’s experiments
As pointed out in the formulation of the problem, the
results of the numerical code depend on eight dimensionless
parameters beside the geometry !thickness and shape" of the
foil. Hence, the high computational costs do not allow an
exhaustive investigation of the phenomenon in the parameter
space. In this work, the phenomenon has been simulated for
values of the parameters equal to those considered in the
experiments carried out by von Ellenrieder et al.8 A detailed
investigation of the velocity and vorticity fields has been
carried out to identify possible differences between the
streak-line behavior described by von Ellenrieder et al.8
and the actual vorticity field and to support or modify the
isometric sketch of the flow structure proposed by the above
authors.
A foil characterized by a value of b equal to 3 has been
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-7
Phys. Fluids 17, 113601 "2005!
Numerical experiments on flapping foils
FIG. 6. Isosurface of %!% = 2.0U0 / C for the region
specified in Fig. 5!c". Side and plan views. In the figure,
the spatial coordinates are made dimensionless using
the foil chord C.
considered and the section of the foil has been taken similar
to that of a NACA0030 profile. The dimensionless distance d
of the pitching axis from the leading edge of the foil is equal
to 1 / 4 and the ratio a between the amplitude of the heaving
oscillations and the foil chord is fixed equal to 0.5. Moreover, the Reynolds number is about 164, as the experimental
value. The pitching amplitude !ˆ and the phase # between the
two oscillatory motions have been fixed equal to 5° and 90°,
respectively, while a vanishing value of ¯! is considered.
Simulations have been performed for two different values of
the Strouhal number, namely St= 0.35 and St= 0.175. The
computational box in the !X , Y , Z" coordinate system is 11.75
foil chords long, 8 chords high, and 5.6 chords wide. Moreover, 311, 188, 325 grid points have been used in the
streamwise, cross-stream, and spanwise directions, respectively. The foil is described using 118, 36, 177 points, respectively. It takes about 24 h to simulate one oscillation
cycle using a PC AMD Athlon XP 3000. Some preliminary
tests have been made to ascertain that the size of the computational grid is small enough to assure a fair accuracy of
the results. For example, Fig. 4 shows, at a certain phase of
the cycle, a comparison between the distribution of the
streamwise vorticity component obtained with the number of
grid points previously given and that obtained decreasing the
number of grid points in the spanwise direction !257 grid
points have been used". The plots in the first row show a
portion of the plane ˜X = 0, while the plots in the second row
show a portion of the plane ˜Z = 1.55C, which is located at a
distance equal to C / 20 from the foil tip. The agreement between the two vorticity fields supports the accuracy of the
obtained results. The streamwise vorticity component and the
region around the foil tips are shown because strong velocity
gradients are expected along the spanwise direction at the
tips of the foil and the code employs a uniform grid spacing
in the Z direction. A run has been also made with a shorter
computational box !7.5, 8 , 5.6 foil chords" and using the
same number of points, so that it has been verified that the
flow close to the foil and in the far wake is not affected by
the length of the computational box. Moreover, at the beginning of the calculations, both the fluid and the foil, which are
initially at rest, are impulsively set into motion. All the simulations have been carried out for a number of cycles large
enough to assure that the flow close to the foil no longer
feels the influence of the initial conditions.
The time development of the flow over half a heave/
pitch cycle is shown for St= 0.35 in Fig. 5, where an isosurface of the modulus of vorticity !%!% = 1.5U0 / C" is plotted at
different phases of the cycle, after four cycles of oscillation.
As in von Ellenrieder et al.,8 the figure displays both the plan
and side views in order to give a sense of the threedimensional structure of the flow. Motion parameters and
cycle phases are the same of von Ellenrieder et al.’s8 Figs.
3!c"–3!f". By looking at Fig. 6, which is an enlargement of
Fig. 5!c", it is possible to recognize the following:
!i"
!ii"
!iii"
!iv"
!v"
The leading-edge vortex !side view" moving along the
lower surface of the foil.
The trailing-edge vortex rolling close to the foil in the
span-wise direction !plan view".
The leading-edge vortex, shed during the previous
half-cycle, which merged with trailing-edge vorticity
and constitutes the head of the developing vortex ring
!side view".
The tip vorticity !plan view".
The vortex rings shed at the previous half-cycles.
These findings are similar to those of Ref. 8 even though the
ring vortex structures identified in the laboratory experiments
are slightly longer and larger. Since the locations of dye
emission, which strongly influence streak-line patterns, are
not specified in Ref. 8 and taking into account that it is quite
difficult to reproduce the details of dye injection !geometry
of the injectors, initial velocity of dye, etc.", no attempt has
been made to make a quantitative comparison between the
computed streak lines and the visualizations of von Ellenrieder et al.8
Since the key aspect of the work is the identification of
vortical structures, we have tried some of the various techniques proposed in the literature. A first, simple, and quite
direct method to identify vortex structures is represented by
the isosurfaces of the strength of vorticity !see Fig. 5", even
if the largest values of the vorticity occur at the foil surface
or close to it, where the flow may exhibit shear but not necessary a swirling motion.24 Moreover, in the experiments,
dissipative effects are rather strong and the vorticity decays
along the wake quite rapidly, so that small values of %!%
should be used if the isovorticity surface should identify the
vortex rings in their completeness. Hence, we cannot avoid
losing the details of the vortex structures in those regions
where vorticity is large.
Another simple and quite direct method to identify vortex structures is represented by the isosurfaces of the pressure distribution !see Blondeaux et al.28". Indeed, local pressure minima are related to the presence of vortex structures.
However, there is an inherent scale difference between a vor-
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-8
Blondeaux et al.
Phys. Fluids 17, 113601 "2005!
FIG. 7. Isosurface of &ci2 #&ci2 = 0.4!U0 / C"2$ for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35, in a perspective view that looks at the foil from the
bottom. !a" . = −" / 2, the foil is drawn in black gray; !b" . = −" / 4, the isosurface is transparent and superimposed to &ci2 contour levels in the plane ˜Z = 0; !c"
. = 0, the isosurface is transparent and superimposed to &ci2 contour levels in the planes ˜Z = 0 and ˜Y = 0.2C; !d" . = " / 4, the isosurface is transparent and
superimposed to /ci2 contour levels in the planes ˜Z = 0 and ˜Y = 0. In !b", !c", and !d", the foil is not displayed.
tex core and the associated low-pressure region, and this
makes the demarcation of vortices using an isosurface of
pressure problematic.24 The so-called “&ci” criterion seems
to be more appropriate to the description of the vortex dynamics both in the region close to the foil and along the
wake. In Fig. 7, a surface characterized by a constant value
2
is shown at four phases of the cycle. Now the structure
of &ci
of the vorticity field close to the foil becomes more understandable, as does the structure of the vortex rings that leave
the foil every half a cycle.
Looking at dye filaments, von Ellenrieder et al.8 distinguished between the vortex structures shed by the leading
edge of the foil, which they labeled Li !the index i indicating
the half-cycle", and the vortices shed by the trailing edge,
which they labeled Ti. However, if the vorticity time development is analyzed in the plane of symmetry of the foil !see
Fig. 8", the present results show that the vorticity shed by the
leading edge merges with that generated all along the foil
surface and that shed by the trailing edge, to produce a
unique vortex structure every half-cycle. An analysis of the
sequential perspective views of the regions characterized by
a constant value of &2ci shows that these unique vortex structures are shed in alternating order near to the extrema of the
heaving motion. These vortices tend to be concave and form
vortex rings being connected to the tip vorticity and the
spanwise vorticity of opposite sign shed during the following
half-cycle. When the vortex rings are convected downstream,
the vortex structures are stretched and tend to assume a hook
shape in the vertical plane. As the cycle continues, the vortex
filaments move toward the midspan plane and simultaneously wrap around and combine with the head of the following vortex. Comparing the results displayed in Figs. 5
and 7 with those of Fig. 3 of von Ellenrieder et al.,8 it can be
˜ , ˜Z"
seen that the vortex wake of the foil is smaller in the !X
˜ , ˜Y " plane.
plane and larger in the !X
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-9
Numerical experiments on flapping foils
Phys. Fluids 17, 113601 "2005!
˜ = 0" for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35, at !a" . = −" / 2, !b"
FIG. 8. Contour levels of -z in the plane of symmetry of the foil !Z
. = −" / 4, !c" . = 0, !d" . = " / 4 , %-z = 1.85U0 / C. In the figure, the spatial coordinates are made dimensionless using the foil chord C.
If we look at the contour plot of the spanwise vorticity
-z in the plane of symmetry of the foil ˜Z = 0 !see Fig. 8", we
find that the vorticity generation along the surface of the foil
is not much different from that characterizing the twodimensional flow around a foil of infinite length !see Fig. 9".
Plots refer to four phases during half-cycle. During the upward motion of the foil, flow separation from the leading
edge is observed, associated with a weak back flow along the
lower surface of the foil. The free shear layer rolls up and
gives rise to a negative vortex structure that is convected
downstream. When the foil reverses its motion, new negative
vorticity is generated around the trailing edge of the foil
because of its downward motion. The evolution of these two
vortices can be followed looking at the plot for . = −" / 2
#Fig. 8!a"$ and considering the mirror image of the positive
vorticity shed during the previous half-cycle at the upper
FIG. 9. Contour levels of -z for
the 2D simulation of a high aspect ratio foil similar to a NACA0030, with
a = 0.5, !ˆ = 5°, ¯! = 0, # = 90°, Re= 164,
St= 0.35, at !a" . = −" / 2, !b"
. = −" / 4, !c" . = 0, !d" . = " / 4 , %-z
& 0.46U0 / C. In the figure, the spatial
coordinates are made dimensionless
using the foil chord C.
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-10
Phys. Fluids 17, 113601 "2005!
Blondeaux et al.
FIG. 11. Isosurface of &2 = −0.7!U0 / C"2 for !ˆ = 5 ° , a = 0.5, # = 90° , St
= 0.175 at . = 3" / 4.
FIG. 10. Isosurface of &2#&2 = −0.3!U0 / C"2$ for !ˆ = 5 ° , a = 0.5, #
= 90° , St= 0.35, at . = 3" / 4 after about four cycles. The isosurface is translucent and superimposed to contour levels of the streamwise vorticity -x at
planes ˜X = −2 . 0C , − 3.3C , − 4.65C !negative values are dashed". Arrows
indicate the approximate direction and versus of vorticity.
surface of the foil. Far from the foil, some differences appear
between the 2D and 3D cases. Indeed, in the 2D case, the
vorticity strength decreases less rapidly as the vortex structures are convected downstream. Moreover, looking at the
isosurface of &ci at the same phases of the cycle in the threedimensional case !see Fig. 7, where the foil is seen in a
bottom view", the same vortices with positive -z generated
at the upper surface of the foil, namely the leading-edge
vortex, which is going to constitute the head of the new
vortex ring, and the trailing-edge vortex, can be recognized.
However, these vortices are connected by the vorticity generated at the tips of the foil: the two arms of the tip vorticity
elongate, lose connection with the tips, and a ring is formed
and shed when the foil is leaving its extreme transversal
position. The leading-edge vortex with negative -z shed at
the lower surface of the foil contributes to the formation of
the ring connecting the tip vorticity at the front.
Looking at the spanwise velocity component along a different vertical plane crossing the foil at different values of
the coordinate ˜X, we can observe a spanwise symmetrical
transport of momentum toward the foil extrema to feed the
strong tip vortices. These results agree with Maxworthy’s29
finding, which shows a strong spanwise transport of vorticity
toward the wing tips of a three-dimensional foil.
The arrows in Fig. 10, lying along the vortex structures
identified by a constant negative value of &2, indicate the
direction of circulation and have been obtained plotting the
different vorticity components provided by the calculation.
Indeed, even though the &ci and &2 criteria are reliable approaches to visualize the regions where vortex structures are
present, they are unable to provide the sense of rotation of
the vortex structures. In accordance with the rule that no
vortex can simply begin or end in the flow, the vortices form
a continuous system of interconnected loops that extends
downstream of the foil.
To give another example, further results are shown in
Fig. 11, where an isosurface of &2 is represented for
a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.175, i.e., for
a smaller value of the Strouhal number. In this case, the shed
vorticity is weaker, the vortex filaments turn out to be less
curved, and the characteristic length of the vortex structures
is approximately twice that for St= 0.35. Moreover, the wake
tendency to become narrower in the ˜Z direction and to open
and diverge in the ˜Y direction is less evident.
V. CONCLUSIONS
In Lighthill,30 it is first explained how vortex lines shed
from the central part of the tail/foil must bend downstream,
then turn inward and close up when they come close to the
vorticity of opposite sign shed by the tail/foil half a period
later. In accordance with Lighthill,30 if the foil is supposed to
flap in the vertical plane moving alternatively up and down,
the wake turns out to consist of a series of horizontal vortex
rings.
The experimental works of Drunker and Lauder31,32
and Lauder,33 as well as the recent flow visualizations of
von Ellenrieder et al.,8 have confirmed the formation of ringlike vortical structures by oscillating tails/foils.
The present numerical simulations also show the generation of vortex rings by an oscillating foil. Indeed, the results
show that, for the chosen parameters, the leading vortices,
shed by the foil for high angle of attack, merge with the
trailing-edge vorticity, giving rise to large vortex structures.
These originate ring-like structures being connected to alternating sign tip vortices and to the spanwise vorticity of opposite sign shed during the following half-cycle.
Numerical simulations show that sometimes the rings
have a complex geometry and are characterized by a complex dynamics. In particular, on the basis of the present results, it is possible to guess the existence of reconnection
phenomena at relatively high Strouhal numbers. As described in Kida and Takaoka,34 first the vorticity tubes tend
to approach each other advected by the mutual- and selfinduction velocity. As the two vortex tubes get closer, the
shape of the vortex core is deformed.35 This can be seen in
Fig. 12, which shows a detail of the vorticity field plotted in
Fig. 7. Then, the generation of vorticity in the direction per-
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-11
Phys. Fluids 17, 113601 "2005!
Numerical experiments on flapping foils
FIG.
12.
Isosurface
of
&2 #&2 = −0.3!U0 / C"2$ for !ˆ = 5 ° , a
= 0.5, # = 90° , St= 0.35 at . = 3" / 4
after about four cycles. The region that
is emphasized in the top figure is enlarged in the two bottom plots. In the
bottom left figure, the isosurface of &2
is colored with the contour levels of
-x. In the right bottom figure, the isosurface of &2 is colored with the contour levels of -z.
' (
pendicular to the original vortex tubes, which can be observed in Fig. 12, gives rise to vorticity reconnection. However, the viscous cancellation of opposite signed vorticity is
not complete in Fig. 12 and the original rings are not definitively broken.
Recent detailed data on three-dimensional foils by
Martin et al.36 show little degradation of propulsive performances for moderate aspect ratio foils, compared with the
results of two-dimensional foils. The evaluation of the forces
acting on the foil should clarify this point and the simulation
of different foil geometries should clarify the apparent evolutionary convergence of the fish tail toward a lunate form
and/or suggest another form to obtain optimal efficiency.
# V d 2!
d!
− 2X+
# T dt
dt
+V
2
Y−
#V
d 2A
d!
U+U
2 cos ! − 2
#X
dt
dt
)
*
#V
#V
1 #P
# 2V # 2V # 2V
+W
=−
+$
+
+
,
#Y
#Z
0 #Y
# X2 # Y 2 # Z2
!A3"
#W
#W
#W
#W
+V
+W
+U
#X
#Y
#Z
#T
=−
)
*
1 #P
# 2W # 2W # 2W
+$
+
+
.
0 #Z
# X2 # Y 2 # Z2
!A4"
ACKNOWLEDGMENTS
Thanks are due to the Office of Naval Research for funding the research under Contract No. 000140310193 monitored by Dr. Thomas F. Swean.
APPENDIX
Continuity and momentum equations in the moving reference frame !X , Y , Z":
#U #V #W
+
+
= 0,
#X #Y #Z
' (
# U d 2!
d!
+
Y−
# T dt2
dt
+V
!A1"
2
X−
#U
d 2A
d!
sin ! + 2 V + U
#X
dt2
dt
)
*
#U
#U
1 #P
# 2U # 2U # 2U
+W
=−
+$
+
+
,
#Y
#Z
0 #X
# X2 # Y 2 # Z2
!A2"
1
M. S. Triantafyllou, A. H. Techet, and F. S. Hover, “Review of experimental work in biomimetic foils,” IEEE J. Ocean. Eng. 29, 585 !2004".
2
R. Mittal, “Computational modeling in biohydrodynamics: Trends, challenges, and recent advances,” IEEE J. Ocean. Eng. 29, 595 !2004".
3
J. M. Anderson, K. Streitlien, D. S. Barret, and M. S. Triantafyllou,
“Oscillating foils of high propulsive efficiency,” J. Fluid Mech. 360, 41
!1998".
4
Z. J. Wang, “Vortex shedding and frequency selection in flapping flight,” J.
Fluid Mech. 410, 323 !2000".
5
G. C. Lewin and H. Haj-Hariri, “Modelling thrust generation of a twodimensional heaving airfoil in a viscous flow,” J. Fluid Mech. 492, 339
!2003".
6
G. Pedro, A. Suleman, and N. Djilali, “A numerical study of the propulsive efficiency of a flapping hydrofoil,” Int. J. Numer. Methods Fluids 42,
493 !2003".
7
L. Guglielmini and P. Blondeaux, “Propulsive efficiency of oscillating
foils,” Eur. J. Mech. B/Fluids 23, 255 !2004".
8
K. D. von Ellenrieder, K. Parker, and J. Soria, “Flow structures behind a
heaving and pitching finite-span wing,” J. Fluid Mech. 490, 129 !2003".
9
D. I. Pullin and A. E. Perry, “Some flow visualization experiments of the
starting vortex,” J. Fluid Mech. 97, 239 !1980".
10
F. R. Hama, “Streaklines in a perturbed shear flow,” Phys. Fluids 5, 644
!1962".
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
113601-12
11
Blondeaux et al.
K. Parker, K. Ellenrieder, and J. Soria, “Using stereo multigrid DPIV
!SMDPIV" measurements to investigate the vortical skeleton behind a
finite-span flapping wing,” Exp. Fluids 39, 281 !2005".
12
K. Parker, K. von Ellenrieder, and J. Soria, “The flow of a threedimensional thrust producing wing,” in Advances in Turbulence, 10th European Turbulence Conference, Fluid Mechanics and Its Applications, edited by P. A. Krogstad and H. I. Andersson !CIMNE, Trondheim, Norway,
2004", Vol. 49, pp. 11–14.
13
J. Katz and A. Plotkin, Low-Speed Aerodynamics from Wing Theory to
Panel Methods !McGraw-Hill, New York, 1991".
14
M. S. Vest and J. Katz, “Unsteady aerodynamic model of flapping wing,”
AIAA J. 34, 1435 !1996".
15
M. Smith, P. Wilkin, and M. Williams, “The advantages of an unsteady
panel method in modelling the aerodynamic forces on rigid flapping
wings,” J. Exp. Biol. 199, 1073 !1996".
16
C. P. Ellington, C. van der Berg, A. P. Willmott, and A. L. R. Thomas,
“Leading-edge vortices in insect flight,” Nature 384, 626 !1996".
17
H. Dong, M. Mittal, M. Bozkurttas, and F. Najjar, “Wake structure and
performance of finite aspect ratio flapping foils,” AIAA Technical Papers
2005-81, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV
!2005".
18
P. Freymuth, “Visualizing the connectivity of vortex systems for pitching
wings,” J. Fluids Eng. 111, 217 !1989".
19
J. Kim and P. Moin, “Application of fractional step method to incompressible Navier-Stokes equations,” J. Comput. Phys. 59, 308 !1985".
20
M. M. Rai and P. Moin, “Direct simulations of turbulent flow using finitedifference schemes,” J. Comput. Phys. 96, 15 !1991".
21
R. M. Beam and R. F. Warming, “An implicit finite-difference algorithm
for hyperbolic systems in conservation law form,” J. Comput. Phys. 22,
87 !1976".
22
J. Mohd-Yusof, “Combined immersed boundaries/B-splines methods for
simulations of flows in complex geometries,” Technical Report, CTR Annual Research Briefs, NASA Ames, Stanford University, 1997.
23
E. A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof, “Combined
Phys. Fluids 17, 113601 "2005!
immersed-boundary finite-difference methods for three-dimensional complex flow simulations,” J. Comput. Phys. 161, 35 !2000".
24
J. Jeong and F. Hussain, “On the identification of a vortex,” J. Fluid Mech.
285, 69 !1995".
25
J. Zhou, R. J. Adrian, S. Balachandar, and T. M. Kendall, “Mechanisms
for generating coherent packets of hairpin vortices in channel flow,” J.
Fluid Mech. 387, 353 !1999".
26
S. J. Sherwin, Ph.D. thesis, Princeton University !1995".
27
T. C. Warburton, Ph.D. thesis, Brown University, Division of Applied
Mathematics !1999".
28
P. Blondeaux, F. Fornarelli, L. Guglielmini, M. S. Triantafyllou, and R.
Verzicco, “Vortex structures generated by a finite-span oscillating foil,”
AIAA Technical Papers 2005-84, 43rd AIAA Aerospace Sciences Meeting
and Exhibit, Reno, NV !2005".
29
T. Maxworthy, “Experiments on the Weis-Fogh mechanism of lift generation by insects in hovering flight. Part 1. Dynamics of the ‘fling’,” J. Fluid
Mech. 93, 47 !1979".
30
M. J. Lighthill, “Hydromechanics of aquatic animal propulsion,” Annu.
Rev. Fluid Mech. 1, 413 !1969".
31
E. G. Drucker and G. V. Lauder, “Locomotor forces on a swimming fish:
Three-dimensional vortex wake dynamics quantified using digital particle
image velocimetry,” J. Exp. Biol. 202, 2393 !1999".
32
E. G. Drucker and G. V. Lauder, “A hydrodynamic analysis of fish swimming speed: Wake structure and locomotor force in slow and fast labriform swimmers,” J. Exp. Biol. 203, 2379 !2000".
33
G. V. Lauder, “Function of the caudal fin during locomotion in fishes:
Kinematics, flow visualization, and evolutionary patterns,” Am. Zool. 40,
101 !2000".
34
S. Kida and M. Takaoka, “Vortex reconnection,” Annu. Rev. Fluid Mech.
26, 169 !1994".
35
S. Kida, M. Takaoka, and F. Hussain, “Formation of head tail structure in
a 2-dimensional uniform straining flow,” Phys. Fluids A 3, 2688 !1991".
36
C. B. Martin, F. S. Hover, and M. S. Triantafyllou, “Propulsive performances of a rolling and pitching wing,” Proceedings of the UUST, International Symposium on Unmanned Untethered Submersible Technology,
Durham, NH !2001".
Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp