## SUMMARY

The hydrodynamics of anguilliform swimming motions was investigated using three-dimensional simulations of the fluid flow past a self-propelled 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 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.

## 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 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

The computational model solves the three-dimensional Navier-Stokes equations for the incompressible, viscous flow past a deforming, anguilliform body. The time dependent motion of the body is defined by the two-dimensional deformation of its mid-line, allowing for yaw, but no pitch or roll motions. The fluid-body interaction includes forces in all three coordinate directions and the torque perpendicular to the plane of the deformation of the mid-line(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

*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:

*w*

_{h}=

*s*

_{b}=0.04

*L,s*

_{t}=0.95

*L*and

*w*

_{t}=0.01

*L*. The height

*h*(

*s*) is set as an elliptical curve with the two half axis

*a*=0.51

*L*and

*b*=0.08

*L*:

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 *I*_{z}, the body is discretized into ellipsoid disks of constant density along the mid-line *s* with a step size of Δ*s*=0.001*L*. 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 mid-line. 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 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.

_{s}as:

*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:

The parameters *K*_{i} 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).

*L*=1 to:

where *y*_{s} 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, *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' (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 *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.

*P̄*

_{total}that is being transformed into forward motion:

*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:

**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:

*Ū*

_{∥}and two terms that penalize high mean and peak input power requirements:

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 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

**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}**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,

*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:

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.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 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 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 fluid-body 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*=2*fA*/*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

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.

. | 2D reference . | 3D reference . | Efficient . | Fast . |
---|---|---|---|---|

Ū_{∥} | 0.54 | 0.40 | 0.33 | 0.47 |

V | 0.73 | 0.73 | 0.60 | 0.74 |

Slip | 0.74 | 0.55 | 0.55 | 0.63 |

A | 0.16 | 0.15 | 0.11 | 0.14 |

St | 0.59 | 0.75 | 0.67 | 0.59 |

C_{P} | 0.084 | 0.161 | 0.153 | 0.124 |

. | 2D reference . | 3D reference . | Efficient . | Fast . |
---|---|---|---|---|

Ū_{∥} | 0.54 | 0.40 | 0.33 | 0.47 |

V | 0.73 | 0.73 | 0.60 | 0.74 |

Slip | 0.74 | 0.55 | 0.55 | 0.63 |

A | 0.16 | 0.15 | 0.11 | 0.14 |

St | 0.59 | 0.75 | 0.67 | 0.59 |

C_{P} | 0.084 | 0.161 | 0.153 | 0.124 |

### 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.1×10^{4} grid cells distributed on the computational domain of radius *R*=3*L* and downstream length *L*_{dstr}=8*L*, 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 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.

### 3D simulations of reference motion pattern

The 3D simulations of the flow past a self-propelled 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 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 *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. 9Aby 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. 12Ashows 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 *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 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

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 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|*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 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 *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 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.

- ^
filtered quantities, cf. Eqn 20

- —
time averaged quantities

- ζ, ζ
position and velocity of computational grid

- η
propulsive efficiency

*K*_{i}curvature amplitude at the

*i*th spline interpolation point- κ
_{s}(*s,t*)time dependent curvature of mid-line

- μ
viscosity

- ω
vorticity

- ρ
density

- \({\bar{{\sigma}}}\)
viscous stress tensor

- τi
time shift of oscillating curvature at interpolation point

*i* - Φ,\({\dot{{\phi}}}\)
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 mid-line of the body

- S
surface of the 3D body (circumference for 2D geometry)

*St=*2*fA*/*Ū*_{∥}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 mid-line description

**′(**

*r**s*) in the 2D Cartesian system(

**′,**

*O***′,**

*x***′)is computed in which the center of mass of the body coincides with**

*y***′ 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:**

*O*and a subsequent transformation to assure a vanishing rotational impulse and the center of mass coinciding with ** O**′ in the system(

**′,**

*O***′,**

*x***′).ξ and ν denote tangent and normal vectors of the curve, respectively. The position and orientation (and their change in time) of the body fixed system(**

*y***′,**

*O***′,**

*x***′)with respect to the inertial system (**

*y***) 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**

*O,x,y,z**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. **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

**ẋ**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*:

The velocity of interior nodes is determined by the grid deformation procedure described below.

(3) Advance the transient flow solution from *n*-½ to *n*+½.

**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:

α was set to 0.2.

*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 fluid-body 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}[½ρf*U*^{2}(π*d*^{2}/4)],and from this we getρ _{s}/ρ_{t}=1+*C*_{D}(3*U*^{2}/4*g*). We choose *g*=20 and thus the density of the sphere is set toρ _{s}=1.041. The grid used is an O-O 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

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*=2*L* and total length of 8*L*. 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 fluid-body interaction and the resulting dynamics.

. | n_{cell}
. | n_{i} × n_{j} × n_{k}
. | n_{i}^{B}
. | h_{head}
. | h_{tail}
. | R
. | L_{dstr}
. |
---|---|---|---|---|---|---|---|

Grid 1 | 70 × 10^{3} | 128 × 29 × 19 | 48 | 3.0 × 10^{−3} | 4.9 × 10^{−3} | 2L | 5L |

Grid 2 | 200 × 10^{3} | 133 × 41 × 37 | 52 | 2.5 × 10^{−3} | 3.0 × 10^{−3} | 2L | 5L |

Grid 3 | 301 × 10^{3} | 157 × 47 × 41 | 60 | 1.8 × 10^{−3} | 2.9 × 10^{−3} | 2L | 5L |

Grid 4 | 423 × 10^{3} | 171 × 53 × 47 | 64 | 1.8 × 10^{−3} | 2.3 × 10^{−3} | 2L | 5L |

. | n_{cell}
. | n_{i} × n_{j} × n_{k}
. | n_{i}^{B}
. | h_{head}
. | h_{tail}
. | R
. | L_{dstr}
. |
---|---|---|---|---|---|---|---|

Grid 1 | 70 × 10^{3} | 128 × 29 × 19 | 48 | 3.0 × 10^{−3} | 4.9 × 10^{−3} | 2L | 5L |

Grid 2 | 200 × 10^{3} | 133 × 41 × 37 | 52 | 2.5 × 10^{−3} | 3.0 × 10^{−3} | 2L | 5L |

Grid 3 | 301 × 10^{3} | 157 × 47 × 41 | 60 | 1.8 × 10^{−3} | 2.9 × 10^{−3} | 2L | 5L |

Grid 4 | 423 × 10^{3} | 171 × 53 × 47 | 64 | 1.8 × 10^{−3} | 2.3 × 10^{−3} | 2L | 5L |

For the optimization procedure, grid 1 was used, and the final simulations where carried out on grid 3.

*n*_{cell} denotes the total number of grid cells, *n*_{i,j,k} the number of cells in the three grid directions,and *n*_{i}^{B} indicates the number of cells along the body. *h*_{head} and *h*_{tail} denote the height of the first cell layer on the surface of the body at head and tail,respectively.

## 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.

## References

**Blickhan, R., Krick, C., Zehren, D., Nachtigall, W. and Breithaupt, T.**(

**Carling, J., Williams, T. L. and Bowtell, G.**(

**Cheng, J. Y. and Blickhan, R.**(

**Ekeberg, O.**(

**Farhat, C. and Lesoinne, M.**(

**Flanagan, P., Hotchkiss, R. and Stock, D.**(

**Gray, J.**(

**Hansen, N. and Kern, S.**(

**Hansen, N. and Ostermeier, A.**(

**Hess, F. and Videler, J. J.**(

*Pollachius virens*): a dynamic analysis of bending moments and muscle power.

**Ijspeert, A. J. and Kodjabachia, J.**(

**Johnson, A. T. and Patel, V. C.**(

**Kern, S., Müller, S. D., Hansen, N., Büche, D.,Ocenasek, J. and Koumoutsakos, P.**(

**Lauder, G. V. and Tytell, E. D.**(

**Milano, M. and Koumoutsakos, P.**(

**Müller, U. K., Smit, J., Stamhuis, E. J. and Videler, J. J.**(

*Anguilla anguilla*).

**Pedley, T. J. and Hill, S. J.**(

**Schmidt, J.**(

**Schultz, W. W. and Webb, P. W.**(

**Taylor, G. K., Nudds, R. L. and Thomas, A. L. R.**(

**Tsai, H. M., Wong, A. S. F., Cai, J., Zhu, Y. and Liu, F.**(

**Tytell, E. D.**(

**Tytell, E. D. and Lauder, G. V.**(

**van Ginneken, V., Antonissen, E., Müller, U. K., Booms, R.,Eding, E., Verreth, J. and van den Thillart, G.**(