SUMMARY
The hydrodynamics of anguilliform swimming motions was investigated using threedimensional simulations of the fluid flow past a selfpropelled body. The motion of the body is not specified a priori, but is instead obtained through an evolutionary algorithm used to optimize the swimming efficiency and the burst swimming speed. The results of the present simulations support the hypothesis that anguilliform swimmers modify their kinematics according to different objectives and provide a quantitative analysis of the swimming motion and the forces experienced by the body.
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 threedimensional fluidbody 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 selfpropelled body. Finally, the optimization procedure helps to identify, in a systematic fashion, links between swimming motion and biological function.
Introduction
Anguilliform swimming is the primary mode of locomotion for numerous aquatic species across a range of diverse taxa. Anguilliform swimmers propel themselves forward by propagating waves of curvature towards the posterior of the body, and this type of locomotion is widespread among species ranging in scale from nematodes to eels (Müller et al., 2001). 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 threedimensional simulations of a selfpropelled eellike body immersed in a viscous fluid. In addition, twodimensional simulations are compared with existing related works (Carling et al., 1998), in order to assess the role of two and threedimensional 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' (CMAES) (Hansen and Kern, 2004; Hansen and Ostermeier, 2001). The CMAES 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 twodimensional 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 selfpropelled 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
The computational model solves the threedimensional NavierStokes equations for the incompressible, viscous flow past a deforming, anguilliform body. The time dependent motion of the body is defined by the twodimensional deformation of its midline, allowing for yaw, but no pitch or roll motions. The fluidbody interaction includes forces in all three coordinate directions and the torque perpendicular to the plane of the deformation of the midline (yaw). The other two components (pitch and roll) of the torque are neglected in order to reduce computational complexity of the system. In the following sections we describe the geometrical model, the parameterization of the body deformation, and the procedure to solve the fluid dynamic problem including the interaction of the flow with the deforming body.
Geometrical model
The threedimensional 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 midline 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 w_{h}=s_{b}=0.04L, s_{t}=0.95L and w_{t}=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 threedimensional (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 I_{z}, the body is discretized into ellipsoid disks of constant density along the midline s with a step size of Δs=0.001L. The deformation of the body is composed from the motions of the individual segments.
Parameterization of the motion
Several works have proposed a description of the anguilliform motion by providing an analytical expression for the lateral displacement of the midline. Pedley and Hill (Pedley and Hill, 1999) 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 midline 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 midline κ_{s} as: (3)
where K(s) is the cubic spline through the m interpolation points K_{i}, 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 K_{i} and τ_{tail} are determined by the optimization procedure. The description of the midline 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 y_{s} describes the lateral displacement of the midline 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 selfpropelled 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, K_{i} 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' (CMAES) (Hansen and Kern, 2004; Hansen and Ostermeier, 2001). The efficiency of CMAES 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 K_{i}, 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 P̄_{total} that is being transformed into forward motion: (6)
In the present optimization, the period T is constant and we integrate P̄_{total} over the undulation cycle in order to obtain the amount of work W_{cycle} 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 P̄_{forward motion} is related to the kinetic energy of the forward motion of the body E_{kin}=½mŪ_{∥}^{2}. Hence, we characterize the swimming efficiency by the ratioη∼ E_{kin}/W_{cycle}, 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 P̄_{max}=2.0×10^{3} and P_{max}=4.8×10^{3} based on preliminary computational experiments. The scaling factors where set to c_{1}=10^{6} and c_{2}=10^{8}. The maximum curvature of the midline 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 NavierStokes 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 u_{f}, and the body velocity u_{c}: (12)
The motion of the selfpropelled body is in turn described by Newton's equations of motion: (13) (14)
F and M_{z} are the fluid force and yaw torque acting on the body surface, x_{c} is the position of the center of mass of the body, its global angular velocity, m the total mass, and I_{z} 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 M_{z} 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 NavierStokes equations were solved using the commercial package STARCD v. 3.15. STARCD computes flows on arbitrary LagrangianEulerian grids by solving an additional space conservation equation. Eqn 10 and Eqn 11 are discretized using a finite volume approach with firstorder discretization in time and secondorder in space. The solution of Newton's equations of motion and the resulting grid deformation and motion are implemented in user subroutines linked to STARCD. 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 arclengthbased transfinite interpolation is used to propagate the deformations of the corner points (and the prescribed surface deformations at the fluidbody boundary) to the interior grid cells. Fig. 1B illustrates the deforming grid. The simulation setup and the flowstructure 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.4×10^{4}, body length L=1, densityρ _{fluid}=ρ_{body}=ρ=1, and undulation frequency f=1. The resulting Reynolds number Re=ρLŪ_{∥}/μ is in the range 24003900, 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 timeaveraged 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 3×10^{5} cells and a time stepΔ t=5×10^{4} (corresponding to 2000 time steps per undulation cycle) was shown to accurately resolve the flow. During the optimization, a grid with 7×10^{4} cells andΔ t=5×10^{3} was used in order to reduce the computational cost while still allowing a satisfactory computation of the fluidbody interaction and the resulting dynamics. The results of the optimization were verified by simulations of high resolution.
The fluid forces acting on the body are measured in nondimensional coefficients C_{∥}=F_{∥}/(½ρŪ_{∥}^{2}S) and C_{ ⊥}=F_{⊥}/(ρŪ_{∥}^{2}S) 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 C_{M}=M_{z}/(½ρŪ_{∥}^{2}LS) and C_{P}=P_{total}/(½ρŪ_{∥}^{3}S). 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 midline 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
We present results of the simulations, including swimming kinematics, hydrodynamic forces, and wake morphology for a reference swimming pattern and for the efficient and fast swimming.
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 K_{i}, 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 K_{i}, 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.
Twodimensional 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 Ctype mesh corresponding to a slice of the mesh used in the 3D setup. With 2.1×10^{4} grid cells distributed on the computational domain of radius R=3L and downstream length L_{dstr}=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 C_{M} 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 C̄_{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 counterrotating 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.
3D simulations of reference motion pattern
The 3D simulations of the flow past a selfpropelled body with the motion pattern defined in Eqn 5 constitute the reference case of this study. It allows direct comparison of 2D and 3D results and a quantification of the relative success of the optimization processes.
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 C_{M} oscillating around zero with amplitudes of 0.04, 0.06 and 0.03, respectively. The input power coefficient C̄_{p} reaches the asymptotic mean value of C̄_{p}=0.161.
The 2D signature of the wake is given in Fig. 4B in terms of contours of the vorticity component normal to the image plane along with the corresponding components of the velocity field. The vorticity shed in every half tail beat cycle breaks up into two vortices, pairing up with vortices produced during the preceding and the succeeding half tail beat cycle, respectively. The resulting wake shows the characteristic pattern observed in experiments (Müller et al., 2001; Tytell and Lauder, 2004) with strong lateral jets, separated by two samesign 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 K_{3} and K_{4} 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.
Fig. 2B depicts the unsteady longitudinal and lateral velocities of the body; the reached asymptotic mean forward velocity is U_{∥}=0.33 and is oscillating with an amplitude of 0.007. The lateral velocity U_{ ⊥} 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 C_{M} 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 C̄_{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 closeup 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 15 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 whiplike 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 C_{M} 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 C̄_{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 closeup 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 segmentwise 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
The 3D simulations provide us with quantitative information regarding the kinematics of the optimized body motions, the vorticity in the wake as well as the distribution of fluid forces acting on the surface of the body.
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 530 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.260.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.20.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 C̄_{⊥}=0.097 and C̄_{p}=0.023. The experimental value for C̄_{⊥} 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 C̄_{p} is a factor of 4 smaller than the present 2D results and a factor 57 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 C_{M} 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 2Dsections of the 3D results at the midplane, 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 breakup of the primary vortices occurs as they are shed at the tail. The early breakup 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 2040% of the mean swimming velocity. In contrast, the jet velocity obtained in our simulations is within 7590% 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 nonnegligible 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.
 List of symbols
 ^
 filtered quantities, cf. Eqn 20
 —
 time averaged quantities
 ζ, ζ
 position and velocity of computational grid
 η
 propulsive efficiency
 K_{i}
 curvature amplitude at the ith spline interpolation point
 κ_{s}(s,t)
 time dependent curvature of midline
 μ
 viscosity
 ω
 vorticity
 ρ
 density
 viscous stress tensor
 τi
 time shift of oscillating curvature at interpolation point i
 Φ,
 mean angle and angular velocity of body with respect to inertial system
 A
 tail undulation amplitude
 C_{∥}=F_{∥}/(½ρŪ_{∥}^{2}S)
 net force coefficient in swimming direction
 C_{ ⊥}=F_{⊥}/(½ρŪ_{∥}^{2}S)
 lateral net force coefficient
 C_{M}=M_{z}/½(ρŪ_{∥}^{2}LS)
 net yaw torque coefficient
 C_{P}=P_{total}/(½ρŪ_{∥}S^{3})
 total input power coefficient
 f
 undulation frequency
 f_{η}
 cost function for optimizing swimming efficiency
 f_{U}
 cost function for optimizing swimming velocity
 F
 net fluid force acting on center of mass
 g
 body force per unit mass
 I_{z}
 inertial moment about yaw axis
 L
 body length
 L_{dstr}
 length of downstream section of computational domain
 m
 mass
 M_{z}
 net yaw fluid torque acting about center of mass
 n
 outer surface normal vector
 P, P̄
 (time averaged) output power
 R
 radius of computational domain
 Re=ρLŪ_{∥}/μ
 Reynolds number
 s
 coordinate along the midline of the body
 S
 surface of the 3D body (circumference for 2D geometry)
 St=2fA/Ū_{∥}
 Strouhal number
 t
 time
 T=1/f
 cycle time
 U_{∥},Ū_{∥}
 (time averaged) forward swimming velocity
 U_{ ⊥}
 lateral velocity of center of mass
 V
 wave velocity
 W
 work per undulation cycle
 x,
 ẋ position and velocity of fish surface
Appendix
Additional materials and methods
(A) Transformation of the midline description
The description of the midline according to Eqn 3 and Eqn 5 has to be transformed to the 3D Cartesian inertial system used by the flow solver to construct the computational grid enclosing the deforming body. First, a representation r′(s) in the 2D Cartesian system (O′,x′,y′) is computed in which the center of mass of the body coincides with O′ and the rotational impulse of the deforming body amounts to zero (taking into account the distribution of mass along the body resulting from Eqn 1 and Eqn 2). This is done by integrating the Frenet equations for a planar curve with arbitrary initial conditions: (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 NavierStokes 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. x_{c} 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,..., n_{end} 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 fluidbody interface S: (A2)
(2) Update the fluid dynamic grid on the fluidbody 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 x_{c} of the fish with the filtered forces F̂^{(n+½)} and the moment M̂_{z}^{(}^{n}^{+½)} according to:
where the inertial moment about yaw axis I_{z} of the deforming body and its temporal derivative İ_{z} can be directly computed from the discretized geometry and the prescribed motion pattern.
(C) Falling sphere test case
The implementation of the fluidbody coupling procedure was validated on a falling sphere test case at Re_{d}=(ρ_{l}dU)/μ=100. A rigid sphere with ρ_{s}>ρ_{l} is released from rest and accelerates until it reaches its asymptotic falling velocity. The sphere diameter d is set to d=1.0 and the density and the kinematic viscosity is chosen as ν=μ/ρ_{l}=0.01 to aim for an asymptotic falling velocity U=1. For Re_{d}=100 the drag coefficient for the flow past a sphere is C_{D}=1.1 (Johnson and Patel, 1999). To determine the gravity constant g and ρ_{s} we write F_{grav}F_{buoy}=F_{D}=C_{D}[½ρfU^{2}(πd^{2}/4)], and from this we getρ _{s}/ρ_{t}=1+C_{D}(3U^{2}/4g). We choose g=20 and thus the density of the sphere is set toρ _{s}=1.041. The grid used is an OO type with radius 15 and 100×40×100 cells with exponential clustering towards the wall, according to the grid used by Johnson and Patel (Johnson and Patel, 1999). The time step was set to Δt=0.001 and the same numerical method and solver settings where used as for the swimming simulations described above. The course of the falling velocity is plotted in Fig. A2A. At time t=20 an asymptotic falling velocity of U=1.006 is reached which matches the predicted value pertaining to the chosen parameters very well. Fig. A2B shows vorticity contours at time t=20 that agree with the results in literature (Johnson and Patel, 1999).
(D) Convergence in space and time
Spatial and temporal convergence was tested using the motion pattern (Eqn 5). The convergence order of the presented computational setup was determined by computing the norm of the relative error of the unsteady longitudinal force with respect to the bestresolved case over the last undulation cycle: (A6)
The convergence order in space was obtained using the grids summarized in Table A1 ranging from 7×10^{4} to 4.2×10^{5} and a time stepΔ t=0.0005. The size of the computational domain was constant for all cases with a radius of R=2L and total length of 8L. Convergence order in time was obtained using the grid with 3×10^{5} cells and Δt={0.0005, 0.001, 0.002, 0.004}. The resulting graphs are given in Fig. A3A,B and show that the coupled simulation setup is second order accurate with respect to the grid size and first order accurate with respect to time. Fig. A3C,D show the longitudinal force F_{∥} and the forward swimming velocity U for the different grids with a timestepΔ t=5×10^{4}. To resolve the flow, a grid with 3×10^{5} cells proved to be suitable. The time step wasΔ t=5×10^{4}, which corresponds to 2000 time steps per undulation cycle. For the optimizations cases, a grid with 7×10^{4} grid cells and a time step ofΔ t=5×10^{3} were used to reduce the computational cost to a viable amount. This setup still allows a satisfactory computation of the fluidbody interaction and the resulting dynamics.
ACKNOWLEDGEMENTS
We wish to acknowledge enlightening discussions with George Triantafyllou (National Technical University of Athens, Greece) and the valuable input of the anonymous reviewers that considerably improved the manuscript. We thank Philippe Chatelain, Jens Walther and Michael Bergdorf for valuable discussions and support regarding this project. Financial support of the Swiss National Science Foundation SNF is acknowledged.
FOOTNOTES

Supplementary material available online at http://jeb.biologists.org/cgi/content/full/209/24/4841/DC1
 © The Company of Biologists Limited 2006