|
| |
|
|
|
|
||||
| Home Help Feedback Subscriptions Archive Search Table of Contents | ||||
First published online December 1, 2006
Journal of Experimental Biology 209, 4841-4857 (2006)
Published by The Company of Biologists 2006
doi: 10.1242/jeb.02526
Simulations of optimized anguilliform swimming
Institute of Computational Science, ETH Zurich, CH-8092, Switzerland
* Author for correspondence (e-mail: petros{at}inf.ethz.ch)
Accepted 10 October 2006
 |
Summary |
|---|
 TOP Summary IntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
The kinematics of burst swimming is characterized by the large amplitude of the tail undulations while the anterior part of the body remains straight. In contrast, during efficient swimming behavior significant lateral undulation occurs along the entire length of the body. In turn, during burst swimming, the majority of the thrust is generated at the tail, whereas in the efficient swimming mode, in addition to the tail, the middle of the body contributes significantly to the thrust. The burst swimming velocity is 42% higher and the propulsive efficiency is 15% lower than the respective values during efficient swimming.
The wake, for both swimming modes, consists largely of a double row of vortex rings with an axis aligned with the swimming direction. The vortex rings are responsible for producing lateral jets of fluid, which has been documented in prior experimental studies. We note that the primary wake vortices are qualitatively similar in both swimming modes except that the wake vortex rings are stronger and relatively more elongated in the fast swimming mode.
The present results provide quantitative information of three-dimensional fluid-body interactions that may complement related experimental studies. In addition they enable a detailed quantitative analysis, which may be difficult to obtain experimentally, of the different swimming modes linking the kinematics of the motion with the forces acting on the self-propelled body. Finally, the optimization procedure helps to identify, in a systematic fashion, links between swimming motion and biological function.
Key words: anguilliform swimming, hydrodynamics, fluid-body interaction, self-propelled, optimization, wake pattern
|
Introduction |
|---|
|
TOPSummaryIntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
). Starting from the pioneering work of Gray in 1933
(Gray, 1933), anguilliform
swimming has attracted the attention of researchers from diverse scientific
fields, ranging from neuroscience to hydrodynamics
(Ekeberg, 1993;
Ijspeert and Kodjabachia,
1999; Tytell and Lauder,
2004).
A number of intriguing aspects regarding this mode of locomotion remain
unknown, including the link between kinematics and the forces that propel the
body, as well as the overall efficiency of anguilliform swimming
(Schmidt, 1923;
van Ginneken et al., 2005).
The quantitative analysis of the flow characteristics are instrumental in
elucidating fundamental swimming mechanisms. A number of recent studies have
presented flow field measurements of living eels using Particle Image
Velocimetry (PIV). PIV has been employed in order to visualize the flow field
around freely swimming juvenile eels
(Müller et al., 2001) and
the authors postulate that the wake structure is determined by a phase lag
between the vorticity shed from the tail and vortices produced along the body.
In addition, the authors suggest that eels can modify their body wave to
achieve different propulsive modes. Recent experiments employed high
resolution PIV to quantify the wake structure and to measure the swimming
efficiency of Anguilla rostrata
(Tytell, 2004;
Tytell and Lauder, 2004). The
results indicate the shedding of two vortices for each half tail beat
resulting in a wake with strong lateral jets and the notable absence of
substantial downstream flow.
The goal of this study was to investigate anguilliform swimming using
three-dimensional simulations of a self-propelled eel-like body immersed in a
viscous fluid. In addition, two-dimensional simulations are compared with
existing related works (Carling et al.,
1998), in order to assess the role of two- and three-dimensional
simulations in quantifying the fluid dynamics of anguilliform swimming.
The present simulations are coupled with an optimization procedure in order
to test the hypothesis that distinct motion kinematics and resulting wake
patterns correspond to different swimming modes
(Müller et al., 2001). In
this work, the fish motion is not a priori specified. Instead, the
motion parameters are obtained by optimizing objective functions corresponding
to fast and efficient swimming, using an evolutionary optimization algorithm.
We employ the `Covariance Matrix Adaptation Evolution Strategy' (CMA-ES)
(Hansen and Kern, 2004;
Hansen and Ostermeier, 2001).
The CMA-ES encodes relations between the parameters of the body motion and the
objective that is being optimized in an efficient manner, without requiring
explicitly the gradients of the objective function.
The simulations provide detailed information of the complete flow field
enabling the quantification of the vortex formation and shedding process and
allow comparisons with related experimental works
(Müller et al., 2001;
Tytell, 2004;
Tytell and Lauder, 2004). In
particular, the present simulations help to clarify discrepancies between
previous two-dimensional computational studies
(Carling et al., 1998) and
experimental work (Müller et al.,
2001; Tytell,
2004; Tytell and Lauder,
2004).
In addition, the simulations enable the identification of the force distribution along the self-propelled body and their link with the kinematics of the body and the vorticity dynamics of the wake. It is demonstrated that distinct swimming kinematics and body force distribution correspond to different swimming modes. At the same time, differences in terms of vortex wake dynamics are more subtle. The primary vortex rings in the fast swimming mode have a higher strength and are relatively more elongated along the swimming direction.
|
Materials and methods |
|---|
|
TOPSummaryIntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
Geometrical model
The three-dimensional geometry of the anguilliform swimmer is constructed
from spatially varying ellipsoid cross sections. The length of the two half
axes w(s) and h(s) are defined as
analytical functions of the arc length s along the mid-line of the
body. Following the work of Carling et al.
(Carling et al., 1998), the
width of the body is set to its length L at the head and at the tail.
The analytical description of is divided into three regions:
 | (1) |
where wh=sb=0.04L,
st=0.95L and wt=0.01L.
The height h(s) is set as an elliptical curve with the two
half axis a=0.51L and b=0.08L:
 | (2) |
Fig. 1 presents the form of
the ellipsoid cross sections and the geometry of the resulting
three-dimensional (3D) body. The body length L is normalized to 1,
resulting in a surface area of S=0.304 for the fish model. In order
to compute the mass m and inertial moment Iz, the
body is discretized into ellipsoid disks of constant density along the
mid-line s with a step size of 
s=0.001L. The
deformation of the body is composed from the motions of the individual
segments.
|
) proposed a generic model by fitting data available
from observations of undulating fish. In this work we propose a
parameterization defined by the instantaneous curvature along the mid-line of
the body. The choice of curvature as a controlling parameter was dictated by
the fact that curvature relates to bending moments and thus to muscle
contractions along the swimming body. This parameterization allows for direct
control over unrealistic body deformations that may emerge during the
optimization procedure.
We define the curvature of the mid-line 
s as:
 | (3) |
where K(s) is the cubic spline through the m
interpolation points Ki, i=1,...,m,
(s) is the phase shift along the body, t is the time
variable and T is the cycle time. In the current study the phase
shift (s) is linearly proportional to the position s
along the body:
 | (4) |
The parameters Ki and tail are
determined by the optimization procedure. The description of the mid-line is
transformed to the 3D Cartesian inertial coordinate system as described in the
Appendix (A).
In order to compare our simulations with the results of Carling et al.
(Carling et al., 1998), we
transformed their motion parameterization for a body of length L=1
to:
 | (5) |
where ys describes the lateral displacement of the mid-line in a local coordinate system.
Optimization of the motion
The motion parameters of the body are obtained through an optimization
process. Two objectives have been maximized, namely the fish speed (burst
swimming velocity) and the swimming efficiency of the self-propelled body.
These objectives may correspond to biological functions of the swimmer, such
as hunting (for the burst velocity) or migrating (for the efficient swimming).
The parameters of the optimization, Ki and
tail, define in turn a realization of a motion pattern
according to Eqn 3.
Optimal parameters were obtained using a stochastic optimization technique
to obviate difficulties related to noise, discontinuities, and multiple optima
in fitness functions usually associated with fluid mechanics problems
(Milano and Koumoutsakos,
2002). We implement an `Evolution Strategy with Adaptation of the
Covariance Matrix' (CMA-ES) (Hansen and
Kern, 2004; Hansen and
Ostermeier, 2001). The efficiency of CMA-ES has been demonstrated
in a number of benchmark optimization problems
(Kern et al., 2004) and is
essential in the present study, as the evaluation of a single motion pattern
is computationally intensive (about 3 h computation time on a workstation with
a 3 GHz Pentium IV processor).
We have tried to maintain a minimum number of parameters and found that
m=4 interpolation points, evenly distributed at
s=0,
,
,1 with a linear phase shift (s)
are sufficient to allow a wide range of motion patterns. The choice of these
five parameters Ki, i=1,...,4 and
tail was verified from the capability to provide a close
approximation of the motion pattern specified in Carling et al.
(Carling et al., 1998).
We note here that the definition of efficiency for steady undulatory
swimming is controversial (Schultz and
Webb, 2002). The difficulty arises because the regions of drag and
thrust production are distributed over the whole body and their locations vary
during an undulation cycle. Hence drag and thrust components could only be
distinguished by summing positive and negative contributions to the net force
separately in every instant in time.
In this paper we relate optimizing swimming efficiency 
to maximizing
the amount of mean total input power
total that is being
transformed into forward motion:
 | (6) |
In the present optimization, the period T is constant and we
integrate total over the
undulation cycle in order to obtain the amount of work
Wcycle used by the swimmer. This work is computed as a
time integral of the power output of the surface S of the deforming
body on the surrounding fluid:
 | (7) |
where
is the viscous stress tensor, n the outer surface normal vector, and
u the velocity of the moving surface. We postulate that the time
integral of forward motion
is related to the kinetic energy of the forward motion of the body
Ekin=
m
||2.
Hence, we characterize the swimming efficiency by the ratio
Ekin/Wcycle, and we define
the objective function for optimizing swimming efficiency as:
 | (8) |
The objective function for burst swimming is defined as the sum of the mean
forward velocity || and two terms that penalize
high mean and peak input power requirements:
 | (9) |
where 
(.) is the Heaviside function. The penalization of high input
power is motivated by natural physiological limits, and was deemed to be
necessary after obtaining rather unnatural motion patterns (though resulting
in a fast swimming velocity) using only || as
the objective function. The constraints for the mean and peak input power
where chosen as
max=2.0x10-3
and Pmax=4.8x10-3 based on preliminary
computational experiments. The scaling factors where set to
c1=106 and
c2=108. The maximum curvature of the mid-line
was constrained in both optimization cases, in the interval [0,2
] (note
that for a curvature of s
2 the body would bend in
a full circle).
Equations and numerical method
The system of the deforming body interacting with the surrounding fluid is
described by the 3D incompressible Navier-Stokes equations:
 | (10) |
 | (11) |
where 
, and
g are the body forces per unit mass. The feedback from the body to the
fluid is realized by imposing no slip boundary conditions on the moving fish
surface relating the fluid velocity at the boundary uf,
and the body velocity uc:
 | (12) |
The motion of the self-propelled body is in turn described by Newton's
equations of motion:
 | (13) |
 | (14) |
F and Mz are the fluid force and yaw torque
acting on the body surface, xc is the position of the center
of mass of the body,
its global
angular velocity, m the total mass, and Iz the
(time dependent) inertial moment about the yaw axis. The feedback of the fluid
torque is limited to the yaw direction to simplify computations. The fluid
force F and the torque Mz are computed as follows:
 | (15) |
 | (16) |
The far field boundary condition is set to a constant static pressure,
modeling the fish propelling itself through an infinite tank of still fluid.
The 3D Navier-Stokes equations were solved using the commercial package
STAR-CD v. 3.15. STAR-CD computes flows on arbitrary Lagrangian-Eulerian grids
by solving an additional space conservation equation. Eqn 10 and Eqn 11 are
discretized using a finite volume approach with first-order discretization in
time and second-order in space. The solution of Newton's equations of motion
and the resulting grid deformation and motion are implemented in user
subroutines linked to STAR-CD. This setup allows for explicit coupling of
fluid and body physics only. The implemented coupling procedure is a staggered
integration algorithm (Farhat and Lesoinne,
2000); a detailed description is given in the Appendix (B).
The fluid flow is computed on a single block structured grid. The
computational domain around the deforming body is a cylinder with a
hemispherical end at the head of the body. The outer domain boundary is
translated according to the mean kinematics of the swimming body while the
deformation of the surface of the body is propagated to the interior grid
cells using a grid deformation approach proposed elsewhere
(Tsai et al., 2001). The grid
is subdivided into smaller blocks, for which a spring network is used to
determine the deformation of the corner points. Subsequently, in each block an
arc-length-based transfinite interpolation is used to propagate the
deformations of the corner points (and the prescribed surface deformations at
the fluid-body boundary) to the interior grid cells.
Fig. 1B illustrates the
deforming grid. The simulation setup and the flow-structure coupling were
validated on the testcase of a sphere falling in liquid due to gravity. The
results of the validation are outlined in the Appendix (C).
All simulations are conducted with a constant viscosity of
µ=1.4x10-4, body length L=1, density
fluid=body==1, and undulation frequency
f=1. The resulting Reynolds number
Re=L||/µ is in the range
2400-3900, depending on the swimming mode, and is comparable to previous
simulations (Carling et al.,
1998). All simulations are started with the body at rest, and the
motion is initialized by gradually increasing the curvature amplitude from
K(s)0 to its designated value during the first cycle.
The optimization cases are run for 12 cycles and the objective function values
are computed using time-averaged values during the last cycle.
The convergence of the results in terms of spatial and temporal resolution
was tested using the motion pattern proposed in Eqn 5 (see Appendix). A grid
with 3x105 cells and a time step
t=5x10-4 (corresponding to 2000 time steps
per undulation cycle) was shown to accurately resolve the flow. During the
optimization, a grid with 7x104 cells and
t=5x10-3 was used in order to reduce the
computational cost while still allowing a satisfactory computation of the
fluid-body interaction and the resulting dynamics. The results of the
optimization were verified by simulations of high resolution.
|
||2S)
and C 
=F/(||2S)
parallel and lateral to the mean swimming direction. The yaw torque exerted by
the fluid and the total input power are measured in the nondimensional
coefficients
CM=Mz/(||2LS)
and
CP=Ptotal/(||3S).
For the 2D case, the forces are computed per unit length and are normalized by
the circumference of the body. The Strouhal number of the body is computed as
St=2fA/U||, where A is the
amplitude of the undulating tail. The local body wave velocity is determined
by tracking the zero crossing of the body mid-line and computing their
transition in the swimming direction. Note that the local velocity of these
zero crossings may vary considerably as they travel downstream along the body.
We define the mean wave velocity V (and the corresponding body slip
||/V) as the average velocity of the
zero crossings of the wave traveling from head to tail. |
Results |
|---|
|
TOPSummaryIntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
Reference motion pattern
The swimming motion is defined by Eqn 5 conforming with the motion pattern
used by Carling et al. (Carling et al.,
1998). The 2D simulations are compared with the results of Carling
et al. (Carling et al., 1998).
Subsequently we present the corresponding 3D simulations that serve as
reference case for the optimized motion patterns and allow for a
straightforward comparison between 2D and 3D simulations.
Efficient swimming mode
The swimming motion is defined by Eqn 3 and Eqn 4, and the free parameters
Ki, i=1,...,4 and tail are
obtained by maximizing the cost function presented in Eqn 7.
Fast swimming mode
The swimming motion is defined by Eqn 3 and Eqn 4, and the free parameters
Ki, i=1,...,4 and tail are
obtained by maximizing the cost function presented in Eqn 8.
Note that our results are normalized by the length L of the
swimming body, the undulation period T, and the fluid density .
Hence, the presented velocities correspond to body length per cycle (stride).
The velocity fields are plotted in the stationary frame where the free stream
velocity is zero. Unless otherwise noted we report results after the body has
reached its asymptotic swimming state.
Table 1 summarizes the
kinematics of the different cases.
|
Two-dimensional simulations of reference motion pattern
The simulation setup described in the preceding section was employed in
simulations of the flow past a 2D configuration. We use the motion pattern
defined in Eqn 5 in order to compare the results of the present simulations
with those reported in Carling et al.
(Carling et al., 1998). The
results of these simulations facilitate the identification of differences in
2D and 3D predictions of swimming velocity and wake patterns.
In the present simulations the thickness distribution of the 2D deformable body corresponds to the 3D cases, except that the thickness reduction from head to tail is linear instead of the quadratic function used in Eqn 1. The solution in 2D is computed on a C-type mesh corresponding to a slice of the mesh used in the 3D setup. With 2.1x104 grid cells distributed on the computational domain of radius R=3L and downstream length Ldstr=8L, the flow was shown to be largely resolved.
Fig. 2A depicts the unsteady
longitudinal and lateral velocities of the center of mass of the body. The
body accelerates from rest to an asymptotic mean forward velocity of
||=0.54 in about seven undulation cycles. The
velocity U|| varies slightly during a cycle while the
lateral velocity U has an amplitude of 0.04. The
tail beat amplitude is A=0.16 and the corresponding Strouhal number
is St=0.59. The wave velocity is V=0.73, which results in a
slip of ||/V=0.74. The force and moment
coefficients C||, C and
CM are plotted in Fig.
3 and converge to oscillation modes with zero mean and a constant
amplitude of 0.03, 0.04, and 0.03, respectively. The input power coefficient
reaches a mean value of
p=0.084.
|
The wake of the two dimensional simulation is visualized in Fig. 4A by velocity vectors and contours of vorticity. The wake consists of two counter-rotating vortices per undulation cycle shed at the tail. The vortex shedding starts and ends with the tail tip changing direction in every half cycle. The centers of the vortices, shed at the tail, have a distance of half a stride length and are well aligned in the downstream direction. Each pair of opposite sign vortices (a vortex dipole) induces a jet of fluid in the lateral direction. As the vortices detach from the tail, the lateral jets have a velocity magnitude of about 0.3.
|
The unsteady velocity of the center of mass of the body is depicted in
Fig. 2A. The body reaches an
asymptotic mean forward velocity of ||=0.40 and
oscillates with an amplitude of 0.01. The lateral velocity U has zero mean and an amplitude of 0.03. The wave velocity
V=0.73 is equal to the 2D case and consequently the slip is
||/V=0.55. The tail beat amplitude is
measured as A=0.15 with St=0.75.
Fig. 7 shows the net force and
moment coefficients C||, C
and CM oscillating around zero with amplitudes of 0.04,
0.06 and 0.03, respectively. The input power coefficient
p reaches the asymptotic
mean value of p=0.161.
|
;
Tytell and Lauder, 2004) with
strong lateral jets, separated by two same-sign vortices. The lateral fluid
velocity in the jet corresponding to the vortex pair shed at the tail is 0.35.
The 3D wake structure resembles a double vortex ring wake as outlined in
Müller et al. (Müller et al.,
2001). Further details of the 3D wake structure are presented in
the following section describing the flow characteristics of the optimized
swimming patterns.
Efficient swimming
A motion pattern corresponding to maximum swimming efficiency was
identified by the present optimization methodology. The optimization was
started at the initial search point K=(4.37,2.22,6.07,3.07) and
tail=1.71 obtained from preliminary investigations. The course
of the optimization process is plotted in
Fig. 5A,B. The optimization
strategy was stopped after 460 function evaluations. The best parameter set
found for Eqn 3, Eqn 4 is K=(3.34,1.67,6.28,6.28) and
tail=1.72. We note that K3 and
K4 converged to the maximum allowed value of 2,
indicating that swimming with higher curvature values towards the tail is
beneficial for efficient swimming. At the same time this tail motion is
combined with significant undulation in the anterior part of the body. The
resulting efficient swimming motion pattern is visualized in
Fig. 6.
|
|
has a
zero mean and an amplitude of 0.02. The wave velocity is V=0.60
resulting in a slip of ||/V=0.55 and
the tail undulation amplitude is A=0.11, resulting in
St=0.67. The net force and moment coefficients
C||, C and
CM are shown in Fig.
7 to oscillate around zero with amplitudes of 0.032, 0.052 and
0.028, respectively. The input power coefficient has an asymptotic mean value
of p=0.153.
Flow field and wake pattern are visualized in Figs
8,
9,
10 and
11. The 2D cross sections of
velocity and vorticity fields plotted in
Fig. 8A show a structure
similar to the one observed in the reference case, with a double row of vortex
pairs forming strong lateral jets. As the vortices detach from the tail, the
lateral jets have a velocity of about 0.3. The 3D structure of the flow
consists of a double row of vortex rings with pronounced secondary structures
that are visualized in Fig. 9A
by plotting isosurfaces of vorticity magnitude
||
||2. The formation of the
secondary structures is apparent in Fig.
10A showing isosurfaces of axial vorticity
|x|1. Axial vorticity is produced
primarily by the lateral motion of the body; it remains close to the surface
and is convected rather undisturbed along the body. As the local lateral
velocity changes sign, the vorticity is detached from the surface and leads to
the formation of complex secondary flow structures at the tail. A close-up
view of the secondary flow structures is given in
Fig. 11A showing isosurfaces
of vorticity magnitude ||||2 colored
by lateral vorticity y. The coloring of y
illustrates the boundary layer of the top and bottom part of the body being
separated at the tail. In addition, signs of a longitudinal jet formed by the
secondary structures can be identified.
|
|
|
|
The unsteady thrust and drag production along the body was investigated by recording the fluid forces acting on segments of the body. The fish surface was split into five segments of equal length along the body numbered from 1-5 from head to tail. Fig. 12A shows the unsteady fluid forces acting on the five segments over one undulation cycle. The anterior part of the body (segments 1 and 2) has no thrust contribution at any time. The major part of the thrust is produced at the tail (segment 5); however, a remarkable contribution comes also from the middle and posterior parts of the body (segments 3 and 4).
|
Fast swimming
The optimization process for maximal swimming velocity was initialized with
K=(1.29,0.52,5.43,4.28) and tail=1.52 obtained from
preliminary investigations. The course of the optimization process is shown in
Fig. 5C,D. The optimization
strategy was terminated after 460 function evaluations. The best parameter set
for Eqn 3, Eqn 4 after the given optimization time is
K=(1.51,0.48,5.74,2.73) and tail=1.44. These
parameters result in a motion pattern where the anterior part of the body
remains practically straight while the tail performs whip-like undulations
(cf. Fig. 6).
Fig. 2B depicts the unsteady
longitudinal and lateral velocities of the body. The asymptotic mean forward
velocity achieved is ||=0.47 and oscillates with
amplitude 0.02; the lateral velocity U has zero
mean and an amplitude of 0.05. The wave velocity is V=0.74 and the
resulting slip is 0.63. The tail undulation amplitude is A=0.14 and
the resulting Strouhal number St=0.59. The force and moment
coefficients C||, C and
CM are plotted in Fig.
7 and oscillate with zero mean and amplitudes of 0.043, 0.051 and
0.032, respectively. The input power coefficient reaches an asymptotic mean
value of p=0.124.
Flow field and wake pattern are visualized in Figs 8, 9, 10 and 11. The 2D signature of the wake is plotted in Fig. 8B and shows a similar structure as the reference case and the efficient swimming mode. The plots indicate stronger secondary flow structures present in the fast swimming resulting from the increased longitudinal and lateral velocities of the body. The increased tail undulation amplitude is also reflected in the velocity of 0.35 of the lateral jets.
The strong secondary flow structures become apparent in plots of the
vorticity magnitude ||||2 given in
Fig. 9B. The vortex rings shed
are stretched in longitudinal direction; the stretching seems to be
proportional to the swimming velocity increased by a factor of 1.4 when
compared to the efficient swimming mode.
Fig. 10B shows isosurfaces of
axial vorticity |x|1; the structures
and the shedding of x resemble the efficient case with
increased vorticity strength at the tail. The close-up view of the flow
structure at the tail in Fig.
11B illustrates the complex interaction of the tail tip vortices
and the separation of the boundary layer of top and bottom part of the
body.
The segment-wise investigation of fluid forces acting on the body is depicted in Fig. 12B. As for the efficient swimming, the first two segments only produce drag. The third and fourth segment contribute both to instantaneous thrust and drag; the thrust contributions only exceed slightly the drag obtained by averaging over an entire undulation cycle. Hence, the tail (segment 5) is almost solely responsible for the thrust produced during the fast swimming mode.
|
Discussion |
|---|
|
TOPSummaryIntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
We compare the results of the present simulations with the experimental
investigations reported elsewhere
(Müller et al., 2001;
Tytell and Lauder, 2004), and
previous 2D simulations (Carling et al.,
1998). We remark that the Reynolds number in the experimental
studies is a factor of 5-30 higher than in our simulations. Nevertheless, the
wake patterns obtained in our simulations show strong similarities with the
experimental flow fields and the swimming velocities and slip values achieved
in the simulations lie well within the range of values reported in the
experimental studies. Our findings provide valuable quantitative information
that further elucidates the mechanisms of anguilliform propulsion and
correlates them with physiological objectives such as migratory or burst
swimming.
Swimming kinematics and hydrodynamic forces
2D reference case versus results by Carling et al.
The swimming velocity and the slip of the present simulations have small
differences from the values reported in Carling et al.
(Carling et al., 1998). At the
same time the flow structures observed in the two simulations are drastically
different and we do not venture into further comparison of kinematics and
fluid forces.
3D cases versus experimental studies
The mean swimming velocity || obtained in the
present 3D simulations lies well within the range of 0.26-0.5 body lengths
cycle-1 reported in experimental studies
(Müller et al., 2001;
Tytell, 2004;
Tytell and Lauder, 2004). The
values for the slip ||/V of 0.55 and
0.63 obtained for the efficient and fast swimming modes, respectively, are
slightly below the experimental values ranging between 0.6 and 0.73. The
largest discrepancy in the kinematics of our 3D simulations when compared to
the kinematics reported in the experimental studies, is found in the values of
the Strouhal number. The values for St obtained in our simulations
(0.67 and 0.59) are clearly above the values reported by Müller et al.
(Müller et al., 2001) and
Tytell (Tytell, 2004) and the
range of 0.2-0.4 identified for efficient oscillatory propulsion of swimming
and flying animals (Taylor et al.,
2003). This discrepancy may be attributed to the reduced Reynolds
number and the fixed undulation frequency used in the present study and will
be addressed in future investigations.
The comparison of experimental and computational results in terms of force
and total input power is difficult due to the limited available data and the
different models used to infer these forces from 2D cross sections of the
wake. In Tytell and Lauder (Tytell and
Lauder, 2004) coefficients for mean lateral force and total input
power are computed based on vortex ring momentum as estimated from the 2D PIV
data. The values provided in (Tytell and
Lauder, 2004) are
||=0.097 and
p=0.023. The experimental
value for || is a factor of 2.4 larger than the present 2D result and
about 1.5 times larger than the values obtained in our 3D simulations; the
experimental value for p
is a factor of 4 smaller than the present 2D results and a factor 5-7 smaller
than the present 3D results. The lower values obtained in the experiments may
be attributed to the employed vortex ring models that may not be able to
account for the effect of secondary flow structures.
2D versus 3D reference motion
The swimming velocity obtained for the 2D simulation of the reference
motion pattern is a factor of 1.35 larger than in the corresponding 3D
simulation, which is also reflected in the slip and the Strouhal number. The
force coefficients C|| and C, and the moment coefficient CM of the 2D and
3D case match surprisingly well. The limited prediction capabilities from 2D
swimming models, however, become evident in the comparison of the wake
structure described below.
Comparison of the optimized swimming modes
The different characteristics of the optimized motion patterns become
evident in Fig. 6. The striking
feature of the efficient swimming mode is that the anterior part (except the
head) has a very small undulation amplitude and is kept well aligned with the
swimming direction. Compared to the reference mode, the tail amplitude is a
factor of 1.3 smaller. In addition, the angle of the tail tip, as the tail
changes direction, is reduced in the efficient swimming motion supporting
smoother vortex shedding. In contrast to the efficient and the reference
motion patterns, in the fast swimming motion the anterior half of the body
hardly bends. Meanwhile, the motion of the tail resembles the pattern of the
efficient swimming case with increased undulation amplitudes. The swimming
velocity of the fast motion pattern exceeds the velocity of the efficient
pattern by 42% while the total input power is increased by a 130%. The
efficient swimming mode therefore needs 1.6 times less energy to swim a given
distance, validating the proposed fitness functions used in the optimization
process.
A comparison of the unsteady drag and thrust production along the body
(Fig. 12) reveals that for the
fast swimming mode thrust is almost exclusively produced at the tail, while
for the efficient swimming mode, the middle part of the body has a
considerable contribution to the thrust. The results for the efficient
swimming pattern support the hypothesis
(Blickhan et al., 1992;
Müller et al., 2001;
Tytell and Lauder, 2004), that
anguilliform swimmers produce thrust not only with their tail but also with
parts of the body anterior to the tail. As it is shown in
Fig. 12, the first two body
segments are largely responsible for drag generated as the flow first
encounters the body. Segments 3 and 4 serve to adjust the vorticity propagated
from upstream associated with a net thrust in the case of efficient swimming;
the shedding of the vorticity at the tail (segment 5) results in thrust.
Wake morphology
2D simulation
The wake structure obtained in our 2D simulation disagrees with reported
results (Carling et al., 1998)
even though we use a similar geometry, the same motion pattern and a Reynolds
number of Re=3857 that closely matches the reported value of
Re=3840 (computed based on body length L). The reasons for
the conflicting results remain unclear. Our results, however, are supported by
similar wake patterns reported in previous 2D fluid dynamic simulations of
swimming motions (Flanagan et al.,
2004; Pedley and Hill,
1999).
The striking difference of the wake obtained in our 2D simulations when compared to 2D-sections of the 3D results at the mid-plane, is that the vortices shed in every half period do not break up and are well aligned in the downstream direction.
3D simulations
The wake patterns obtained for our 3D simulations are in good qualitative
agreement with 2D PIV measurements of real anguilliform swimmers. The
experimental works of Müller et al.
(Müller et al., 2001) and
Tytell and Lauder (Tytell and Lauder,
2004) report wakes with the characteristic lateral jets observed
in our simulations. In contrast to the shedding mechanism described in Tytell
and Lauder (Tytell and Lauder,
2004), our results indicate that the break-up of the primary
vortices occurs as they are shed at the tail. The early break-up of the
primary vortex observed in our results can be recognized in the flow fields
provided in Müller et al.
(Müller et al., 2001).
The wake structure and the shedding mechanism are similar for all three motion
patterns that have been investigated in the current study. The main
differences in the vortical structures of the different swimming patterns are
concentrated in the strength of the secondary vortices.
The 3D vortex ring structures in the wake obtained in our simulations (Figs
9 and
11) are in good agreement with
the predictions sketched in Müller et al.
(Müller et al., 2001) and
Lauder and Tytell (Lauder and Tytell,
2006). The shedding of the double vortex rings is accompanied by
complex secondary flow structures that are generated by the transversal motion
of the body, convected downstream, and eventually shed in the wake (Figs
10,
11). We presume that the
continuous caudal fins often present in natural anguilliform swimmers play an
important role in modifying these secondary flow structures in order to
increase swimming efficiency.
In the fast swimming case the secondary structures are much stronger, and the vortex rings are stretched in the swimming direction, forming jets that are approximately 1.4 times wider than the ones observed in the efficient swimming case. In experiments, the velocity magnitude of the lateral jets, just after being detached from the tail, is 20-40% of the mean swimming velocity. In contrast, the jet velocity obtained in our simulations is within 75-90% of the swimming velocity. These values are significantly higher than the ones reported in the experimental results and they are also reflected in an increased Strouhal number.
Regarding the direction of the lateral jets, Tytell and Lauder
(Tytell and Lauder, 2004)
report that jets with a non-negligible downstream component are an indication
of acceleration of the anguilliform swimmer. This observation is consistent
with the results of the present simulations where in the initial acceleration
phase of the simulation the jets point downstream with an angle of about
30° to the lateral direction. The jet angle smoothly decays close to zero
when the asymptotic swimming velocity is reached.
For the present simulation setup, we find that the optimization of swimming
efficiency does not result in a wake pattern with a reverse von
Kármán vortex street as observed in Gray (Gray, 1968), plate I
p. 228, and postulated by Müller et al.
(Müller et al., 2001). We
suppose that in addition to the motion pattern, the body geometry and the
Reynolds number may have a considerable influence on the formation and
shedding of the wake. Continuous caudal fins observed in many species of
anguilliform swimmers presumably serve to control the formation of the
secondary vortex structures. Simulations with increased Reynolds number and
investigations of the influence of the body geometry (back and anal fins,
etc.) and body motion will help to clarify this issue.
Concluding remarks
We have investigated optimized patterns of anguilliform swimming using 3D
simulations. The wake structures of the present simulations are consistent
with several experimental observations. The fast and the efficient swimming
mode both shed a double row of vortex rings responsible for the strong lateral
jets observed in the wake. The results provide quantification of the vortex
formation and shedding processes and enable the identification of the portions
of the body that are responsible for the majority of thrust in anguilliform
swimming. In burst swimming the tail is responsible for the majority of the
thrust, while in efficient swimming the anterior part of the body also
contributes to the thrust.
,
s(s,t)
i
, 
||2S)
=F/(||2S)
(||2LS)
||S3)
L||/µ
||
||
position and velocity of fish
surface
|
Appendix |
|---|
|
TOPSummaryIntroductionMaterials and methodsResultsDiscussionAppendixReferences |
|---|
 | (A1) |
and a subsequent transformation to assure a vanishing rotational impulse
and the center of mass coinciding with O' in the system
(O',x',y').
and 
denote tangent and normal vectors of the curve, respectively.
The position and orientation (and their change in time) of the body fixed
system
(O',x',y')
with respect to the inertial system (O,x,y,z) is obtained by
solving coupled system. The above steps to ensure that the net lateral force
and torque provided by hydrodynamic forces on the whole body is equal to the
rate of change of its lateral transitional and angular momentum is usually
called recoil correction (Cheng
and Blickhan, 1994; Hess and
Videler, 1984). The two coordinate systems are illustrated in
Fig. A1.
|
(B) Coupled solution procedure
In the following we give a detailed description of the procedure used to
solve the coupled system of the Navier-Stokes equations for the fluid and
Newton's equations for the undulating body. The method is based on the
Improved Serial Staggered (ISS) procedure
(Farhat and Lesoinne, 2000).
The solution of the fluid phase is shifted by half a time step compared to the
solution of the body dynamics. We use to denote the position of fluid
grid nodes and x for the position of the surface of the deforming body.
xc denotes the position of the center of mass of the body
and c the global angle in yaw direction. The discrete time is
denoted in bracketed superscripts, and n is the current time step.
The iteration steps are as follows:
(0) Initialize the simulation with the body at rest floating in the fluid
at rest by setting
(-)=x(0) according to the
starting position and geometry of the body and
c(0)=0.
FOR n=1,..., nend DO
(1) Compute the deformation velocity
of the surface S of the deforming body (prescribed by the deforming
centerline). Set the fluid grid velocity
(n)
equal to the velocity of the body
(n)
on the fluid-body interface S:
 | (A2) |
(2) Update the fluid dynamic grid on the fluid-body interface:
 | (A3) |
The velocity of interior nodes is determined by the grid deformation procedure described below.
(3) Advance the transient flow solution from n- to
n+.
(4) Compute the fluid force
F(n+) and torque
M (n+)z
acting on the deforming body with respect to its center of mass according to
Eqn 15 and Eqn 16. Smooth the computed fluid forces with a low pass filter
implemented as:
 |
 | (A4) |
was set to 0.2.
(5) Advance the center of mass xc of the fish
with the filtered forces
(n+) and
the moment
z(n+)
according to:
 |
 |
 |
 |
 |
where the inertial moment about yaw axis Iz of the
deforming body and its temporal derivative 
z can be
directly computed from the discretized geometry and the prescribed motion
pattern.
|
