The hydrodynamics of ribbon-fin propulsion during impulsive motion

Weakly electric knifefish have been studied for several decades to gain insights into how vertebrates process sensory information (for reviews, see Bullock, 1986; Turner et al., 1999). Knifefish continuously emit a weak electric field, which is perturbed by objects that enter the field and distort the field because their electrical properties differ from the surrounding fluid. These perturbations are detected by thousands of electroreceptors on the surface of the body. MacIver and coworkers (Snyder et al.,2007) have shown that the knifefish is able to both sense and move omnidirectionally. The primary thruster that weakly electric fish use in achieving their remarkable maneuverability is an elongated anal fin(Fig. 1), which we generically refer to as a ribbon fin.

Approximately 150 species of South American electric fish (family Gymnotidae) swim by using a ribbon fin positioned along the ventral midline. This fin is also used for swimming in one weakly electric species in Africa, Gymnarchus niloticus, where the fin is positioned on the dorsal midline, and in species of the non-weakly electric family Notopteridae, where the fin is positioned along the ventral midline. In the present study, we address the mechanical principles of force generation by the ribbon fin in the context of the South American weakly electric black ghost knifefish, Apteronotus albifrons (Linnaeus 1766)(Fig. 1A).

Knifefish swim by passing traveling waves along the ribbon fin. The waveform is often similar in overall shape to a sinusoid(Fig. 1B). The body is typically held straight and semi-rigid while swimming(Fig. 1B), i.e. body deformations are small compared with fin deformations. This may facilitate sensory performance as the body houses the electric field generator, and movement of the tail causes modulations of the field that are more than a factor of 10 larger than prey-related modulations(Chen et al., 2005; Nelson and MacIver, 1999). Knifefish frequently reverse the direction of movement without turning by changing the direction of the traveling wave on the fin and are as agile swimming backward as they are swimming forward(Blake, 1983; Lannoo and Lannoo, 1993; MacIver et al., 2001; Nanjappa et al., 2000).

The ability to switch movement direction rapidly (in ≈100 ms)(MacIver et al., 2001) is integral to several behaviors. Previous work by MacIver et al.(MacIver et al., 2001) has shown that prey are usually detected while swimming forward, after the prey has passed the head region, and the fish then rapidly reverses the body movement to bring the mouth to the prey during the prey strike. During inspection of novel objects, the fish are observed to engage in forward–backward scanning motions(Assad et al., 1999), which may be important for increasing spatial acuity(Babineau et al., 2007).

Ribbon-fin-based swimming is commonly referred to as the gymnotiform mode by Breder (Breder, 1926). In addition to being agile, prior research on ribbon-finned swimmers has also suggested that they are highly efficient for movement at low velocities(Blake, 1983; Lighthill and Blake, 1990). This claim is supported by the discovery that these fish use half the amount of oxygen per unit time and mass as non-gymnotid teleosts(Julian et al., 2003).

Two goals motivate the current study. First, in order to advance from the mature understanding we have of sensory signal processing in weakly electric knifefishes to an understanding of how these signals are processed to control movement, we need to characterize the hydrodynamics of ribbon-fin propulsion. Second, artificial ribbon fins may provide a superior actuator for use in highly maneuverable underwater vehicles for applications such as environmental monitoring (Epstein et al.,2006; MacIver et al.,2004).

We use computational fluid dynamics to examine the flow structures and forces arising from a sinusoidally actuated ribbon fin. We compare the computed flow structures with those measured from a robotic ribbon fin using digital particle image velocimetry (DPIV) and compare the computed surge force with the drag force measured from towing a cast of the fish. Whereas tow drag can provide a useful estimate of the thrust needed during steady swimming, for the impulsive motions modeled in this study, the thrust needed to undergo typical accelerations is more directly relevant. Thus, we also compare computed forces with the thrust that we estimate is needed for two different types of swimming direction reversals: reversals that occur during prey capture strikes from kinematic data collected in a previous study(MacIver et al., 2001) and reversals that occur during refuge tracking behavior, where fish placed in a sinusoidally oscillating refuge will move to maintain constant position with respect to the refuge.

For the present study, we idealize the fin kinematics as a traveling sinusoid on an otherwise stationary (i.e. non-translating, non-rotating)membrane (Fig. 1C). As a consequence, the top edge of the fin remains fixed at all times, and all points on the fin below this edge move in a sinusoidal manner. The fin deformation is specified in Eqn 1below. As indicated in Fig. 1C,positive surge is defined as the force on the fin from the fluid in the direction from the tail to the head. If the traveling wave passes from the tail to the head, then the force on the fin from the fluid would be from the head to tail, corresponding to negative surge. Positive heave is vertically upward.

We chose to characterize the hydrodynamics of a fixed fin in a stationary flow because this is most relevant for understanding forces arising from maneuvering movements that occur when the body is at near-zero velocity with respect to the fluid far away from the body. In future work, using this approach will also allow us to compare our simulated force estimates with those obtained empirically from a robotic ribbon fin placed on a linear track pushing against a force sensor (Epstein et al., 2006). In subsequent studies, we will be examining the hydrodynamics of a stationary fin under imposed flow conditions and when the fin and an attached body are allowed to self-propel through the fluid.

Flow visualizations from computational simulations and DPIV indicate that the mechanism of thrust generation is a streamwise central jet and associated attached vortex rings. We show that, despite the lack of cylindrical symmetry in the morphology of the fin, its peculiar deformation pattern – the traveling wave – produces a jet flow often found in highly symmetric animal forms, such as jellyfish and squid. Whereas previous research focused exclusively on the surge force (Blake,1983; Lighthill and Blake,1990), we find that the ribbon fin is also able to generate a heave force, which pushes the body up. This arises from the generation of counter-rotating axial vortex pairs that are shed downward and laterally from the bottom edge of the fin. We hypothesize that the slanted angle of the fin base with respect to the spine observed in many gymnotids (e.g. Fig. 1A) leverages this heave force for forward translation. We also find that, in certain cases, as the number of waves on the fin decreases to below approximately two-thirds, the heave force surpasses the surge force. This switchover from an undulatory parallel thrust mode to an oscillatory normal thrust mode may provide insight into how the position and orientation of median fins varies with the length of the fin and the number of wavelengths that can be placed on it.

We show how the surge force from the fin scales as a function of a few key parameters. We found that for a stationary fin without imposed flow: (1) the surge force is proportional to (frequency)²× (angular amplitude)^3.5×(fin height/fin length)^3.9×Φ(wavelength/fin length), where Φ is a function that approximates the variation of surge force as a function of wavelength normalized with fin length; (2) for angular deflections aboveθ=10 deg., where θ is defined in Fig. 1C, previous analytical work (Lighthill and Blake,1990) underestimated the magnitude of surge force; (3) surge force shows a peak when the wavelength is approximately half the fin length, similar to what is observed biologically (Blake,1983) and contrary to the monotonic increase in force with wavelength predicted by the analytical results of Lighthill and Blake(Lighthill and Blake, 1990);and (4) the computed surge force compares well with empirical estimates based on fin kinematics and accelerations during swimming direction reversals, as well as with body drag measurements.

This completes the two-step algorithm at a given time step. The same procedure is repeated at each time step for the desired total simulated time,typically the time course of several traveling waves passing down the length of the fin.

Force on the fin is computed as the rate of change of momentum of the fluid in the correction step (Eqn 5). The first step of obtaining û (Eqns 3 and 4), as well as Eqns 7 and 8, conserve momentum. Hence, all of the change in the total fluid momentum comes from the immersed boundary condition of the correction step in Eqn 5. This change is computed as the difference between the total momentum before and after the correction step.

In the immersed boundary approach, a single regular grid can be used to solve the fluid and pressure equations at all times. This grid need not be body conforming. An alternative approach is to use a body-conforming method in which the grid conforms to the solid body surface. The immersed boundary approach has certain advantages over the body-conforming grids approach. In the latter case, as the fin shape is changing with time, remeshing is often necessary at every time step. This adds to the computational cost. Also, the use of fast elliptic solvers (e.g. fast Fourier transform-based) for pressure is not straightforward for complex meshes. Lastly, in the immersed boundary approach, the velocity correction step does not require significant computer time.

The correction step (Eqn 5) is implemented by using Lagrangian marker particles representing the fin. The prescribed fin deformation velocity is known at these marker particle locations. The velocity correction Eqn 5 would be exact if the grid points on which the fluid equations are solved coincided with the fin marker particle locations. In general, this is not true because the fin geometry is usually complex(Fig. 2). Therefore, when implementing Eqn 5, we interpolate the fin velocity from the marker particle location to its neighboring fluid grid points. The interpolation is carried out using a top-hat interpolation function as follows. For each Lagrangian marker particle, we determine which eight Eulerian grid nodes surround it. The velocity u_fin at the marker particle is then assigned to those Eulerian nodes. For nodes that have contributions from multiple marker particles, the arithmetic mean of the values from all contributing marker particles is assigned to that node. For Eulerian nodes in the fluid domain, no particle contribution results, hence the values at these nodes are left unchanged during the correction step.

The non-conformity of the fin geometry with the fluid grid and the resulting need for interpolation lead to some smearing of the velocity field near the fin. Refined grids would likely be required if the flow is turbulent near the fin. This increases the computational cost. In this study, we will consider Reynolds numbers (Re) corresponding to typical adult black ghost knifefish, which are approximately 5600 based on the wavelength and the wave speed. This Re is moderate and the flow is anticipated to be laminar, although in some cases the flow may allow transition to turbulence.

The above numerical scheme gives the velocity and pressure fields in the fluid domain at each time step. This scheme is used to set up the ribbon-fin problem where the fin is fixed along the top edge. We consider a fin size that matches that of an adult black ghost knifefish(Fig. 1A).

The physical parameters involved are the geometric and kinematic parameters of the fin, and the fluid properties. They are: L_fin, h_fin, f, θ_max, λ, ρ,μ, which are, respectively, the length of the fin, the height of the fin(see Fig. 1), the frequency of the traveling wave, the maximum angular deflection of the fin from the midsagittal plane, the wavelength of the traveling wave, the fluid density and the dynamic viscosity. Typical values of the geometric parameters for an adult(15 cm in length) black ghost knifefish (Apteronotus albifrons) are L_fin≈10 cm and h_fin≈1 cm. These values are taken from photographs of live specimens in our laboratory. These geometric values were the ones used for all simulations performed in this study, except for one series of numerical experiments where we varied h_fin. Typical kinematic parameters are f≈3 Hz,θ _max≈30 deg. and λ≈5 cm(Blake, 1983). The density and viscosity of water are taken as ρ=1000 kg m^–3 andμ=8.9×10^–4 Pas.

Dynamic similarity tells us that the force on the fin from the fluid does not independently depend on these parameters but is a function of the dimensionless groups of these parameters. Using dimensional analysis, these dimensionless groups can be written as:(ρfλ²)/(2πμ), θ_max, h_fin/L_fin,λ/L_fin.

In the following sections, we study the effect of these nondimensional parameters by changing the physical parameters as follows. We vary the Re = (ρfλ²)/(2πμ), by changing f and keeping all the other physical variables fixed. Similarly, we vary h_fin/L_fin by varying h_fin, and vary the specific wavelengthλ/L_fin by choosing different values of λ. The range of kinematic and geometric parameters numerically investigated is shown in Table 1. The Re for the baseline case is 894. Note that in all cases, the fin is simulated without a fish body (see Materials and methods) and is approximated as an infinitesimally thin membrane, within the smearing over 2–3 grid cells caused by interpolation on the computational grid (0.7–1.1 mm).

All of the simulations were performed using the San Diego Supercomputer Center's IA-64 Linux Cluster, which has 262 compute nodes each consisting of two 1.5 GHz Intel Itanium 2 processors running SuSE Linux. For interprocessor communication, the cluster uses the Myrinet 2000 gigabit ethernet interconnect network. The computational fluid dynamics code was written in Fortran 90 and C and uses the FFTW Library(www.fftw.org)for fast Fourier transforms.

We used DPIV to visualize the flow field in the coronal (horizontal) plane of a robotic ribbon fin (Fig. 3). The ribbon fin was 23.5 cm long and 7.0 cm high. The experimental details are described in Epstein et al.(Epstein et al., 2006).

The water in the tank was seeded with ≈3 mg l^–1silver-coated glass spheres with a mean diameter of 16 μm and a mean density of 1.4 g cm^–3. In the water tank, the working area was 120×35×34 cm. A 25 mJ Nd:YAG laser (New Wave Research,Fremont, CA, USA) was used to illuminate the particles. The laser beam was synchronized with a 1 megapixel charge couple device (CCD) camera (TSI,Shoreview, MN, USA) to record the sequence of images. For a better flow field resolution, the field of view for the images covered only the trailing 50% of the fin.

The time interval between image pairs was 750 μs, with image pairs obtained at a rate of 15 Hz. The laser sheet was approximately 5 mm below the bottom edge of the ribbon fin. The velocity field was obtained by analyzing each image pair with commercial particle velocimetry software (Insight, TSI). These velocity fields were then post-processed and analyzed using MATLAB (The Mathworks, Natick, MA, USA).

For the flow visualization, the robotic ribbon fin was operated at f=1.5 Hz, θ_max=25 deg.,λ/L_fin=1.4 and h_fin/L_fin=0.3.

An accurate urethane cast of a 185-mm long Apteronotus albifronswas used from a previous study (MacIver and Nelson, 2000). It was bolted to a rigid rod suspended from a custom force balance that used three miniature beam load cells (MB-5-89,Interface, Scottsdale AZ, USA). For force balance and calibration details, see Ringuette et al. (Ringuette et al.,2007). The fish cast was then towed through a tow tank that was 450×96×78 cm in length, width and depth (GALCIT towtank, Caltech,Pasadena, CA, USA) using a gantry system driven by a speed-controlled DC servomotor above the tank. Details of the towtank and gantry system can be found in Ringuette et al. (Ringuette et al., 2007). Trials were conducted at three speeds: 10, 12 and 15 cm s^–1. Drag was measured with the cast at two different orientations: (1) head-first, with the long axis of the body (spine) parallel to the oncoming flow (pitch=0 deg.), and (2) tail-first, with the spine parallel to the oncoming flow (pitch=0 deg.). Only the data collected after the startup force transient had settled were analyzed, until just before the end of the towing distance (300 cm). The data were filtered with a digital Butterworth low-pass filter (cutoff at 5 Hz) to remove transducer transients prior to further statistical analysis.

We estimated the thrust needed for two types of swimming direction reversals: (1) rapid reversals made during prey strike(MacIver et al., 2001) and (2)slower reversals that occur during refuge tracking behavior(Cowan and Fortune, 2007). These thrust estimates were based on motion capture data from a prior study(MacIver et al., 2001),measurements of refuge tracking movements and fin kinematics, and the effective mass of the fish (mass and added mass).

During prey-capture reversals, the body is initially moving forward at approximately 10 cm s^–1 while searching for prey. Following prey detection, the traveling wave on the fin reverses and the fish rapidly decelerates to reverse the direction of body movement and bring the mouth to the prey. We analyzed reversal accelerations across 116 prey-capture trials using methods described previously(MacIver et al., 2001). Trials within one-half of one standard deviation (number of trials N=45)were selected for quantification of the ribbon-fin traveling wave frequency immediately after the moment of body direction reversal.

In refuge tracking behavior, a fish within a refuge (such as a transparent plastic tube) will attempt to maintain a fixed relationship with respect to the refuge. Thus, if the tube is moved sinusoidally, the fish will too. The kinematics of this behavior have been analyzed previously(Cowan and Fortune, 2007) but without measurement of traveling wave frequency and for a different species of knifefish. Thus, to collect preliminary data for comparison with our computed forces, we placed an adult Apteronotus albifrons underwater in a tube suspended from a custom XY robot described elsewhere(Solberg et al., 2008). Digital video recordings were made of the refuge and fish while the refuge was moved sinusoidally with an amplitude of 5 cm at frequencies of 0.25–0.4 Hz. The video was manually inspected to estimate the fin frequency immediately following body reversals.

The geometry of the fin simulated in the present study approximates that of the fish examined in the prey capture and refuge tracking behaviors. We estimate that the ribbon fin was moved with amplitudes similar to 30 deg.,with a λ/L_fin of ∼0.5 and a L_fin of ∼10 cm. To approximate the thrust from the fin, we use Eqn 9 (to be introduced later) using these parameters and the traveling wave frequency measured from the video.

To estimate the thrust needed to accomplish the reversal, we multiply the effective mass previously estimated for the knifefish (10.4 g)(MacIver et al., 2004) by the accelerations extracted from the motion capture data of MacIver et al.(MacIver et al., 2001).

For our mixed Eulerian–Lagrangian immersed boundary approach, an array of validation tests was performed to ensure that both the Eulerian (for the fluid) and Lagrangian (for the fin) grid resolutions are sufficient to obtain accurate results. The flow solver was validated using test cases that allow comparison with an analytical solution or experimental data [for details, see Shirgaonkar and Lele(Shirgaonkar and Lele, 2007)]. Sensitivity tests with respect to the number of Lagrangian particles per grid cell (N_p) representing the fin showed that although the velocity field away from the body is relatively insensitive to N_p, N_p=8 is needed to obtain accurate fluid-solid surface forces. This is consistent with the previous observation that two particles are needed in each direction per grid cell to prevent`leakage' of fluid through the membrane(Peskin, 2002).

To validate the immersed boundary implementation, we simulated an impulsively started flat plate normal to the flow. For such a plate, flow separation occurs past the plate at both of its edges. This leads to the formation of twin vortices behind the flat plate. The development of such vortices was studied experimentally by Taneda and Honji(Taneda and Honji, 1971). For validation, we choose their case with Re=126 to compare with their experimental results. In this case, the Re is defined as Ud/v, where d is the length of the plate, U is the constant translational velocity of the plate and v is the kinematic viscosity. Fig. 4 depicts how the vortex length normalized by the plate length, s/d, grows as a function of time and compares the simulation results with those of Taneda and Honji (Taneda and Honji,1971). The vortex length is defined as the distance of the rear stagnation point from the plate (see Fig. 4). Although this validation test was performed at Re=126, as a part of verification of the code, we performed a grid convergence study at the actual fish Re of 5600 based on the wavelength and the wave speed. This ensures that the effect of higher Re, namely the thinning of boundary layers, is captured numerically.

For all simulations with the ribbon fin, we used a domain size of L_x, L_y, L_z=15 cm, 4 cm, 4 cm where the x, y and z axes are shown in Fig. 1C. A uniform Cartesian grid was used for simulations. For the grid sensitivity study, three different grid sizes, N_x, N_y, N_z=300, 80, 80(coarse), N_x, N_y, N_z=400, 100, 100(nominal) and N_x, N_y, N_z=440, 120,120 (fine), were used and the comparison of forces is shown in Fig. 5. It is seen that the nominal grid gives reasonably converged values of forces on the fin and,hence, it is the grid size chosen for all simulations reported later. It is noted that the forces can have some high-frequency noise, which is inherent to the immersed boundary method (Uhlmann,2003). To remove the noise and obtain physical values of the forces up to the grid-timescale, we filter the modes in the forces, which vary faster than the grid-timescale. Here, the grid-timescale is defined as the time required by a representative point on the fin to travel one grid cell. As the spatial discretization can only capture length scales up to the grid cell size, the above smoothing retains the forces down to the temporal scale consistent with the smallest resolved spatial scale.

The domain size used can potentially affect the force values due to the periodic boundary conditions used. To examine whether the above domain is sufficiently large, two different simulations with domain sizes L_x, L_y, L_z=15 cm, 4 cm, 4 cm and 17 cm, 6 cm, 6 cm were carried out. The forces on the fin with the two domains differed by less than 8%. Hence, we chose the shorter of the two domains above to minimize cost while maintaining adequate accuracy.

The vortex structure around the ribbon fin is shown in Fig. 6A and B where the isosurface of the total vorticity magnitude 30 s^–1 is shown and it is colored by the x-velocity. Movies 1 and 2 in supplementary material show the time evolution of the central jet underlying these vortex ring structures. Fig. 6Cdepicts the same vorticity isosurface together with the x-velocity isosurface u=2 cm s^–1. This shows the correlation of the high x-velocity regions and the vortex tubes. Fig. 7 shows x-velocity contours on a horizontal slice slightly below the lower edge of the fin and also the velocity vectors in the plane of the slice. The central axial jet along the surge direction is evident from the velocity vectors. The horizontal rolls of surge vorticity are seen in Fig. 8, which in later sections will be shown to be related to the successive upward and downward y-velocity near the fin surface(Fig. 9).

Fig. 10 shows the DPIV flow field of the robotic ribbon fin. Fig. 10 B1, B2 and B3 shows the flow field at t=2.67, 2.80 and 3.00 s, respectively. The color map represents the x-velocity normalized by the wave speed. The x-velocity is higher (red in the color map) in the concave part of the ribbon fin and close to the inflection points. These high x-velocity regions are advected to the posterior part of the ribbon fin at approximately the wave speed(V=fλ). These results confirm the presence of the streamwise jet found in the simulation results(Fig. 7).

Fig. 11 shows the temporal behavior of the surge, heave and sway forces (see Fig. 1 for definitions) from the simulation of the baseline case where f=2 Hz,θ _max=30 deg., λ/L_fin=0.5 and h_fin/L_fin=0.1. Surge and heave forces oscillate at twice the fin frequency, and sway force varies with the fin frequency. This is because at any given axial cross-section of the fin, the fluid is being pushed by the fin in the downward and streamwise directions at both extremities of the fin oscillations. The two peaks of surge forces in one fin period are unequal. This is a result of the asymmetric wake structure observed in Fig. 6C. Asymmetric wakes have also been shown to exist in an even simpler system consisting of an oscillating cylinder in cross-flow(Williamson and Roshko, 1988). Williamson and Roshko show that the wake structure (and hence forces) can significantly change with changes in system parameters (frequency and amplitude) or the flow history (Williamson and Roshko, 1988). Hence, to examine the sensitivity of our results, we varied the initial phase of the fin displacement from 0 deg. to 180 deg. in steps of 90 deg. It was found that the time-averaged mean forces are robust against the changes in the initial phase. In addition, the forces we report later do not show any sudden change in behavior with system parameters, leading us to believe that the nature of the vortex shedding in the wake is not sensitive to the system parameters for the range of parameters studied here.

For the characterization and scaling of the ribbon fin forces, we calculated the mean values of the fin forces over at least one cycle after the initial transient has passed, i.e. after the forces show approximately quasi-steady oscillatory behavior. It was found that at later times in the simulations, the finiteness of the domain causes a small drift in the mean forces due to the boundary effects. The averaging duration for mean force computation was selected to lie before the beginning of this late regime. The mean forces from all the ribbon fin simulations are shown in Fig. 12A–D, which show the results for the simulation Set 1, Set 2, Set 3 and Set 4, respectively(from Table 1). The mean sway forces are zero and are not shown in the plots. Surge force ranges from 0 to 1.85 mN over the parameter space examined.

Fig. 13 shows the drag measurements for the black ghost knifefish cast as a function of time for both backward and forward motion through the fluid. The drag force varied between 1 and 2 mN over the velocity range examined, 10–15 cm s^–1. This range is behaviorally relevant to black ghost knifefish, e.g. the knifefish have been observed to exhibit a mean swimming velocity of 10 cm s^–1 prior to prey detection(MacIver et al., 2001).

For the prey-capture study, the reversal acceleration was found to be 144±63 cm s^–2 (mean ± s.d.). Trials within one-half of one standard deviation (N=45) were selected for quantification of the ribbon-fin traveling wave frequency through frame-by-frame inspection of the video. Of these, 13 could not be quantified due to the shortness of post-reversal movement and video quality. For the remaining trials (N=32), the traveling wave frequency was 7±1 Hz.

The calculation of available thrust for this case, using the scaling laws described below, gave a value of 5.5±0.1 mN. The calculation of needed thrust using effective mass and measured acceleration gave 12.6±4 mN. We address the gap between these values in the Discussion.

For the preliminary refuge tracking data, we found reversal accelerations of 106±72 cm s^–2 (N=7). Across these accelerations, the fin traveling wave frequency was 4±0.4 Hz. The calculation of available thrust gave 2±0.4 mN. The calculation of needed thrust using effective mass and measured acceleration gave 1±0.7 mN.

First, we study the baseline case to examine the fluid dynamics near the fin. The flow features observed in the vicinity of the ribbon fin indicate certain mechanisms for surge and heave force generation.

Fig. 6A shows an instantaneous isosurface of vorticity magnitude colored by the axial velocity. Two distinct flow features are manifested here. (1) The fin generates a series of vortex rings on both of its surfaces. These rings are attached to the fin on both sides, and their vortex axes are located near the fin's midsagittal plane, a small distance below the fin, as depicted in Fig. 6B. (2) Associated with these crab-shaped vortex rings, there is a central jet along the fin axis,which represents the momentum imparted to the fluid by the fin along the direction of wave motion. The central jet is more prominently seen in the slice of x-velocity (Fig. 7). The central jet is also observed in the flow visualization from the DPIV results for the robotic ribbon fin (see Fig. 10).

It is generally known that vortex shedding behind the posterior of a swimming animal is the core signature of thrust generation, as found in previous studies of anguilliform swimming(Hultmark et al., 2007; Kern and Koumoutsakos, 2006; Tytell and Lauder, 2004). These shed vortices can be of different types, e.g. the `2-S' vortex tubes[i.e. two single, counter-rotating vortices for each periodic cycle (see Koochesfahani, 1989; Williamson and Roshko, 1988)]shed at the trailing edge of anguilliform swimmers such as eel(Tytell and Lauder, 2004) or vortex rings shed by jellyfish (Dabiri et al., 2005a; Dabiri et al.,2005b). Traditionally, the generation of vortex rings as a means of self-propulsion has been considered prominent mainly in animals that swim using a backward jet directly emanated from their interior surfaces (e.g. jellyfish and squid). Despite lacking the high degree of cylindrical symmetry that is found in the jet-generating surfaces of such organisms, the ribbon fin is able to produce a significantly strong central propelling jet. This is observed in Fig. 6C, where the x-velocity has a very strong spatial correlation with the vortex rings. Movies 1 and 2 in supplementary material show the temporal behavior of the central jet. Such vortex ring structure was observed in experiments by Drucker and Lauder for pectoral fins of black surfperch and bluegill sunfish,where the rings are generated by cyclic flapping motion of the pectoral fins(Drucker and Lauder, 1999). Although the rings observed in our fin simulations seem similar, they are generated by a different motion pattern – a traveling wave that is perpendicular to the flapping velocity of the fin. This suggests that vortex rings can be observed along the body surface in different modes of swimming,when a predominantly axial jet is present along the body. This ring structure is distinct from the `2-S' vortex structures observed in eels(Tytell and Lauder, 2004) or the `2-P' structures (i.e. two counter-rotating vortex pairs, one on each side of the body – left and right with reference to the swimming direction)computed in eel simulations (Kern and Koumoutsakos, 2006). The difference between the current ribbon-fin flow field and the eel wakes is that the ribbon fin generates a prominent backward central jet along the length of the body.

In general, for a swimming animal in the Eulerian regime (i.e. Re1, wherein vortex shedding is expected to occur), the dominant vortex structure depends on the morphology and the swimming mode. For instance, in the case of eels,experiments indicate that 2-S vortex tubes are the dominant structures(Tytell and Lauder, 2004)whereas numerical simulations show that eels may also shed 2-P vortex rings under certain circumstances (Kern and Koumoutsakos, 2006). For gymnotiform (ventral anal fin) or ammiform (dorsal fin) swimming using a ribbon fin, our results indicate that vortex rings and the corresponding central jet are the dominant mechanism of surge force generation.

Fig. 6A also indicates that there are secondary, smaller vortex rings that are not attached to the fin but are shed into the surrounding fluid. Associated with them are the smaller jets located at the centers of these secondary rings (henceforth referred to as the secondary jets, as opposed to the primary central jet). Both the primary and secondary jets exhibit some angular deflections from the surge direction. This means that some momentum is lost in the lateral direction. Some comments can be made regarding the biological implications of this observation. For animals needing rapid maneuverability in order to capture prey (e.g. black ghost knifefish), the ability to exchange lateral momentum with the fluid is essential. The morphological characteristics required for this purpose can entail an undesirable energy loss in the lateral direction even in the cruising mode. This is also observed in numerical simulations of eel swimming with a 2-P vortex structure (Kern and Koumoutsakos, 2006), where the wake has an angular shape,reflecting energy loss in the lateral direction. By contrast, for low-maneuverability jellyfish, the jet is almost unidirectional in the surge direction (Dabiri et al.,2005a; Dabiri et al.,2005b). This would minimize the lateral energy loss but, as a consequence, the ability to generate rapid lateral forces would be diminished.

Heave generation can be understood from the streamwise vorticity(ω_x) shown in Fig. 8. The opposite surfaces of the fin produce opposite-signed vorticity. It then separates from the fin surface, rolls into streamwise vortex tubes at the bottom edge of the fin and is advected in the lateral direction (z direction in Fig. 1) by the fin motion. The induced velocity of such axial vortex pairs is vertically downward, imparting a downward momentum to the fluid. The resultant upward reaction on the fin causes the heave force. The downward jet from these rolls is depicted in the vertical velocity contours (Fig. 9).

The mechanisms of surge and heave generation that are inferred from Figs 6 and 8 are schematically depicted in terms of the streamwise jet and attached vortex tubes in Fig. 14A and B, respectively. The downward motion of the vortex tubes was observed only when the heave generation is significant, i.e. above θ_max=20 deg. Below this angle, the fin was found not to generate any significant heave(Fig. 12B). Flow atθ _max=20 deg. showed that these tubes are attached to the lower end of the fin, with little or no downward shedding, indicating they are weak compared with those in the higher heave case.

The total forces on the fin are shown in Fig. 11. The trends shown by these forces will be discussed in the next section.

Our computational results show that the surge force from the simulated ribbon fin varies between 0 and 1.85 mN over the parameter range examined. According to the measurements by Blake(Blake, 1983), a frequency of 3Hz with θ_max of 30 deg. and a specific wavelength of approximately 0.5 is typical of gymnotids. In our results, this regime results in a force of 1 mN (Fig. 12A). We compare this result with three empirical findings: estimated thrust needed during two kinds of body reversals and the measured drag force on the body. The most relevant comparison for our simulations, occurring with a non-translating fin in stationary fluid, is with body reversals where the body is momentarily non-translating.

Our estimate of the magnitude of surge available during the rapid body reversals typical of prey capture behavior, based on the 7 Hz traveling wave frequency that was measured (see Results), is ≈6 mN. Our estimate of the necessary surge force to achieve the measured accelerations, using an estimate of effective mass, is ≈13 mN. We hypothesize that the gap between these estimates is due to the presence, during the prey capture rapid reversals, of 1–3 rapid bilateral pectoral fin strokes. It should also be noted that our simulations do not account for the acceleration of the surrounding fluid during traveling wave reversals because we are assuming a stationary fluid.

During the refuge tracking behavior body reversals, no such bilateral pectoral fin strokes were observed. In this case, the estimated available surge force is ≈2 mN whereas the estimated necessary surge force based on measured acceleration and effective mass is ≈1 mN. In this case, the measured fin frequency was 4 Hz. This compares favorably with our computed thrust for similar kinematics.

Finally, the body tow drag measure enables us to estimate the necessary thrust capacity of the fin under steady swimming conditions. Our drag measurement, at typical steady-state swimming velocities, was ≈1–2 mN(Fig. 13), similar to our computed surge force. Prior theoretical estimates of the thrust of a similar ribbon-finned fish, Xenomystus nigri, are also in this range (gray line of Fig. 13)(Blake, 1983).

Fig. 12A shows the mean surge and heave forces on the ribbon fin for different frequencies. It shows that both surge and heave forces have a quadratic dependence with frequency. As the velocity scale is V=fλ, the force varies asρ V²A, where A is the characteristic area of the fin, which can be interpreted as the total wetted area. This indicates that the force generation is essentially inertial in nature at these moderate Re. Nevertheless, viscous effects may not be negligible as they change the flow in the vicinity of the fin substantially so as to affect the overall flow dynamics.

The fin mean surge and heave forces as a function of θ_maxare shown in Fig. 12B. In this case, the surge force follows a power law curve, whereas heave force becomes significant beyond around θ_max>20 deg. Similar trends have been shown in an oscillating airfoil study(Koochesfahani, 1989). In this work, Koochesfahani shows that an airfoil could produce drag (negative heave)or thrust (positive heave) depending on the frequency and amplitude of oscillation. In addition, he shows that the frequency at which the airfoil produces thrust depends on the amplitude of oscillation. Moreover, in previous work (see Blake, 1983), it has been reported that gymnotiform swimmers do not substantially change wave amplitude with swimming speed. Based on the results of the present study, we speculate that this may be to ensure an upward heave force. The empirically observed wave amplitude from a black ghost knifefish during swimming is approximately 30 deg. (Blake,1983) above the 20 deg. threshold found above.

In the numerical simulations for Set 3, the height of the ribbon fin was varied while keeping the rest of the parameters fixed. As in all the previous cases, surge force is greater than the heave force and both forces follow a power law relationship with fin height(Fig. 12C). For the shortest h_fin, the heave force obtained was insignificant compared with the surge force, and for all other fin heights the heave force was positive.

In the numerical simulations used to study the effect of specific wavelength, f and θ_max were kept constant. Increasing the wavelength produces two competing effects in force production:(1) an increase in the wave speed (V=fλ) and (2) a decrease in the total wetted area of the fin caused by a reduced number of waves on the fin (L/λ). Whereas increasing wave speed produces higher force, a decrease in the wetted area of the fin produces smaller forces. The interaction between these competing effects may be the basis of the peak in the surge force with specific wavelength shown in Fig. 12D.

For the range of traveling wave parameters that we considered, the optimal specific wavelength lies between 0.5 and 0.8. This is close to the empirically observed value of 0.4 in gymnotids (Blake,1983). It should be noted that this optimal range is only with respect to maximization of surge force. It does not consider hydrodynamic energy loss due to dissipation or turbulence or physiological energy consumption in the muscles/appendages that generate the fin deformation. Our current goal is to consider the hydrodynamic effects only and relate the force generation mechanism both qualitatively and quantitatively to various kinematic parameters of the fin motion.

At λ/L_fin=1.5, and other variables set to the baseline case, the surge force component is equaled (and for larger values ofλ/L_fin, surpassed) by the heave force component. This crossover point marks a distinct shift in the mode of propulsion. Below this critical value, the motion of the fin is predominantly undulatory and surge force dominates whereas above this value the motion of the fin is predominantly flapping (or oscillatory) and heave force dominates. The transition point between the undulatory and flapping propulsion modes will vary with the values of the other kinematic parameters.

Although we have not observed λ/L_fin near the crossover point in normal swimming of knifefish, we have observed nearly vertical translations during rare startle events (M.A.M. unpublished observations). Our results indicate that this motion may be accomplished with large values of λ/L_fin.

During forward cruising, the sum of the force components from a fish's propulsors should ideally be parallel to the body axis. Given the relationship between surge and heave as a function of fin deformation kinematics found in the present study, we can ask what the resultant angle is and compare this with the insertion angle of the fin. For the case of f=4 Hz,λ/L_fin=0.5 and θ_max=30 deg., the angle between the base (i.e. the upper edge) of the ribbon fin and the resultant force is approximately 11 deg. This is indeed similar to the angle of insertion of the fin with respect to the spine that we observe in Apteronotus albifrons (Fig. 1A).

The surge–heave crossover point may also provide insight into the position and orientation of short median paired fins in species, such as the triggerfish (family Balistidae). As the fin becomes shorter, for a given fin ray density, the number of fin rays is reduced. This results in large values of λ/L_fin. Surge and heave forces will then be comparable in magnitude. In this regime, to maximize translational thrust, the fin should have an angle of insertion similar to 45 deg. This is the case for some species of triggerfish with short median paired fins [see fig. 1 in Lighthill and Blake (Lighthill and Blake,1990)].

Scaling analysis was conducted to understand the relationship between surge force generation and the traveling wave variables. For the force scaling with f, θ_max and h_fin/L_fin, we assumed that the surge force and the independent variable have a power–law relationship of the form g(χ)=aχ^b, where g(χ) is the computed force, χ is one of the independent variables listed above and a and b are constants. aand b were found, respectively, from the intercept and the slope of a linear regression fit of the log-transformed data. For specific wavelength variation, the force did not follow a power law (see Fig. 12D), most likely due to the two competing effects contributing to the force behavior as explained above in `Effect of specific wavelength'. Hence, forλ/L_fin we used a curve fit of another form(described below).

In a series of papers, Lighthill and Blake(Lighthill and Blake, 1990)and Lighthill (Lighthill,1990a; Lighthill,1990b; Lighthill,1990c) analyzed the fluid dynamics and force production for ribbon-fin swimmers. Lighthill and Blake did a theoretical analysis of a ribbon fin to obtain an analytical expression for the fin's mean propulsive force assuming a two-dimensional irrotational flow in the plane transverse to the fin surface. The expression for the mean propulsive force consists of two terms: (1) the rate of backward momentum transport and (2) the pressure force.

The rate of backward momentum transport is proportional to Umvcalculated at the trailing edge of the fin, where the overbar denotes a mean over one period of the traveling wave on the fin, U is the swimming velocity, m is the added mass per unit length of an infinitely thin fin and v is the lateral fin velocity component that is perpendicular to the fin surface. Because in our simulations the fin is stationary(U=0), this term vanishes. This is exactly what will occur in a swimming motion starting from rest, or during the switch from forward to reverse swimming direction, a stereotypical maneuver for a black ghost knifefish during prey capture (MacIver et al., 2001). In these cases, there are instants when the force from the momentum transport term will tend to zero and the propulsive force generated by the pressure term will be dominant.

The frequency scaling from simulations and Eqn 11 is f². This contribution comes from the v² dependence of the pressure force as discussed by Lighthill and Blake (Lighthill and Blake,1990).

The scaling for the fin height parameter, h_fin/L_fin, is also in good agreement and is close to the fourth power. The h_fin scaling arises from two distinct factors in the pressure force: (1) the added mass of the fin, proportional to

For θ_max, our results show a scaling power of 3.5 vs 1.9 of Lighthill and Blake(Lighthill and Blake, 1990). Note that θ²_max scaling arises from the v² contribution to the pressure force. The difference between this scaling and the 1.9 scaling shown in Fig. 15B (gray line) is due to the effect of the fin twisting term (cos³α_av). Our simulations give an extra factor of θ^1.5_max, which we speculate more accurately resolves the effect of the fin twist onθ _max scaling. This could also be because of the local two-dimensional assumption by LB, which does not capture the full three-dimensionality of the flow-field around the fin.

For the specific wavelength variation, the estimation of the propulsive force with Lighthill and Blake's theory (Eqn 11) shows an initial rapid increase in surge force with specific wavelength followed by an asymptotic value at large specific wavelengths (see Fig. 15D). However, our computational results show a decrease in the propulsive force after a certain specific wavelength, contrary to the results of Lighthill and Blake(Lighthill and Blake, 1990). Such a peak in the surge force has been observed experimentally for a robotic ribbon fin mimicking a gymnotiform swimmer(Epstein et al., 2006). We attribute this mismatch between our results and those of LB to the fact that their inviscid theory does not capture the effect of decreased thrust due to decreased wetted area for increasing wavelength.

For very small angular deflections of the fin, the analytical estimates of Lighthill and Blake (Lighthill and Blake,1990) are larger than our computed propulsive forces. But for the majority of the parameter space, Lighthill and Blake's(Lighthill and Blake, 1990)forces are lower than our results by as much as a factor of five(Fig. 15). There are several simplifying assumptions in LB's theory that may account for this mismatch: (1)a lack of hydrodynamic interaction between different parts of the fin along the surge direction, (2) the two-dimensional flow approximation and (3) the inviscid flow assumption. Our results show a clear flow interaction and axial propagation of vorticity along the ribbon fin (see Fig. 8). This interaction will affect the force-generation capability of the ribbon fin.

Furthermore, it is possible that in start-up motions, or when the swimming speed is very small, viscous forces may contribute to thrust, something that is neglected in the analytical solution of the force estimate, which relies on the inviscid flow assumption.

According to Lighthill and Blake's theory, in cruise mode, the pressure term is small compared with the overall thrust produced(Lighthill and Blake, 1990). However, in regimes with U∼0 (e.g. start-up motions and during rapid forward-to-reverse maneuvers), the propulsive force is dominated by the pressure contribution. Our computations show that this pressure contribution is large enough to provide sufficient thrust for refuge tracking body reversals and can play a significant role in pectoral-fin enhanced reversals in prey capture strikes. This suggests that real swimmers could in fact use this mechanism during impulsive motion.

The goal of this study was to understand the mechanical basis of propulsive force generation by the anal ribbon fin of a black ghost knifefish during impulsive starting movements. Using computational fluid dynamics and DPIV, we studied the key flow mechanisms underlying surge and heave force generation by an undulating but non-translating ribbon fin. We found that the mechanism of thrust generation is a streamwise central jet associated with a crab-like vortex ring structure attached to the fin. Heave force is generated for certain wave parameters by vortex tubes attached to the tip of the ribbon fin. The simulated surge force shows good agreement with the estimates for refuge tracking body reversals, as well as our drag measurements from a cast model of the knifefish. The flow features from simulations reproduce those from DPIV. Our results also provide quantitative scaling of the surge force with amplitude, frequency and wavelength of the traveling wave, and the fin height. The present high-fidelity simulations account for many of the fluid dynamical aspects not captured in the inviscid theory of Lighthill and Blake(Lighthill and Blake, 1990)and provide force mechanisms based on flow structure, qualitative trends in forces and better quantitative accuracy desirable for subsequent modeling and engineering efforts.

LIST OF ABBREVIATIONS

 
• a, b

  constants in the power–law equation

 
• A

  fin area

 
• C₁

  constant

 
• d

  flat plate length

 
• f

  frequency of the fin

 
• F_surge

  surge force

 
• h_fin

  fin height

 
• i, j, k

  unit vectors in the x, y and zdirections

 
• L_fin

  fin length

 
• L_x, L_y, L_z

  computational domain size in x, y and z

 
• m

  added mass

 
• N

  number of trials

 
• N_p

  number of Lagrangian particles per grid cell

 
• N_x, N_y, N_z

  number of grid points in x, y and z

 
• p

  pressure

 
• r_m

  vector from the rotational axis to a point on the fin

 
• r_m

  distance of a point on the fin from the axis of rotation

 
• Re

  Reynolds number

 
• s

  vortex length

 
• t

  time

 
• u

  velocity vector

 
• u_fin

  velocity at any point on the fin

 
• û

  intermediate velocity field

 
• ũ

  intermediate velocity solution after velocity correction at the fin location

 
• u^n–1

  velocity field at previous time step

 
• U

  constant translational velocity of the plate

 
• v

  fin velocity perpendicular to the fin surface

 
• V

  wave speed

 
• x, y, z

  Cartesian coordinate directions

 
• αav

  mean angle of a fin section with respect to the sagittal (vertical)plane

 
• Δt

  time step

 
• ζ

  lack of enhancement coefficient

 
• θ

  angular position of any point on the fin

 
• θmax

  maximum angular deflection of the fin from the midsagittal plane

 
• λ

  wavelength

 
• μ

  dynamic viscosity

 
• ν

  kinematic viscosity

 
• χ

  independent variable

 
• ρ

  fluid density

 
• ϕ

  function that approximates the variation of surge force with respect toλ/L_fin

 
• ωx

  vorticity in the axial direction

 
• \({\nabla}{\cdot}\)

  divergence operator

N.A.P. led the computational and theoretical fluid analysis effort; M.A.M. led the experimental fluids, robotics and all other efforts. O.M.C. and A.A.S. contributed equally. A.A.S. wrote the computational fluid dynamics code. We thank Richard Lueptow and Nick Pohlman for the DPIV setup and help with data collection and analysis. The drag measurements were made by M.A.M. with Ann Marie Polsenberg, using Joel Burdick's and Mory Gharib's facilities at Caltech, with Matt Ringuette's generous assistance. We thank two anonymous reviewers for their very helpful suggestions. This work was supported by NSF grant IOB-0517683 to M.A.M. Partial support for N.A.P and A.A.S. was provided by NSF CAREER CTS-0134546 to N.A.P. Computational resources were provided by the San Diego Supercomputer Center through NSF's TeraGrid project grant CTS-070056T to A.A.S., N.A.P. and M.A.M., and in part by a developmental grant from the Argonne National Laboratory to N.A.P.