## SUMMARY

We employ numerical simulation to investigate the hydrodynamics of carangiform locomotion as the relative magnitude of viscous and inertial forces, i.e. the Reynolds number (*Re*), and the tail-beat frequency,i.e. the Strouhal number (*St*), are systematically varied. The model fish is a three-dimensional (3D) mackerel-like flexible body undulating with prescribed experimental kinematics of carangiform type. Simulations are carried out for three *Re* spanning the transitional and inertial flow regimes, *Re*=300 and 4000 (viscous flow), and ∞ (inviscid flow). For each *Re* there is a critical Strouhal number, *St*^{*}, at which the net mean force becomes zero, making constant-speed self-propulsion possible. *St*^{*} is a decreasing function of *Re* and approaches the range of *St* at which most carangiform swimmers swim in nature (*St*∼0.25) only as *Re* approaches infinity. The propulsive efficiency at *St*^{*} is an increasing function of *Re* while the power required for swimming is decreasing with *Re*. For all *Re*, however, the swimming power is shown to be significantly greater than that required to tow the rigid body at the same speed. We also show that the variation of the total drag and its viscous and form components with *St* depend on the *Re*. For *Re*=300, body undulations increase the drag over the rigid body level, while significant drag reduction is observed for *Re*=4000. This difference is shown to be due to the fact that at sufficiently high *Re* the drag force variation with *St* is dominated by its form component variation, which is reduced by undulatory swimming for *St*>0.2. Finally, our simulations clarify the 3D structure of various wake patterns observed in experiments –single and double row vortices – and suggest that the wake structure depends primarily on the *St*. Our numerical findings help elucidate the results of previous experiments with live fish, underscore the importance of scale (*Re*) effects on the hydrodynamic performance of carangiform swimming, and help explain why in nature this mode of swimming is typically preferred by fast swimmers.

## INTRODUCTION

The majority of fishes use body/caudal fin (BCF) undulations for propulsion. Only about 12% of 450 extant fish families discussed by Nelson(Nelson, 1984) use undulations of median or pectoral fins (MPF) as their routine propulsive mode(Videler, 1993). However, many more species use median or pectoral fins for stability and maneuverability,especially at low speeds, while many MPF swimmers use BCF mode to reach high speeds during escape reactions (Videler,1993). Within the general MPF and BCF propulsion modes there are several other sub-modes of steady straight-line swimming, as initially classified (Breder, 1926) and later expanded (Lindsey, 1978)(for recent reviews, see Sfakiotakis et al., 1999; Lauder and Tytell,2006).

*Scomber scombrus*(Teleostei: Scombridae). Fishes in the Scombridae family are characterized by their streamlined body with a homocercal caudal fin and can achieve high swimming speeds (Jordan et al.,1930). The two non-dimensional parameters that characterize the steady inline performance of a carangiform swimmer are the Reynolds number(

*Re*) of the flow and the Strouhal number (

*St*) of the undulatory body motion, which can be defined as follows(Triantafyllou et al., 2000; Lauder and Tytell, 2006):

*L*is the fish length,

*U*is the steady inline swimming speed, ν is the kinematic viscosity of the water,

*A*is the maximum lateral excursion of the tail over a cycle and

*f*is the tail beat frequency.

As carangiform swimmers typically achieve high swimming speeds, their motion is characterized by very high Reynolds numbers, *Re*>10^{4}(Triantafyllou et al., 2000). This range of *Re* values is well within the so-called inertial regime where viscous forces are negligible and inertial forces dominate the dynamics of the motion. As for the Strouhal number *St*, most fishes have been shown to swim near a `universal' optimal value *St*_{opt}=0.3(Triantafyllou et al., 2000). Experimental data with flapping hydrofoils have suggested that fishes prefer this specific value of *St* because the propulsive efficiency is indeed maximized near this optimal *St* value(Triantafyllou and Triantafyllou,1995; Triantafyllou et al.,2000). However, data for Pacific salmon swimming show that this optimal *St* value is not constant, depending strongly on the swimming speed and, by extension, on the *Re* of the flow(Lauder and Tytell, 2006). For instance, at low swimming speeds Pacific salmon swims at *St* much greater than *St*_{opt} (with values as high as *St*=0.6), with the optimal value *St*_{opt}=0.3 being approached only as swimming speed increases(Lauder and Tytell, 2006). Lauder and Tytell comment (Lauder and Tytell, 2006) that such data suggest that fishes at low speeds may either choose for some unclear reason to swim inefficiently or that the *St* alone may not be adequate to explain the intricacies of fish swimming at low speeds. Clearly available data point to a complex relationship between swimming *St* and *Re*, which is far from being understood.

Numerous recent experiments with the state-of-the-art particle image velocimetry (PIV) techniques (Muller et al., 1997; Muller et al.,2000; Nauen and Lauder,2001; Drucker and Lauder,2002; Nauen and Lauder,2002; Tytell and Lauder,2004) have provided a wealth of data in terms of both swimming kinematics and wake flow field. Such experiments cannot, however, clarify the aforementioned relationship between *St* and *Re*, mainly because carrying out controlled experiments in which governing parameters can be systematically varied with live fish is difficult, if not impossible. Another related issue stems from the difficulties in quantifying the swimming efficiency and the locomotive forces from experiments alone. It has been recently shown (for example, Dabiri,2005) that wake velocity or vorticity fields alone (the quantities that are typically measured using PIV) are not sufficient to calculate the locomotive forces, and a pressure-like measurement is also required(Dabiri, 2005). Furthermore,calculating the swimming efficiency based only on wake measurements of a steady swimming fish is not possible as the net momentum in the wake is zero[for a detailed discussion, see Schultz and Webb(Schultz and Webb, 2002)]. Therefore, even in the most recent experiments(Muller et al., 1997; Tytell and Lauder, 2004) the efficiency is calculated based on hydrodynamic models with the kinematic data as input (Lighthill, 1960; Wu, 1960; Lighthill, 1970; Lighthill, 1971; Wu, 1971; Weihs, 1972; Weihs, 1974). Such models,however, are inviscid and should work well for *Re* in the inertial regime but should not be expected to work at lower *Re* when viscous effects play a significant role.

The above discussion underscores how difficult it is for experiments alone to provide conclusive insights into the complex relationship between *Re* and *St* and to explore the energetics of various modes of aquatic swimming under controlled conditions and over a wide range of governing parameters. Such insights can be obtained by combining experimental observations with numerical simulation to design and carry out controlled numerical experiments. However, numerical investigations of fish swimming are relatively scarce, especially when compared with the exploding number of experimental papers dedicated to the same subject. Perhaps the most comprehensive numerical studies are by Wolfgang et al.(Wolfgang et al., 1999) and Zhu et al. (Zhu et al., 2002),who employed an inviscid method to study the wake structures of a straight-swimming giant danio. Their work shed new insights into the vorticity dynamics of the flow, but because of the inviscid assumption of their numerical model their findings are inherently limited in the inertial flow regime. Two-dimensional viscous simulations have been performed in simulated tadpole swimming (Liu and Wassersug,1996) and simulations of a self-propelled eel(Carling and Williams, 1998). Three-dimensional (3D) viscous simulations have been reported for tadpole swimming in a grid with about 4×10^{5} points(Liu and Kawachi, 1999; Liu and Wassersug, 1997). Nevertheless, these simulations were at a fixed *Re*=7200 and could not explore *Re* and *St* effects. Similarly, 3D viscous flow simulations were carried out to investigate the mechanisms of thrust production associated with the flapping aquatic flight of a bird wrasse at a fixed swimming speed and flapping rate, i.e. fixed *Re* and *St*in a computational grid with about 1.5×10^{5} points(Ramamurti et al., 2002). More recently, 3D viscous self-propelled anguilliform swimming was simulated and optimized at several *St* on a mesh with about 3×10^{5}grid nodes (Kern and Koumoutsakos,2006). Numerical simulations for swimming and flying in nature have recently been reviewed by Liu (Liu,2005). These studies have produced important results and shed new light into the hydrodynamics of aquatic swimming. For the most part, however,all these studies focused on simulating a specific aspect or flow regime of aquatic swimming, and as such systematic parametric investigations of the hydrodynamic performance of various modes of aquatic swimming have yet to be reported in the literature.

In the present study we carried out a systematic investigation of carangiform swimming over a range of *Re* and *St*, spanning the transitional and inertial flow regimes. We employed an anatomically realistic model of a mackerel body reconstructed from detailed measurements of an actual fish body. All minor fins were neglected due to lack of detailed kinematical data and only the caudal fin was retained in the model.

The BCF kinematics is prescribed using the experimental data of Videler and Hess (Videler and Hess, 1984). The fish is assumed to be swimming along a straight line at constant speed in a uniform ambient flow. The flow induced by the body undulations is calculated by solving the unsteady 3D Navier–Stokes equations using the sharp-interface, hybrid Cartesian/Immersed-boundary method described elsewhere(Gilmanov and Sotiropoulos,2005; Ge and Sotiropoulos,2007). Calculations are carried out on fine computational meshes(5×10^{6}–10^{7} grid nodes) to ensure grid-independent results and accurate resolution of the viscous region near the fish body. Viscous flow simulations are carried out for two *Re,Re*=300 and *Re*=4000. Inviscid calculations are also carried out representing the flow in the limit of infinite *Re*(*Re*=∞). For all three cases, the *St* is varied systematically, starting from zero (rigid body case), while the swimming speed *U* (i.e. the *Re*) is held constant. Note that in order to be able to vary the *St* while maintaining the swimming speed *U*constant, we simulate the flow induced by a model fish that is attached to and towed by a rigid tether that translates the fish in a stagnant fluid at a given constant velocity *U*. By fixing the speed of the tether *U* we can obtain the desired value of *Re*. The *St* is adjusted by changing the fish tail beat frequency *f*, i.e. by assuming that our virtual swimmer is trained to always undulate its tail at the desired constant frequency. For any given combination of the so-obtained *Re*and *St*, the simulated flow field is used to calculate the force *F* exerted on the fish body by the flow. If *F*≠0 the excess force is absorbed by the hypothetical tether so that the net force acting on the fish is always zero and the constant swimming velocity assumption is satisfied. In such cases, if the hypothetical tether is instantaneously severed, the fish will either accelerate forward or decelerate backward under the action of the excess force *F*. For a given *Re* we vary the *St* until the net mean force acting on the fish is zero, *F*=0. In such a case the hypothetical tether obviously has no effect on the fish,since if it is severed the fish will continue swimming at constant speed *U*. Using this procedure we are able to find for a given *Re*the *St* for which steady, inline swimming is possible. The computed results are analyzed to elucidate several important aspects of carangiform swimming. These include, among others, the ability of carangiform kinematics to produce thrust as a function of *Re*, the swimming efficiency and propulsive power requirements in the transitional and inertial regimes, and the 3D structure of the wake as a function of *Re* and *St*.

The paper is organized as follows. First, we briefly describe the numerical method and present the details of the fish model and prescribed kinematics. We then discuss the numerical experiments and the results in terms of drag increase/reduction, swimming efficiency, and the 3D vortical structures in the wake. Finally we summarize our findings, present the conclusions of this work,and outline areas for future research.

## MATERIALS AND METHODS

### The governing equations and boundary conditions

*U*and the fish length

*L*, the governing equations read as follows:

*u*

_{i}are the non-dimensional Cartesian velocity components of the fluid,

*p*is the pressure divided by the density,and

*D*/

*Dt*is the material derivative defined as(

*D*/

*Dt*)(

^{.})=(∂/∂

*t*)(

^{.})+

*u*

_{j}(∂/∂

*x*

_{j})(

^{.}). The inviscid (Euler) equations, which are also solved in this work, are recovered from Eqn 3 by letting

*Re*→∞. We are interested in solving these equations in a domain containing an arbitrarily complex 3D flexible body moving with prescribed kinematics. Therefore, the governing Eqn 3 needs to be supplemented with appropriate boundary conditions at the outer boundary of the flow domain,which could be either occupied by ambient fluid or enclosed by a solid surface, as well as the inner moving immersed boundaries.

*t*). Γ(

*t*) is discretized with K material points, which lie on it at all times and can be tracked with their global Lagrangian position vectors

**r**

*(*

^{k}*t*):

*k*

^{th}material point on Γ(0). Here, the motion ofΓ(

*t*) is prescribed with known velocities

**U**

*(*

^{k}*t*) – the prescribed swimming kinematics. Therefore, the shape of Γ(

*t*) at time

*t*can then be obtained by solving the advection equation for all material points on the surface (for

*k*=1,K):

*t*) known for time

*t*, boundary conditions for the Eulerian fluid velocity vector

**u**(

**r**,

*t*) must be prescribed at all points ofΓ(

*t*). For viscous flow, the no-slip and no-flux boundary conditions need to be satisfied as follows:

*u*

_{n}is the fluid velocity normal to the body,

**u**is the fluid velocity vector tangential to the body and

_{t}**n**is the normal vector to the surface.

### The numerical method

The flexible fish body is handled as a sharp interface immersed in the background Cartesian grid using the hybrid Cartesian/immersed-boundary (HCIB)method, which has been described in detail(Gilmanov and Sotiropoulos,2005) and so only a very brief description of the technique is given herein. The method employs an unstructured, triangular mesh to discretize and track the position of a fish body. Boundary conditions for the velocity field at the Cartesian grid nodes that are exterior to but in the immediate vicinity of the immersed boundary (IB nodes) are reconstructed by linear or quadratic interpolating along the local normal to the boundary. No explicit boundary conditions are required for the pressure field at the IB nodes due to the hybrid staggered/non-staggered mesh formulation(Gilmanov and Sotiropoulos,2005). The HCIB reconstruction method has been shown to be second-order accurate on Cartesian grids with moving immersed boundaries(Gilmanov and Sotiropoulos 2005). The IB nodes at each time step are identified using an efficient ray-tracing algorithm (Borazjani et al., 2008)

The method has been validated extensively(Gilmanov and Sotiropoulos,2005) for flows with moving boundaries and has also been applied to simulate fish-like swimming using the mackerel model we employ in this work. More specifically, inviscid simulations for various slip ratios *U*/*V* (where *V* is the body undulatory wave phase velocity) showed that the simulated wake structures for various slip ratios are in good agreement with experimental observations: for *U*/*V*<1 a reverse Karman street was obtained, in agreement with the observations (Gilmanov and Sotiropoulos, 2005). Gilmanov and Sotiropoulos also carried out viscous flow simulations for *Re*=3000 and showed that on a sufficiently refined mesh (3×10^{6} grid nodes) the method can simulate a thrust-producing wake with a reverse Karman street (for details,see Gilmanov and Sotiropoulos,2005).

In Gilmanov and Sotiropoulos' study(Gilmanov and Sotiropoulos,2005), the governing equations were solved with an explicit dual-time stepping artificial compressibility method. To further enhance the efficiency of the numerical method, which is a crucial prerequisite for carrying out highly resolved viscous flow simulations on grids with 10^{7} grid nodes, we modified the flow solver by incorporating a recently developed fractional step method(Ge and Sotiropoulos, 2007). The Poisson equation is solved with FGMRES(Saad, 2003) and multigrid as preconditioner using parallel libraries of PETSc(Satish Balay et al., 2001). For more details the reader is referred elsewhere(Ge and Sotiropoulos, 2007; Borazjani et al., 2008).

### Fish body kinematics and non-dimensional parameters

We employ a fish body shape closely modeled after the actual anatomy of a mackerel. The emphasis in our work is on the body/caudal fin mode of aquatic swimming and for that only the caudal fin is retained in the model fish. The model is meshed with triangular elements as needed by the HCIB method(Fig. 1).

*L*):

*z*is the axial (flow) direction measured along the fish axis from the tip of the fish head;

*h*(

*z,t*) is the lateral excursion at time

*t*;

*a*(

*z*) is the first Fourier coefficient (Videler and Hess,1984) defining the amplitude envelope of lateral motion as a function of

*z*;

*k*is the wave number of the body undulations that corresponds to a wavelength λ; and ω is the angular frequency.

The four important non-dimensional similarity parameters in fishlike swimming are: (1) the Reynolds number *Re*, based on the fish length *L*, the swimming speed *U* and the fluid kinematic viscosityν: *Re*=*LU*/ν; (2) the Strouhal number *St*, based on the maximum lateral excursion of the tail *A*=2*h*_{max}, and the tail beat frequency *f*: *St*=2*fh*_{max}/*U*; (3) the non-dimensional wavelength λ/*L*; and (4) the non-dimensional amplitude envelope *a*(*z/L*)/*L.* Sometimes the so-called slip velocity or slip ratio, defined as slip=*U*/*V=U*/(ω/*k*), is used instead of non-dimensional wavelength. Using either parameter is correct. However, the slip velocity changes if the tail beat frequency is changed,while the wavelength and the tail beat frequency are independent.

*L*and the amplitude envelope

*a*(

*z*) parameters, named shape parameters hereafter, are specified such that the fish body motion is similar to the typical carangiform swimmers' body motion. The amplitude envelope

*a*(

*z*) can be well approximated by a quadratic curve of the form:

*a*

_{0}=0.02,

*a*

_{1}=–0.08 and

*a*

_{2}=0.16 to match the experimental curve of Videler and Hess(Videler and Hess, 1984). With this amplitude envelope the maximum displacement at the tail will be

*a*

_{max}=0.1, which gives

*h*

_{max}=0.1

*L*, which is similar to that of an average mackerel. The wave number

*k*in all the simulations is based on the non-dimensional wavelength λ/

*L=*95%, which is in the range of 89–110% observed in most carangiform swimmers(Videler and Wardle,1991).

In all the simulations, as explained above, the shape parameters are kept constant, similar to the carangiform swimmers, while the *Re* and *St* are varied. The non-dimensional angular frequency used in Eqn 9 is calculated based on the *St* asω=2π*fL/U*=2π*St*/2*a*_{max}.

The above non-dimensional angular frequency ω is used along with the non-dimensional time *tU/L* in Eqn 9. Fig. 2 shows the midlines of the fish calculated for one tailbeat cycle using Eqn 9 with the carangiform shape parameters, and the quadratic amplitude envelope calculated by Eqn 10, which was fitted through the experimental curve of Videler and Hess(Videler and Hess, 1984).

In this work we have assumed that the fish tail motion (1) exactly follows the body's backward traveling wave; and (2) is completely symmetric. The first assumption is supported by the work of Gibb et al., who studied in detail the tail kinematics of a chub mackerel and showed that the traveling wave does indeed pass through the tail and as such the tail does not move like a rigid plate attached to the body (Gibb et al.,1999). Gibb et al. also showed, however, that contrary to our second assumption, the tail of a chub mackerel undergoes a dorso-ventral(top–bottom) asymmetry during each tail beat(Gibb et al., 1999). This effect has not been considered in the present study but could be examined in a future study using the detailed kinematics reported(Gibb et al., 1999).

### Computational grid and other details

As already explained in the Introduction, in all our simulations it is assumed that the fish is attached to a rigid tether that tows the fish at constant velocity *U*. Therefore, all the equations are solved in the inertial frame moving with constant velocity *U* attached to the fish. The computational domain is a 2*L*×*L*×7*L*cuboid similar to that used by Gilmanov and Sotiropoulos(Gilmanov and Sotiropoulos,2005), but discretized with a much finer grid including 5.5 million grid nodes. The domain width 2*L* and height *L* are ten times the fish width 0.2*L* and height 0.1*L*, respectively. A uniform mesh with constant spacing *h*=0.008*L* is used to discretize a smaller cuboid enclosing the fish. The mesh is stretched from the faces of this smaller cuboid to the boundaries of the computational domain using a hyperbolic tangent stretching function. The fish is placed 1.5*L* from the inlet plane in the axial direction and centered in the transverse and the vertical directions. The boundary conditions on the domain outer boundaries are uniform flow at the inlet, slip walls on the side boundaries and convective boundary conditions at the outlet.

To test the sensitivity of the computed solutions to the size of the computational domain, a set of simulations was also carried out for a longer domain with dimensions 2*L*×*L*×11*L*. Similar to the shorter (7*L* long) domain, a uniform mesh with constant spacing *h*=0.008*L* is used to discretize a smaller cuboid enclosing the fish and the mesh is stretched from the faces of this smaller cuboid to the boundaries of the computational domain using a hyperbolic tangent stretching function. This results in a grid with 11.4 million nodes. Simulations were carried out for *Re*=4000 and ∞ (inviscid) for various *St* and the computed results were found to be in excellent agreement with those obtained in the shorter domain both in terms of instantaneous and time averaged forces and flow structures. Based on this study it was concluded that the 7*L* long domain is sufficient for carrying the parametric studies reported in the remainder of this paper.

A grid refinement study is also carried out for the *Re*=4000 case using a series of successively finer meshes. The results of this study are reported in the Appendix, where more details on the validation and verification of the numerical method are provided. Here it suffices to mention that based on this grid sensitivity study we concluded that the 5.5 million node grid is adequate for obtaining grid insensitive results.

The tail beat period τ is divided in 120 time steps, i.e.Δ *t*=τ/120. A time refinement study withΔ *t*=τ/1000 was also carried out and showed no appreciable differences in the computed flow patterns and the time averaged net force on the fish body.

### Calculation of swimming forces and efficiency

*F*force is zero for steady inline swimming when the thrust force is exactly equal to the hydrodynamic drag force. It is useful, however, to define a Froude propulsive efficiency based on the thrust force for constant speed inline swimming as follows (see Tytell and Lauder, 2004):

*T̄*is the average over the swimming cycle thrust force,

*U*is the steady swimming speed, and

*P̄*

_{L}is the average over the swimming cycle power loss due to lateral undulations. As previously pointed out (Tytell and Lauder, 2004)the efficiency defined by Eqn 11expresses the percentage of the total power that is used to produce thrust. The main problem with Eqn 11,however, stems from the fact that the thrust force

*T*cannot be measured directly in experiments (Schultz and Webb, 2002). However, the thrust and lateral power loss can be calculated either using a mathematical model, such as Lighthill's elongated body theory (EBT) for steady swimming(Lighthill, 1971), or directly from the results of 3D numerical simulations. In the remainder of this section we discuss these two approaches for calculating swimming efficiency.

*Re*. In addition, the EBT efficiency can be compared with the efficiency obtained

*via*direct calculation from numerical simulations (see below) to evaluate the range of validity of the theory. The Froude efficiency based on EBT for steady swimming is given as follows:

*U*/

*V*is the slip velocity defined as the ratio of the swimming speed

*U*to the speed

*V*of the backward undulatory body wave. An improved EBT efficiency formula(Cheng and Blickhan, 1994)(denoted herein as EBT-2) takes into account the slope of the fish tail where all the mean quantities are computed:

*h*(

*L*) is the undulation amplitude and

*h*′(

*L*) is its derivative (slope) at the tail.

The Froude propulsive efficiency given by Eqn 11 can also be calculated directly from the results of 3D computations; we refer to this approach as CFD. To accomplish this, however, we first need to define and develop an approach for calculating the thrust and drag forces. Note that for fish-like swimming, such a definition is not straightforward since the propulsor in this case is the fish body itself, which produces thrust while producing drag.

*x*

_{3}direction. The component of the instantaneous hydrodynamic force along the

*x*

_{3}direction (which for simplicity will be denoted as

*F*) can be readily computed by integrating the pressure and viscous forces acting on the body as follows(where repeated indices imply summation):

*n*is the

_{i}*i*th component of the unit normal vector on d

*A*and τ

*is the viscous stress tensor. Depending on whether*

_{ij}*F*(

*t*) is negative or positive,it could contribute to either drag

*D*(

*t*) or thrust

*T*(

*t*). To separate the two contributions we propose to decompose the instantaneous force as follows:

*F*(

*t*) is of thrust or drag type, respectively. By construction the decomposition given by Eqn 16 and Eqn 17conserves identically the net force:

*ḣ*is the time derivative of the lateral displacement (

*i*=2 direction), i.e. the velocity of the lateral undulations.

It is important to note that the Froude efficiency equation(Eqn 11) can only be applied under inline, constant-speed swimming when the thrust force is balanced exactly by the drag force and the net force acting on the fish body is zero. If this equilibrium condition is violated, the fish will either accelerate or decelerate, the velocity *U* will no longer be constant and Eqn 11 is not meaningful. In our subsequent simulations with the previously described tethered fish model, we will only apply Eqn 11 to compute the propulsive efficiency at the critical Strouhal number *St*^{*} for which the net force *F* acting on the fish body is zero and constant-speed inline swimming is possible.

## RESULTS AND DISCUSSION

### Strouhal number and Reynolds number effects

To systematically quantify the effect of governing parameters on the energetics and associated flow patterns of carangiform swimming we carry out simulations for three *Re* and various *St*. Viscous flow simulations are carried out for *Re*=300 and 4000 and inviscid simulations are carried out to simulate the flow in the limit of *Re*=∞. For *Re*=300 and 4000, the *St* is varied systematically from zero (rigid body case) in increments of 0.1 until the mean net force on the fish body becomes greater than zero (see below for details). For *Re*=∞ simulations are carried out over a narrower range of *St* centered around the value at which the net force on the fish crosses zero.

*C*

_{F}as a function of

*St*for

*Re=*4000. The force coefficient is defined as follows:

*F*(

*t*) is the instantaneous net hydrodynamic force given by Eqn 15, ρ is the density of the fluid,

*U*is the swimming speed, and

*L*is the length of the fish body. Recall that in our simulation the fish cannot move, and thus the net hydrodynamic force is absorbed by the hypothetical tether that holds the fish in place. In other words, the force shown in Fig. 3 is the net force that would be available to accelerate the fish either forward or backward (depending on its sign) at the instant when the hypothetical numerical tether is removed. Given the sign convention we introduced in the previous section,

*C*

_{F}>0 when

*T*>

*D*, i.e. when the thrust force exceeds the drag force and the net force on the body is in the direction of the fish motion. To facilitate our discussion we shall refer to this situation as the net force being of thrust type. Similarly the situation with

*C*

_{F}<0 will be referred to as the net force being of drag type. Such notation is used herein to characterize the direction of the net force and should not be confused with the terms thrust or drag force, which refer only to the thrust or drag portions of the instantaneous net force (see Eqn 15, 16, 17). The values of

*C*

_{F}in Fig. 3 and in all subsequently presented figures have been scaled with the force coefficient calculated for the rigid body fish (

*St*=0) at the same

*Re*. The line corresponding to the resulting rigid body force

*C*

_{F}=–1 is marked in Fig. 3 to help gauge the level of the net force for each

*St*relative to the rigid body drag. The most important findings that follow from Fig. 3 can be summarized as follows:

For all simulated

*St*, the force coefficient in each cycle shows two peaks corresponding to forward and backward tail strokes. This finding is in agreement with experimental observations(Hess and Videler, 1984).As the

*St*increases from zero, the net force remains of drag type(*C*_{F}<0) throughout the entire swimming cycle up to a threshold*St*(*St*≈0.5) at which the first excursions into the thrust-type regime (*C*_{F}>0) are observed.Further increase of the

*St*above the*St*≈0.5 threshold leads to longer and larger amplitude excursions into the thrust-type regime,ultimately yielding a positive mean net force.For small

*St*, ⩽0.3, the undulations of the body cause a net force of drag type with magnitude greater than the drag force of the rigid fish at the same*Re*. That is, low*St*body undulations cause the magnitude of the drag-type net force to increase over that of the rigid body. For higher*St*(*St*>0.3), the undulations of the body cause a net force that is also of drag type but of lower magnitude than the corresponding rigid body net force.

The specific *St* at which the net force changes sign is not universal but depends on the *Re* of the flow. A plot similar to that shown in Fig. 3 for *Re*=300, for instance, exhibits essentially all qualitative trends as those observed for *Re*=4000, but the net force sign transition occurs at different *St*. To illustrate the dependence of the net force variation on the *St* as a function of the *Re*, we plot in Fig. 4 the variation of the mean net force coefficient *C̄*_{F} (averaged over several swimming cycles) for all three simulated *Re*. As before, the values of *C̄*_{F} for each *Re* are scaled by the corresponding value for the rigid body at the same *Re*,i.e. the line *C̄*_{F}=–1 marks the rigid body case (*St*=0). For the *Re*=∞ case, *C̄*_{F} is scaled with the value for the *Re*=4000 since for inviscid, irrotational flow the net force on the body is zero. It is observed from Fig. 4 that for *Re*=300 and 4000 at low *St* the mean net force is of drag type and its magnitude initially increases relative to that of the rigid body. As the *St* increases, however, the mean net force, while remaining of drag type, is gradually diminishing in magnitude and ultimately its magnitude becomes smaller than the rigid body force. The *St* at which this transition occurs appears to be the same for both *Re*(*St*∼0.25). As the *St* is further increased, above a certain *St* threshold (denoted by *St*^{*}) the force becomes positive in the mean, which indicates the transition from a mean net force of drag type to a mean net force of thrust type. This general trend is observed for all three simulated *Re*. The point when the *C̄*_{F} curve crosses zero is the point when the mean drag force is balanced exactly by the mean thrust force. As we have already discussed, at this point the fish will swim at constant velocity if the hypothetical tether in our simulations is severed. Therefore, the *St* at which *C̄*_{F} crosses zero is the *St* at which constant speed, inline swimming is possible for the given *Re*; we shall denote this *St* as *St*^{*}.

A striking finding from Fig. 4 is that *St*^{*} is a decreasing function of *Re*; *St*^{*}=1.08, 0.6 and 0.26 for *Re*=300,4000 and ∞, respectively. Moreover, *St*^{*} approaches the range of *St* at which most carangiform swimmers swim in nature(*St* ∼0.2–0.35) (see Fish and Lauder, 2006) only in the limit of *Re*=∞. Recall that this is also the range of *St* at which optimal thrust production has been reported in flapping foil experiments (Triantafyllou et al.,1991), which led to the hypothesis that fishes select this range of *St* to optimize their propulsive efficiency(Triantafyllou and Triantafyllou,1995; Triantafyllou et al.,2000). Nonetheless, it is clear from Fig. 4 that for each *Re* there is a unique *St, St*^{*}, at which steady inline swimming is possible. Therefore, our results suggest that in addition to efficiency considerations, which we will further elaborate on below, for a given *Re* carangiform swimmers select the *St* at which they will undulate their body because this is the only *St* at which they can produce enough thrust to cancel the drag they generate and swim steadily. This finding also suggests a possible explanation for the data for Pacific salmon swimming reported (Lauder and Tytell, 2006). As discussed already in the Introduction, Lauder and Tytell (Lauder and Tytell,2006) presented data showing that the swimming *St*increases with decreasing swimming speed (i.e. decreasing *Re*). Based on these data they pondered whether fish for some unknown reason choose deliberately to swim inefficiently at low speeds and also wondered whether the *St* number is the appropriate parameter to quantify the intricacies of aquatic swimming across a wide range of swimming speeds. The results shown in Fig. 4 suggest that at lower swimming speeds Pacific salmon swims at higher *St* simply because *St*^{*} is a decreasing function of *Re*.

### Swimming efficiency

As we discussed in a previous section, the Froude efficiency given by Eqn 11 is meaningful to calculate for a given *Re* only at *St*^{*}, when the assumption of constant swimming speed is valid. In Table 1 the Froude efficiency is given for different *Re* at the corresponding value of *St*^{*} using the EBT (Eqn 12), EBT-2 (Eqn 13)and direct (CFD) calculation (Eqn 11). It is clearly evident from this table that regardless of the calculation method, the swimming efficiency is an increasing function of the *Re*. That is, even though carangiform kinematics can achieve constant-speed swimming, i.e. self-propulsion, at all simulated *Re*,this mode of swimming is very inefficient at low *Re*. Higher propulsive efficiency can only be achieved in the limit of *Re*=∞, which is the range of *Re* at which typical carangiform swimmers swim in nature. To the best of our knowledge, this is the first that time that the effects of scale (*Re*) on propulsive efficiency are so clearly demonstrated. Note that the fact that the efficiency of the carangiform kinematics is higher at higher *Re* is entirely consistent with the fact that *St*^{*} is decreasing with *Re*. Body undulations at higher *St* generally imply faster lateral undulations, i.e. higher lateral velocities, which result in higher lateral power loss and lower efficiency.

. | Froude efficiency (%) . | . | . | ||
---|---|---|---|---|---|

Re
. | EBT . | EBT-2 . | CFD . | ||

300 | 59.3 | 59.5 | 18.86 | ||

4000 | 67.25 | 66.67 | 22.95 | ||

∞ | 92 | 89.47 | 47.55 |

. | Froude efficiency (%) . | . | . | ||
---|---|---|---|---|---|

Re
. | EBT . | EBT-2 . | CFD . | ||

300 | 59.3 | 59.5 | 18.86 | ||

4000 | 67.25 | 66.67 | 22.95 | ||

∞ | 92 | 89.47 | 47.55 |

*Re,* Reynolds number; EBT, elongated body theory; EBT-2, EBT efficiency formula (Cheng and Blickhan,1994); CFD, computational fluid dynamics

It is important to also comment on the consistency of the results obtained from the three methods used to calculate the Froude efficiency in Table 1. It is apparent that qualitatively the three methods yield identical results, all three predicting the same trend of increasing efficiency with increasing *Re*. As expected, the EBT and EBT-2 methods yield similar efficiency values. The slope correction in the EBT-2 approach becomes more important at higher *Re*(i.e. lower *St*^{*}), which accounts for the increasing discrepancy between the EBT and EBT-2 with *Re*. The efficiency values from the two EBT methods differ significantly from the values calculated by the CFD method with the discrepancy becoming greater as the *Re* is decreased. Such discrepancy, however, is entirely consistent with the fact that the EBT theory is an inviscid, slender-body theory, it does not consider 3D effects, and is not at all applicable at *Re* sufficiently low for the viscous forces to play a significant role.

The efficiency values we obtain in this work using CFD are somewhat lower than those reported in the literature in previous numerical studies. For example, 45% efficiency was reported for a tadpole swimming at *Re*=7200 (Liu and Kawachi,1999) while efficiencies ranging from 62% to 72% have been reported for inviscid simulations of tuna swimming(Zhu et al., 2002). It is important to point out, however, the aforementioned studies simulated tethered swimmers at conditions that did not correspond to the self-propulsion state,i.e. in these simulations *St*≠*St*^{*} and the net force on the body was not zero (*T*≠*D*). This is a significant difference between our work and previous studies, which could very well account for the differences in the levels of the calculated swimming efficiencies.

### Power requirements of undulatory locomotion

In this section we employ the results of our simulations to calculate the power requirements of undulatory swimming and compare them with the power requirement of towing the rigid fish at the same speed based.

At *St*^{*} the average axial power is zero since the average axial force per cycle is zero. Therefore, the total power required at *St*^{*} is only the side power calculated by Eqn 19. The power requirement for the rigid fish to be towed at velocity *U* is simply the drag force multiplied by the velocity *U*. The power requirement has been calculated and non-dimensionalized by the factorρ *U*^{3}*L*^{2} and the values are reported in Table 2.

Re
. | P_{rigid}
. | P_{St*}
. |
---|---|---|

300 | 2.6165×10^{–2} | 1.4157×10^{–1} |

4000 | 4.8921×10^{–3} | 1.3671×10^{–2} |

Inviscid | 0 | 3.9598×10^{–4} |

Re
. | P_{rigid}
. | P_{St*}
. |
---|---|---|

300 | 2.6165×10^{–2} | 1.4157×10^{–1} |

4000 | 4.8921×10^{–3} | 1.3671×10^{–2} |

Inviscid | 0 | 3.9598×10^{–4} |

*Re,* Reynolds number; *P,* power; *St*^{*},critical Strouhal number

The results in Table 2clearly show that the power requirement of undulatory swimming decreases as *Re* increases. In conjunction with the conclusion reached in the previous section, this finding shows that at very high *Re* carangiform kinematics not only become very efficient but also require less power for propulsion. Given the previously discussed ambiguities in the definition of a meaningful and objective measure of swimming efficiency, our results reinforce the recommendation (Schultz and Webb,2002) that the power requirement might be a better measure for quantifying the efficiency of fish swimming.

Table 2 also shows that for a given *Re* the power requirement of the undulatory swimming is higher than towing the rigid fish. All the kinematic and computational models to date have shown the same trend (for a review, see Schultz and Webb, 2002). However, this finding is in contrast to other results(Barrett et al., 1999), which showed through experimental measurements with a robotic fish that the power required for the tethered fish to move at constant speed *U* with undulatory body motion is less than that for the rigid body. It is important to point out, however, that whether body undulations increase or decrease, the power required for swimming at all *Re* cannot be deduced with certainty from the results of Table 2 or other kinematic and computational models. This is because previous models were either inviscid or did not consider *Re* in the range bridging the transitional and inertial regimes. It still might be possible that at some sufficiently high *Re*(4000<*Re*<∞) the undulatory swimming requires less power than the rigid fish. To demonstrate this possibility, consider the following argument based on our results. As the *Re* is increased, the *St*^{*} and the power requirement for undulatory swimming are decreased. The rigid body drag coefficient, on the other hand (which multiplied by the velocity gives the power required for towing), tends to asymptote toward a constant value, similar to other bluff bodies(Panton, 1996). Clearly this important point cannot be conclusively resolved by our work and simulations at much higher *Re* will be required for definitive conclusions.

### Is undulatory locomotion drag-reducing or drag-increasing?

Barrett et al. (Barrett et al.,1999) concluded that the undulatory motion is drag-reducing since the upper bound drag estimate was found to be less than that of the corresponding rigid body drag. However, Lighthill's analytical results(Lighthill, 1971), which were verified by simulations (Liu and Kawachi,1999) and experiments (Anderson et al., 2001), indicate that the friction drag increases by the undulatory motion. In addition, Fish et al.(Fish et al., 1988; Fish, 1993) suggested that the pressure component of the drag increases by the undulatory motion in dolphins and seals, thus concluding that the undulatory motion is drag-increasing. In summary, the issue of whether undulatory motion is increasing or decreasing the drag force is still a subject of debate in the literature, as often the reported results appear to contradict each other. It is evident from the above discussion that the main reason for these contradictory conclusions is the inherent difficulty of experimental studies and simplified theoretical models to calculate simultaneously both the total force and its two components. Obviously our numerical simulations do not suffer from this limitation and can be interrogated in detail to elucidate the issue of the impact of undulatory body motion on the hydrodynamic drag force. Using the simulated flow fields we compute the total drag force *D* (Eqn 17) and its form and friction drag components given by *D*_{p} and *D*_{v} in Eqn 17 as functions of *St* for *Re*=300 and 4000. The results are normalized by the rigid body total drag and plotted in Fig. 5. This figure reveals several important trends.

For both *Re* the total drag force initially increases above that of the rigid body drag up to *St*=0.1, with the overall increase level being higher for the higher *Re*. At higher *St* the drag force starts to decrease and at *St*≈0.25 it decreases below the rigid body drag for both *Re*. Beyond that point, however, a distinctly different behavior is observed for the two *Re*. For *Re*=300 the drag starts increasing again above the rigid body threshold, while for *Re*=4000 the drag is reduced monotonically, asymptoting toward a constant value of approximately 75% of the rigid body drag at *St*≈0.6.

The friction drag force increases monotonically with *St* for both *Re* while the form drag initially increases and then decreases asymptoting toward zero at *St*>0.6. As one would anticipate the friction drag is the major contributor to the total drag force at the low *Re* (*Re*=300) and is responsible for the monotonic increase of the total drag force for *St*>0.3. For *Re*=4000 the friction drag is higher than the form drag but varies only mildly with *St*,increasing from 0.66 for the rigid body to an asymptotic limit of 0.75 for *St*>0.5. Consequently the variation of the total drag for this case is dominated by the non-monotone variation of the form drag, which as mentioned above initially increases up to *St*=0.1 and then asymptotes to zero for *St*>0.6.

The above results provide new insights that help reconcile the previous,often conflicting, reports about the effect of body undulations on the drag force. Clearly for both *Re* the friction drag is the major contributor to the total drag and it is increased by the undulatory motion, a trend that is consistent with previous simulations(Liu and Kawachi, 1999) and experiments (Anderson et al.,2001). The form drag on the other hand is initially increased by body undulations but then decreases toward zero as the *St* increases. As the *Re* is increased, however, the importance of viscous stresses diminishes and the friction drag tends to become fairly insensitive to the *St* and the variation of the total drag mimics essentially that of the pressure drag. Even though in our simulations, due to limitations in computational resources, we were not able to reach the range of *Re* at which the experiments of Barrett et al.(Barrett et al., 1999) were carried out, our results suggest that the total drag reduction observed in these experiments is mainly due to the ability of the undulatory body motion to drastically reduce the form drag. Note that the simplified inviscid hydrodynamic models (Lighthill,1971; Fish et al.,1988; Fish, 1993)could not explain the drag reduction observed in experiments since they did not considered form drag and the assumption was that the thrust overcomes the drag only due to friction.

As we conclude this section it is appropriate to discuss the physical mechanisms that lead to the observed reduction in form drag. It has been hypothesized (Triantafyllou et al.,2000) that in undulatory swimming the travelling body wave contributes to a decreased drag force by eliminating separation and suppressing turbulence. This hypothesis was supported by early experiments(Taneda and Tomonari, 1974),which visualized the flow past a waving flat plate and showed that when the wave phase velocity *V* is smaller than the flow velocity *U*,the boundary layer separates at the back of the wave crest, while when *V*>*U* the boundary layer does not separate. A more recent study (Shen et al., 2003)carried out direct numerical simulation of flow past a waving plate and confirmed the earlier findings (Taneda and Tomonari, 1974) by showing that for *V*>*U*separation is indeed eliminated and drag is reduced, relative to the stationary wavy wall, monotonically with increasing *V*. Shen et al. also emphasized the relevance of their waving flat plate results in understanding drag-reduction mechanisms in fish-like swimming(Shen et al., 2003).

Our computed results also show that the ratio of the undulatory wave phase velocity *V* to the swimming speed *U* is indeed a critical parameter insofar as drag reduction is concerned. In our simulations *V* exceeds *U* when *St*>0.22 and it is evident from Figs 4 and 5 that for the viscous flow simulations the drag force is first reduced below that of the rigid body for *St*>0.25, i.e. when the condition *V*>*U* has been satisfied. In Fig. 6 we show instantaneous streamlines and pressure contours at the mid-plane of the virtual swimmer in the frame of reference moving with the undulatory wave phase velocity *V* for *Re*=300 and *St*=0 (rigid body),0.1 (*U*/*V*=2.11) and 0.3 (*U*/*V*=0.7),respectively. Note that the moving frame of reference is selected because in the case of a swimming fish flow separation occurs relative to the undulating body and can only be visualized clearly in the frame of reference that moves with the body wave velocity *V* [see Shen et al.(Shen et al., 2003) for a detailed discussion of this issue]. As seen from Fig. 6A the flow around the rigid body (*St*=0) does not separate, due to the streamlined body shape. As the body begins to undulate, and as long as *V* is less than *U*, the flow separates at the posterior of the body (see results in Fig. 6B for *St*=0.1)due to the fact that the undulatory body wave is such that it acts to retard the near-wall flow relative to the free stream. The onset of separation explains the initial increase of the form drag force relative to the rigid body drag observed in Figs 4and 5. At *St*sufficiently high for the condition *V*>*U* to be satisfied(*St*>0.22 in our case), however, separation is eliminated (see Fig. 6C for *St*=0.3)and the drag force is reduced below that of the rigid body drag at the same *Re*. Under these conditions (i.e. *V*>*U*), the motion of the undulating fish body is piston-like and acts to accelerate the slower moving ambient fluid, thus creating a positive (stagnation) pressure region at the posterior portion of the fish body, which reduces the form drag;this is clearly evident in the pressure contours shown in Fig. 6C. It is of course important to note that the reduction in the form drag at higher *St* is not for free. The fish has to beat its tail faster to achieve drag reduction at higher *St* and thus needs to expend more power to accomplish this. As discussed in the previous section the lateral power cost of such drag reduction was found to be higher than the inline power gained by the undulatory motion in our simulations.

### Three-dimensional wake structure

The wake of carangiform swimmers has been studied extensively in the laboratory using particle image velocimetry (PIV), which can provide the velocity field in several 2D planes(Muller et al., 1997; Wolfgang et al., 1999; Nauen and Lauder, 2002). These experiments showed the vortices in the wake of free swimming carangiform fishes organize in a single row such that a jet flow is formed between the vortices, which has been dubbed a reverse Karman street(Rosen, 1959).

In our simulations we also find a reverse Karman street wake consisting of a single row of vortices for the self-propelled inviscid flow case(*St*^{*}=0.26), which is the case that corresponds more closely (both in terms of *Re* and *St*) to the available in the literature experiments with live carangiform swimmers. The simulated near-wake velocity and vorticity fields in the horizontal and vertical planes for this case are shown in Fig. 7. The flow patterns shown in this figure are very similar to those obtained experimentally (Nauen and Lauder,2002) using PIV on the horizontal and vertical planes near the caudal fin of a swimming mackerel [see figs 3 and 4 in Nauen and Lauder(Nauen and Lauder, 2002),corresponding to Fig. 7A and Bin the present study].

Our simulations have also revealed different wake patterns depending on the *Re* and *St* numbers, as shown in Fig. 8, which depicts three such representative wake patterns. To facilitate the classification of the various wake patterns observed in our simulations let us introduce the following wake characterization convention. (1) Depending on the direction of the common flow between the wake vortices, the wake can be characterized as being of thrust (reverse Karman street) or drag (regular Karman street) type;(2) depending on the layout of the vortical structures the wake can be characterized as consisting of single or double row vortices. A single row wake can be either of drag or thrust type, as shown in Fig. 8A (*Re*=4000 and *St*=0.2) and Fig. 8B(*Re*=∞ and *St*=0.26), respectively. The main characteristic of this wake pattern is that it remains confined within a relatively narrow parallel strip that is centered on the axis of the fish body and consists of Karman-street like vortices. A double row wake that is of thrust type is shown in Fig. 8C(*Re*=4000, *St*=0.7). This wake pattern is distinctly different than the single row wake as it is characterized by the lateral divergence and spreading of the vortices away from the body in a wedge-like arrangement.

Our computed results show that for fixed *Re* both the single and double row wake structures can emerge, depending on the *St* number. Typically at low *St* the single row wake structure is observed (see Fig. 8A,B), while at high *St* the wake splits laterally and the double row pattern emerges(Fig. 8C). The dependence of the wake structure on the *St* is to be expected since by definition(see Eqn 2) the *St* can be viewed as the ratio of the average lateral tail velocity to the axial swimming velocity. Therefore, at high *St* the vortices shed by the tail tend to have a larger lateral velocity component, which advects them away from the centerline causing them to spread in the lateral direction.

The exact value of *St* at which the transition from the single to the double row wake structure occurs depends on the *Re*. Due to limitations in the computational resources at our disposal we did not attempt to precisely calculate the wake transition *St* as a function of *Re*. Our results, however, do suggest that for *Re*=300 the wake transition occurs within 0.3<*St*<0.6 while for *Re*=4000 and ∞ it occurs at somewhat lower *St* in the range 0.3<*St*<0.5.

To illustrate the effect of *Re* on the wake structure, we plot in Fig. 9 the instantaneous vorticity field and streamlines on the mid plane for all three *Re* for *St*=0.3. As one would anticipate, for lower *Re* the thickness of the viscous regions around the body and overall width of the wake become greater as the diffusive effects of the viscous forces begin to dominate. For the *Re*=300 case, the wake is of drag type with a single row of vortices. At *Re*=4000 the wake is still of drag type but it is clearly more disorganized and complex than the *Re*=300 case. The wake pattern is of single row type but its emerging complexity signals its upcoming (for *St*>0.3) transition from the single to double row structure. For the inviscid case the wake also consists of single row vortices but it is now of thrust type. In comparison with the inviscid wake shown in Fig. 8B for *St*=0.26,the wake for the *St*=0.3 case has become more complex. The vortices have intensified, adjacent layers of positive and negative vorticity are observed in the wake, and two cores of high vorticity emerge within the primary wake vortices especially at some distance downstream of the tail. Once again, the emergent complexity of the wake is also suggestive of the transition to the double row structure that will occur at somewhat higher *St*.

To visualize the 3D structure of the wake we show in Fig. 10, Fig. 11 and Fig. 12 instantaneous iso-surfaces of the *q*-criterion(Hunt et al., 1988) for *Re*=300, 4000 and ∞, respectively. The quantity *q* is defined as *q*=½(∥Ω∥^{2}–∥*S*∥^{2}),where *S* and Ω denote the symmetric and antisymmetic parts of the velocity gradient, respectively, and ∥** ^{.}**∥ is the Euclidean matrix norm. According to Hunt et al.(Hunt et al., 1988), regions where

*q*>0, i.e. regions where the rotation rate dominates the strain rate, are occupied by vortical structures. For each

*Re*we show two

*St*corresponding to the single and double row vortex patterns. The 3D wake structure of carangiform swimmers has been approximately reconstructed in the past from the results of PIV experiments using the velocity data in conjunction with the Helmhotz theorem for inviscid barotropic fluids. For carangiform swimmers with a single vortex row wake structure, a series of connected vortex rings has been suggested(Lighthill, 1969), while for anguilliform (eel-like) swimmers, which exhibit a double vortex row, two disconnected vortex rings have been hypothesized(Muller et al., 2001). As seen in Figs 10, 11, 12, both types of 3D vortical structures are observed in the simulations depending on the

*St*. Nevertheless, the rather simple wake structure that was hypothesized in previous experiments is observed in our simulations only for the lowest

*Re*case (Fig. 10). For this case, the double row wake consists of laterally dislocated vortex loops while the single row wake consists of a train of inverted hairpin-like vortices braided together such that the legs of each vortex are attached to the head of the preceding vortex. This wake structure resembles the 3D wake structure observed in laboratory experiments with vibrating spheres(Govardhan and Williamson,2005) and flapping foils(Buchholz and Smits, 2006). For

*Re*=4000 (see Fig. 11), the wake structure becomes significantly more complex. The single row wake consists of braided hairpins as for the

*Re*=300 case,but the hairpins in this case have longer and more stretched legs and more slender heads. In addition, smaller-scale vortical structures attaching to the hairpin vortices are beginning to emerge in the wake. The double row pattern is no longer made up of two rows of simple vortex loops but consists of very complex and highly 3D coherent structures connected together through complex columnar structures. To further appreciate the enormous complexity of the wake structure at this combination of

*Re*and

*St*, we show in Fig. 13 different views of the

*q*iso-surfaces superimposed to out-of-plane vorticity contours. Finally, for the inviscid case the single row wake consists of connected vortex loops, which are more flat in shape and stretched in the streamwise direction. The double row, on the other hand, exhibits smaller structures than the

*Re*=4000 case and the complexity of the wake has increased further.

In addition to the highly 3D structure of the wake vortices, Figs 10, 11, 12 further underscore the effect of *St* on the wake patterns. For *St* outside of the *Re*-dependent range within which the wake transition occurs, the *St* is clearly the dominant parameter that governs the 3D wake structure regardless of the *Re*.

While single row, thrust type wakes have been observed in all experiments involving carangiform swimmers, double row wake patterns have never been observed experimentally for this mode of swimming. Recall, however, that carangiform swimmer in nature undulate (for reasons we have already clarified above) their bodies at *St* in the range *St*=0.2–0.35,which is well within the range of *St* for which single row wakes are observed in our simulations. The double row wake structure has been observed in experiments with anguilliform swimmers(Muller et al., 2000; Tytell and Lauder, 2004). Muller et al. tested eels swimming at *St*=0.31, 0.33 and 0.56, and postulated that the double vortex row is produced by two consecutively shed ipsilateral body and tail vortices, which combine to form a vortex pair that moves away from the mean path of motion(Muller et al., 2000). However, Tytell and Lauder, who studied eel swimming at *St*=0.314 using high-resolution PIV data, suggested that the body anterior to the tail tip produces relatively low vorticity and the wake structure results from the instability of the shear layers separating the lateral jets, reflecting pulses of high vorticity shed at the tail tip(Tytell and Lauder, 2004). Moreover, Tytell and Lauder hypothesized that the difference in the wake structure of the carangiform and anguilliform swimmers comes from their difference in body shape and not in kinematics, i.e. if a mackerel was going to swim like an eel, the wake patterns would still be different(Tytell and Lauder, 2004). This hypothesis can be easily checked by numerical simulations using our present model, by imposing anguilliform kinematics on a mackerel body. Such computations are currently under way and will be reported in a future communication.

Nevertheless, our results do suggest that the wake structure (double row *vs* single row) is not that dependent on the body shape, as both wake structures have been observed from identical body shapes. Therefore, the wake structure is expected to depend on flow parameters such as *St* and *Re* and primarily on the *St*. This conclusion is also supported by a number of recent experimental and computational studies with flapping foils. For example, the flapping foil flow visualization experiments(Buchholz and Smits, 2005) for *Re*=640 also pointed to the *St* as the key parameter governing the wake structure. Buchholz and Smits observed that for *St* between 0.2 and 0.25 the wake consists of a single row of vortices, while for higher *St* the wake splits laterally, forming two separate trains of vortex structure, i.e. double vortex row (Buchholz and Smits, 2005). Experimental results(von Ellenrieder et al., 2003)and computations (Blondeaux et al.,2005) for a rectangular flapping foil at *Re*=164 and *St* in the range 0.175<*St*<0.4 show that as the *St* increases within this range, vortex-to-vortex interactions intensify and the wake becomes more complex but remains a single row structure. Numerical simulations (Dong et al., 2006) for a flapping ellipsoidal foil revealed the emergence of both single and double row wake structures. It is also important to point out that the single row structure has been produced by the vortex induced vibrations of a sphere (Govardhan and Williamson, 2005) and, as pointed out by Buchholz and Smits(Buchholz and Smits, 2006),this is an excellent example of the ubiquity of this type of wake structure,despite great differences in geometry.

### Concluding remarks

In the present study we constructed a virtual carangiform swimmer and employed it to elucidate the hydrodynamics of this type of locomotion and clarify and reconcile the results of laboratory experiments with live fish. The virtual tethered swimmer allowed us to perform controlled numerical experiments by systematically varying the Reynolds and Strouhal numbers while keeping the swimming kinematics fixed. As such, we were able to pose and answer questions that cannot be tackled experimentally due to the inherent difficulties in performing and analyzing the results of controlled experiments with live fish. The most important findings of our work are summarized as follows.

For a given *Re* there is a unique *St, St*^{*}, at which body undulations produce sufficient thrust to exactly cancel the hydrodynamic drag, making constant-speed self-propulsion possible. This is an important finding as it suggests that fish may not be selecting the *St* at which they undulate their bodies solely based on efficiency considerations but also because this is the only *St* at which they can swim steadily at a given speed.

*St*^{*} is a decreasing function of *Re* and approaches the range at which carangiform swimmers swim in nature(∼0.2–0.35) only in the limit *Re*=∞.

The Froude efficiency based on the thrust force at *St*^{*}increases with *Re*, suggesting that carangiform kinematics becomes a more efficient mode of aquatic locomotion only in the inertial regime. This finding is consistent with the fact that carangiform kinematics is the preferred mode of locomotion for fast (high *Re*) swimmers.

The power required for undulatory swimming, which for steady swimming is equal to the power of the lateral body undulations, was found to significantly exceed the power required to tow the rigid body at the same speed. Furthermore, the swimming power was found to be a decreasing function of *Re*. Therefore, the swimming power can be used instead of the Froude efficiency to explain why carangiform kinematics is preferred in nature by fast swimmers. Given the ambiguities involved in the definition and computation of the Froude efficiency, the swimming power provides an objective and unambiguous measure for quantifying the energetics of different modes of swimming and should be used for this purpose, as also suggested elsewhere(Schultz and Webb, 2002).

At a given *Re*, undulatory motion is shown to increase the friction drag above the rigid body level with *St* while only initially increasing the form drag. At *St*∼0.25, the form drag falls below the rigid body form drag level and monotonically decreases with *St*afterwards. The friction drag was found to be the dominant portion of the total drag for all *Re* in our simulations. At lower *Re*(∼300) the variation of the total drag with *St* mimics that of the friction drag i.e. increases with *St*. However, at sufficiently high *Re* (∼10^{3} and higher) the total drag mimics that of the form drag, i.e. initially increases then decreases with *St* because at the higher *Re* the increase of friction drag by the body undulations is very moderate, with the friction drag increasing mildly and eventually asymptoting toward a constant value above a threshold *St*. Consequently, at sufficiently high *Re* the total drag force mimics the form drag, which is effectively reduced below the rigid body level by the undulatory motion for *St*>0.2.

The initial increase of the form drag with *St* above the rigid body level occurs within the range of *St* for which the phase velocity of the undulatory body wave *V* is less than the flow velocity *U*(*St*<0.22). In this case the undulatory body motion acts to retard the near-wall flow relative to the free-stream and leads to the onset of separated flow in the posterior of the body, which accounts for the increased drag force. As *St* is increased further the body wave phase velocity ultimately exceeds the flow velocity (*St*>0.22). Under these conditions the fish body acts like a piston that acts to accelerate the flow backward relative to the free stream flow. Separation is eliminated and a pocket of positive (stagnation) pressure forms in the posterior of the body as a result of the transfer of energy from the beating tail to the flow, which explains the observed monotonic reduction of form drag for *St*>0.2.

The 3D structure of the wake is shown to depend primarily on the *St*. At all *Re* a wake with a single row of vortices resulted at low *St*. At higher *St* a more complex and laterally diverging wake structure with a double row of vortices was observed. The *St* range within which the transition from the single to the double row wake structure occurs was found to depend on the *Re*. The double row wake structure has not been observed before for carangiform swimmers because such fishes tend to swim at low *St* (∼0.2–0.35), for which the single vortex row wake structure dominates. However, our results are entirely consistent with numerical simulations and experiments with flapping foils.

Finally, it important to recognize and comment on the fact that, in addition to the morphological and kinematical parameters we considered herein,there are other parameters that could potentially affect the functional dependence of *St*^{*} on *Re*. First, the *St*^{*} should be expected to depend on the morphological characteristics of a fish and as such it should vary among different fishes. Second, the shape parameters (wavelength and the amplitude envelope), which were fixed in this study, can also affect the *St*^{*}. The wavelength λ/*L* of carangiform swimmers is in the range of 89%to 110% (Videler and Wardle,1991); λ/*L=95*% was used in this study. Higher wavelengths can increase the traveling wave speed, which in turn could slightly reduce *St*^{*}. Higher tail beat amplitudes(amplitude envelope) may also lead to lower *St*^{*}. The assumption that the tail beat amplitude stays constant and does not change with *Re* is valid for many carangiform swimmers, but certain fishes may increase the tail beat amplitude with speed (*Re*), e.g. the tail beat amplitude increases with speed in chub mackerel but not in kawakawa tuna(Donley and Dickson, 2000). The effect of dorso-ventral asymmetry (top–bottom) of the tail movement(Gibb et al., 1999) has not been considered in this study either. Another aspect of carangiform kinematics that we did not address in this work is the detailed characterization of the motion of the tail fin, which in certain fishes such as tuna could introduce additional kinematical parameters (Zhu et al., 2002). While the amplitude of the root point of the tail is the same as that of the body at the point of junction, the tail can form a different angle than the tangent to the body. Also, this angular motion may have a non-zero phase angle with respect to the body motion. These parameters have been considered in the inviscid numerical simulations of tuna swimming(Zhu et al., 2002) and were shown to play an important role in the dynamic interactions of vortices shed by the body and the tail fin. Additional studies will be needed to investigate and quantify the effect of the aforementioned parameters on the hydrodynamics of carangiform swimmers. These studies are beyond the scope of this paper and will be pursued as part of our future work.

## APPENDIX

### Additional materials and methods

In this Appendix we report results aimed at validating the ability of our numerical method to accurately predict the forces acting on a moving body. We also report the results of a grid refinement study for the fish simulations to quantify the effect of numerical resolution on the accuracy of the results we reported in this paper.

### Validation study: forced inline oscillations of a cylinder in a fluid initially at rest

The ability of our numerical method to predict forces acting on a moving body has already been demonstrated(Gilmanov and Sotiropoulos,2005). To further validate the ability of our method to predict the hydrodynamic force given by Eqn 15 and its pressure and viscous contributions, we consider herein the case of a circular cylinder starting to oscillate in the horizontal direction in a fluid initially at rest. Benchmark experimental and computational results for this case have been reported(Dutsch et al., 1998).

*x*is the location of the center of the cylinder,

_{c}*f*is the oscillation frequency, and

*A*

_{m}is the oscillation amplitude. The flow induced by such oscillations is governed by two non-dimensional parameters: (1) the Reynolds number

*Re*=

*U*

_{m}

*D*/ν, based on the maximum oscillation velocity

*U*

_{m}, cylinder diameter

*D*, and the fluid kinematics viscosity ν; and (2) the Kuelegan–Carpenter number KC=

*U*

_{m}/

*fD*. According to Eqn A1, the Kuelegan–Carpenter number is equal to KC=2π

*A*

_{m}/

*D.*The computations are performed at

*Re*=100 and KC=5, for which both experimental and numerical results have been reported (Dutsch et al.,1998). The domain is discretized with a mesh consisting of 721×481 nodes in the inline (oscillatory) and transverse directions,respectively. 300×100 nodes are distributed uniformly in a 3

*D*×

*D*box, which contains the cylinder during the oscillations. The domain outer boundaries are placed 50

*D*from the initial position of the cylinder, and Neumann boundary condition(∂

*u*/∂

_{i}*n*=0), where

_{j}*n*is the normal to the outer boundary surface) has been used.

_{j}Fig. A1 compares the calculated (in red) inline hydrodynamic force (solid lines) and its pressure(dotted lines) and viscous (broken lines) components with published computations (Dutsch et al.,1998) (in black). It is clear that the calculated forces are in excellent agreement with Dutsch et al.'s results(Dutsch et al., 1998).

Fig. A2 compares the inline velocity profiles at *x*_{1}=–0.6*D* at three different phase angles (ϕ=2π*ft*) calculated by our method and measured by Dutsch et al. (Dutsch et al.,1998). As for the force comparisons, the velocity profiles are in excellent agreement with the measurements.

Finally, Fig. A3 shows the calculated instantaneous vorticity field at four different phase angles, which is dominated by two counter-rotating vortices. The computed results are identical to Dutsch et al.'s computational results(Dutsch et al., 1998), which are not shown herein but are reported in the same format in their paper.

### Numerical sensitivity studies for the fish swimming simulations

To investigate the sensitivity of the fish swimming simulations on the grid spacing, we performed simulations on four, successively finer meshes. The uniform spacing in the cuboid containing the fish for the four different meshes is *h*=0.016, 0.012, 0.008 and 0.0004, corresponding to 0.7,2.0, 5.5 and 14.7 million total grid nodes, respectively; we shall refer to these grids as A, B, C and D, respectively. The grid sensitivity study is carried out for *Re*=4000 since this is the most challenging case relative to other two *Re* we simulate in this work due to the need to accurately resolve the boundary layer along the fish body. For each grid, the domain size, boundary conditions and the time step are the same as those used for grid C, which is the grid used for all simulations reported so far in this paper (see Materials and methods). Fig. A4 shows the effect of grid size on the time averaged force coefficient. The figure includes: (1) the results obtained on grid C for all simulated *St* reported in Fig. 4; (2) the results obtained on grids A, B and D for *St*=0.5; and (3) the results obtained on the finest grid D for *St*=0.6. It is evident from this figure that as the grid is refined for the *St*=0.5 case the computed results converge monotonically toward the grid independent solution. Grids A and B are too coarse to obtain accurate results but the results obtained on the two finest meshes (C and D)are very close to each other. Furthermore, the results for *St*=0.6 obtained on the two finest meshes show that the values of *St*^{*} determined on these two meshes are very close to each other. Based on these results we conclude that the results obtained on grid C,even though they are not strictly grid independent, are insensitive to further grid refinement and little is to be gained in terms of overall accuracy by adopting mesh D.

To demonstrate the rate at which the numerical method converges toward the grid independent solution, we plot in Fig. A5 the error in the total force coefficient as a function of the grid spacing in log–log scale. The error for each grid is calculated using the solution on grid D as the `exact' solution. It is evident that the rate of convergence of our method is slightly better than second order. This finding confirms again that the HCIB method is 2nd order accurate, as previously shown (Gilmanov and Sotiropoulos, 2005).

To demonstrate the effect of grid refinement on the instantaneous force, we plot in Fig. A6 the total force time history on grids B, C and D for *St*=0.5. As already concluded above, grid B is too coarse to obtain accurate results. The computed force time series on this grid exhibits spurious high frequency oscillations,which are due to interpolation errors induced by the very coarse mesh in the vicinity of the immersed boundary. As expected, the amplitude of these high frequency oscillations is seen to diminish with grid refinement. The oscillations are drastically reduced on grid C and practically eliminated on grid D.

## Acknowledgements

This work was supported by NSF Grants 0625976 and EAR-0120914 (as part of the National Center for Earth-surface Dynamics) and the Minnesota Supercomputing Institute. The authors greatly appreciate the anonymous referees for their insightful comments, which helped improve this paper.

## References

**Anderson, E., McGillis, W. and Grosenbaugh, M.**(

**Barrett, D. S., Triantafyllou, M. S. and Yue, D. K. P.**(

**Blondeaux, P., Fornarelli, F., Guglielmini, L., Triantafyllou,M. S. and Verzicco, R.**(

**Borazjani, I., Ge, L. and Sotiropoulos, F.**(

**Breder, C. M.**(

**Buchholz, J. H. and Smits, A. J.**(

**Buchholz, J. H. J. and Smits, A. J.**(

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

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

**Dabiri, J. O.**(

**Dong, H., Mittal, R. and Najjar, F. M.**(

**Donley, J. M. and Dickson, K. A.**(

*Euthynnus affinis*) and chub mackerel (

*Scomber japonicus*).

**Drucker, E. G. and Lauder, G. V.**(

**Dutsch, H., Durst, F., Becker, S. and Lienhart, H.**(

**Fish, F. E.**(

*Tursiops truncatus*).

**Fish, F. E. and Lauder, G. V.**(

**Fish, F. E., Innes, S. and Ronald, K.**(

**Ge, L. and Sotiropoulos, F.**(

**Gibb, A. C., Dickson, K. A. and Lauder, G. V.**(

*Scomber japonicus*: testing the homocercal tail model of fish propulsion.

**Gilmanov, A. and Sotiropoulos, F.**(

**Govardhan, R. N. and Williamson, C. H. K.**(

**Gray, J.**(

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

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

**Hunt, J. C. R., Wray, A. A. and Moin, P.**(

**Jordan, D. S., Evermann, B. W. and Clark, H. W.**(

*US Commission for Fish and Fisheries*.

**Kern, S. and Koumoutsakos, P.**(

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

**Lighthill, M. J.**(

**Lighthill, M. J.**(

**Lighthill, M. J.**(

**Lighthill, M. J.**(

**Lindsey, C. C.**(

**Liu, H.**(

**Liu, H., Wassersug, R. J. and Kawachi, K.**(

**Liu, H. and Wassersug, R.**(

**Liu, H. and Kawachi, K.**(

**Muller, U. K., Van Den Heuvel, B. L. E., Stamhuis, E. J. and Videler, J. J.**(

*Chelon labrosus risso*).

**Muller, U. K., Stamhuis, E. J. and Videler, J. J.**(

**Muller, U. K., Smit, J., Stamhuis, E. J. and Videler, J. J.**(

*Anguilla anguilla*).

**Nauen, J. C. and Lauder, G. V.**(

*Scomber japonicus*.

**Nauen, J. C. and Lauder, G. V.**(

*Scomber japonicus*(Scombridae).

**Nelson, J. S.**(

**Panton, R. L.**(

**Ramamurti, R., Sandberg, W. C., Lohner, R., Walker, J. A. and Westneat, M. W.**(

**Rosen, M. W.**(

**Saad, Y.**(

**Satish Balay, Kris Buschelman William, D., Gropp Dinesh Kaushik Matthew, G., Knepley Lois, Curfman McInnes Barry, F., Smith. and Zhang, H.**(

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

**Sfakiotakis, M., Lane, D. M. and Davies, J. B. C.**(

**Shen, L., Zhang, X., Yue, D. K. P. and Triantafyllou, M. S.**(

**Taneda, S. and Tomonari, Y.**(

**Triantafyllou, M. S. and Triantafyllou, G. S.**(

**Triantafyllou, M. S., Triantafyllou, G. S. and Gopalkrishnan,R.**(

**Triantafyllou, M. S., Triantafyllou, G. S. and Yue, D. K. P.**(

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

**Videler, J. J.**(

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

*Pollachius virens*) and mackerel (

*Scomber scombrus*): A kinematic analysis.

**Videler, J. J. and Wardle, C. S.**(

**von Ellenrieder, K. D., Parker, K. and Soria, J.**(

**Weihs, D.**(

**Weihs, D.**(

**Wolfgang, M. J., Anderson, J. M., Grosenbaugh, M. A., Yue, D. K. and Triantafyllou, M. S.**(

**Wu, T. Y.-T.**(

**Wu, T. Y.-T.**(

**Zhu, Q., Wolfgang, M. J. and Yue, D. K. P.**(