Numerical investigation of the hydrodynamics of carangiform swimming in the transitional and inertial flow regimes

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

As carangiform swimmers typically achieve high swimming speeds, their motion is characterized by very high Reynolds numbers, Re>10⁴(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⁵ 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 Stin a computational grid with about 1.5×10⁵ 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⁵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⁶–10⁷ 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 Uconstant, 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 Reand 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 Rethe 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.

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⁶ 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⁷ 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).

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

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=2h_max, and the tail beat frequency f: St=2fh_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.

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/2a_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).

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 2L×L×7Lcuboid 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 2L and height L are ten times the fish width 0.2L and height 0.1L, respectively. A uniform mesh with constant spacing h=0.008L 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.5L 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 2L×L×11L. Similar to the shorter (7L long) domain, a uniform mesh with constant spacing h=0.008L 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 7L 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.

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.

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.

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.

1. 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).

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

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

4. 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 Stincreases 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.

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.

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.

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³L² and the values are reported in Table 2.

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.

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>Useparation 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 Stsufficiently 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.

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=½(∥Ω∥²–∥S∥²),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.

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 Stafterwards. 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³ 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.

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.

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

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₁=–0.6D 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.

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.

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.