Undulating fins produce off-axis thrust and flow structures

A wide variety of fish undulate elongated fins to generate thrust. Knifefish undulate either an elongated dorsal fin, as is the case of the African knifefish Gymnarchus niloticus, or an anal fin, as is the case of species in the order Gymnotiformes and in the family Notopteridae. Triggerfish (family Balastidae) undulate both dorsal and ventral fins. Also, many species such as cuttlefish, skates and rays undulate elongated fins along their lateral margins. Undulatory motion, in contrast to oscillatory motion, implies that these fins produce a wave that travels along the fin with at least one full undulation present along the fin. These fish propel primarily by median and/or paired fin (MPF) propulsion, in contrast to body and/or caudal fin (BCF) propulsion where the fish propels by the actuation of the whole body, part of the body or only the caudal fin (for a review, see Sfakiotakis et al., 1999).

The biological model system used for this research is the elongated anal fin (hereafter termed the ribbon fin) of the weakly electric knifefish, such as the fin present on the black ghost knifefish Apteronotus albifrons (Linnaeus 1766) (Fig. 1A). These knifefish are important model systems for studying sensory neurobiology (Turner et al., 1999; Krahe and Fortune, 2013) and sensory motor integration (Snyder et al., 2007; MacIver et al., 2010; Cowan and Fortune, 2007). Furthermore, research on the unique biomechanics and sensory modality of electric knifefish has inspired the design of numerous underwater robots capable of agile movement and omnidirectional close-range sensing (reviewed by Neveln et al., 2013).

Blake pioneered the research of ribbon fin propulsion in his study of knifefish (Blake, 1983), in which he described the kinematics of forward swimming fish. Subsequent work has quantified body kinematics (MacIver et al., 2001) as well as fin kinematics (Ruiz-Torres et al., 2013). Lighthill and Blake first modeled the hydrodynamics of the undulating ribbon fin using elongated-body theory (Lighthill and Blake, 1990). Work in our laboratory has examined the hydrodynamics of a non-translating undulating fin (Shirgaonkar et al., 2008), as well as measured the forces and free-swimming velocities generated by a biomimetic robotic ribbon fin, shown in Fig. 1B (Curet et al., 2011b). Results show that a single traveling wave along the fin is capable of producing both surge and heave components of thrust (see Fig. 1C for axes definitions). Knifefish can also partition their ribbon fin into two counter-propagating waves that originate on either end of the fin (Curet et al., 2011a; Ruiz-Torres et al., 2013; Sefati et al, 2013). Knifefish control the position of the point where the two waves meet along the fin to hover in still water and swim at slow speeds (Ruiz-Torres et al., 2013; Sefati et al., 2013). These counter-propagating waves on a single ribbon fin can be used to vary the proportion of surge and heave components of thrust (Curet et al., 2011a; Sefati et al., 2013).

Because of the unique properties of the propulsive mechanisms of the ribbon fin, weakly electric knifefish exhibit high maneuverability, including backward, upward and rolling motions. Even for time periods of the order of the sensorimotor neuronal delay time (~100 ms), they can maneuver to any location within their omnidirectional sensing range for small prey (Snyder et al., 2007). This tight coupling of sensory, neuronal and mechanical systems presents an interesting case study in the co-evolution of these complex, interdependent systems and the inherent trade-offs that occur when trying to maximize a relevant objective function such as the ratio of energy intake to expended energy (MacIver et al., 2010). A better understanding of the intricate mechanical properties of the ribbon fin is essential for further decoding the interdependency of sensory, neuronal and mechanical systems in this model system.

This study aims to further our understanding of hydrodynamic principles underlying cruising ribbon-fin swimmers through experimentation on a robotic model using digital particle image velocimetry (DPIV). We also use high-fidelity computational fluid dynamics of a virtual model of our robot and an accurate knifefish model with measured body and fin morphology combined with measured ribbon fin kinematics. We first describe the flow features present in the wake that are responsible for generating thrust along both the surge and heave axes of the fin. We then show the consistency of this flow structure across a range of ribbon fin sizes and kinematics with Reynolds numbers spanning over one order of magnitude. Finally, we introduce a schematic of the ribbon-fin wake vortex structure and compare MPF undulatory swimming to well-known BCF swimming mechanics under a framework for comparing flow structures using non-dimensional parameters such as wave efficiency and the Strouhal number.

The free-swimming pitched robot achieved a swimming speed of 26 cm s⁻¹ (0.8L_fin s⁻¹, where L_fin is fin length) for the kinematic parameters tested [frequency (f)=3 Hz, number of undulations=2, maximum angular amplitude (θ_max)=30 deg]. Simulation FishSim1, which used the mean kinematic parameters of an actual free-swimming A. albifrons [14.9 cm s⁻¹ data (from Ruiz-Torres et al., 2013)], obtained a steady-state swimming speed of 12 cm s⁻¹ (1.33L_fin s⁻¹). This speed underestimates the actual swimming speed of the fish by ~20%. Simulation FishSim2, which is the same morphology and kinematics as FishSim1 except that the amplitude envelope was constant at 30 deg, obtained a steady-state swimming speed of 17.7 cm s⁻¹ (1.97L_fin s⁻¹). The simulation of the robot, GhostbotSim, obtained a steady-state swimming speed of 26 cm s⁻¹ (0.80L_fin s⁻¹). The swimming speed matched the experimental robotic results very closely, providing validation of the simulation at the higher Reynolds number of the robot (Re=80,000) over that of the fish (Re=13,000). Finally, FinSim, which simulated an undulating rectangular fin of the same length as the FishSim fins but without a body, swam 9 cm s⁻¹ (0.89L_fin s⁻¹). These results, as well as important kinematic parameters and non-dimensional numbers such as wave efficiency and Strouhal number, are reported in Table 1. More details about the simulations can be found in Materials and methods.

DPIV data were collected by imaging seven different planes of the swimming robot while submerged in a flow tunnel (Fig. 2). More details of the data collection and analysis can be found in Materials and methods.

Fig. 3 shows both snapshots of the velocity vector field with corresponding vorticity as well as the average jet generated by the fin in the midsagittal plane M (Fig. 2B). In a particular time instant (Fig. 3A), velocity vectors (which have the incoming flow velocity subtracted) have varying directions along the length of the fin, alternating from more aligned in the surge direction to more aligned in the heave direction. A line of clockwise-rotating vortices are shed at an angle to the fin similar to the pitch angle of the fin (horizontal line of blue areas that intersects the front of the fin). Behind the fin, there are areas of vorticity with alternating sign. Fig. 3B shows the phase-averaged velocity vector field at the same instant in the cycle as Fig. 3A. Supplementary material Movie 1 shows the time evolution of the velocity vector field with associated vorticity over one cycle.

On average, the jet generated by the fin is aligned with the flow direction, as expected if thrust is aligned with swimming direction (Fig. 3C). This jet is strongest and deepest just downstream of the caudal end of the fin. The depth of the jet increases approximately linearly from the rostral edge of the fin to the caudal edge, and tapers off further downstream. There is a smaller secondary jet just downstream of the fin aligned along the surge axis of the fin.

Fig. 4 shows four snapshots of the phase-averaged velocity vector field and associated vorticity for the horizontal plane H (Fig. 2B). These snapshots split the fin cycle into quarters. Velocity vectors (which have the incoming flow subtracted) in the area of the plane transecting the fin are oriented largely laterally. Just downstream, vortices are shed off of the fin into the plane. These vortices alternate in rotational direction and are shed twice per fin cycle. Fig. 4D outlines these vortices, where the area encompassed by a loop corresponds to the strength of the associated vortex. The velocity vectors associated with these vortices have both downstream and lateral components. Further downstream, the strength of the vorticity decreases, and the velocity vectors orient more along the direction of the flow. Supplementary material Movie 2 shows the time evolution of the velocity vector field with associated vorticity over one cycle.

Fig. 5 shows four snapshots of the velocity vector field associated with the five transverse planes T1–T5. As in Fig. 4, the snapshots split the fin cycle into quarters. Because there are two undulations along the fin, the fin rays in T1, T3 and T5 are all roughly in phase with each other while in antiphase with the fin rays in T2 and T4. Supplementary material Movie 3 shows the time evolution of these planes over one cycle.

The velocity vectors in these transverse planes have the smallest magnitude in the anterior T1 plane, with magnitudes increasing moving posteriorly to T5. Vortices of alternating direction are present in each plane upon reversal of the fin ray, as indicated in Fig. 5D. These vortices are weakest in T1 and increase in magnitude moving posteriorly to T5. Consequently, the flow through these vortices is strongest in T5, while non-existent in T1. These flows are oriented mostly downward from the fin, with smaller lateral components that oscillate from left to right over a cycle. Also, vortices remain closely associated with the fin edge in T1 as the fin moves laterally, rather than being shed downwards off the fin edge as is evident in T5.

To visualize the 3D wake structure of the simulations, we visualized iso-surfaces of the q-criterion (Dubief and Delcayre, 2000). The q-criterion subtracts the strain rate from the rotation rate of a velocity vector field, so that regions of positive q indicate vortical structures. Fig. 6 shows the bottom and side views of these iso-q surfaces for the four simulations. Vortical structures that are consistently present across the three simulations are indicated by the magenta arrows. These structures are vortex tubes that are initially shed by the rostral half-wave of the fin and then continually shed at the peaks of the waveform as the wave travels caudally. Supplementary material Movies 4–6 show the evolution of these structures over time. As the shed vortex tubes are left by the moving fin, they become entangled downstream to form a complex wake structure, until the strength of the vorticity dies out.

The fourth simulation, labeled FinSim, gives the clearest view of the wake structure. Supplementary material Movie 7 shows the evolution of the iso-q surfaces, and a snapshot view is included in Fig. 6D for comparison with the other simulations. The flow structures in FinSim are less entangled than those in FishSim1 and FishSim2, likely because the traveling wave frequency is less by more than a factor of three (Table 1), giving a clearer view of the wake. However, the manner in which these vortex tubes are shed is similar to the shedding of vortex tubes from the other simulations.

Both the experimental and computational results indicate that the wake structure across both fish and robot body and fin morphologies maintain the same qualitative features. First, we discuss how exactly the flow features of the 2D DPIV planes map to the overall 3D vortical structure of the wake. Second, we compare the 3D wake of free-swimming ribbon fins with prior results from a non-translating undulating fin (Shirgaonkar et al., 2008). Third, we introduce important relationships between the kinematic parameters of the undulating fin and non-dimensional numbers associated with free swimming such as wave efficiency (sometimes referred to as slip) and Strouhal number. Fourth, we relate these non-dimensional numbers to the wake structure of the undulating ribbon fin. Last, we compare the wake of the ribbon fin with the well-characterized wake of carangiform swimming, showing specifically how knifefish have maintained similar thrust-producing mechanisms with a fin of drastically different morphology and kinematics.

There has been some effort to show how experimental 2D DPIV data correlate with hypothesized 3D wake structures in various forms of swimming. For example, the wake of a carangiform swimmer with an oscillating caudal tail is thought to generate a series of connected vortex rings (Nauen and Lauder, 2002). A horizontal slice through these vortex rings illuminates the classic reverse von Kàrmàn vortex street visible in planar DPIV measurements (Nauen and Lauder, 2002). Similar reverse von Kàrmàn vortex streets are present in the transverse planes measured for the Ghostbot fin, especially in T3–T5 (Fig. 5). The vortices in the transverse planes must be sections of a 3D vortex structure. Indeed, a vortex tube is shed as the wave travels along the length of the fin, as seen in the iso-q surfaces of the simulations (Fig. 6, supplementary material Movies 4–7).

The vortex tubes shed from the fin edge also account for the other vortices present in the horizontal and midsagittal planes. In the horizontal plane, strong vortices are present at the posterior point of the intersection between the fin and the plane, indicated by the solid white line in Fig. 4. These vortices are lengthwise slices through the vortex tube as it is being shed, orthogonal to the slices in the transverse planes. In the midsagittal plane, a series of similarly rotating vortices are shed in a line at an angle to the fin closely matching the pitch angle of the robot (Fig. 3). Again, 2D vortices are slices through the vortex tubes shed off the fin. As seen in Fig. 6, the vortex tubes bend medially towards the opposing side and cross the midsagittal plane. Successive vortex tubes have opposite rotation in the transverse plane (Fig. 5), but are also shed on opposite sides of the fin. Therefore, when the vortex tubes bend medially and intersect the midsagittal plane, every vortex tube has the same rotational direction such that higher downstream flow (indicated by the light blue color) is closer to the fin. Over the course of a cycle, this high downstream flow averages to the jet visible in Fig. 3C.

The shedding and distortion of the fin edge vortex tubes are schematized in Fig. 7. Vortex tubes are color-coded based on the side of the fin from which they are shed. This schematic encapsulates many of the main features present in the vortex tubes shed by the fin edge. Towards the rostral portion of the fin, where the tube is just starting to be shed, the tube is oriented laterally due to being shed at a single time by the rostral half wave of the fin (Fig. 7A). Each vortex tube should be connected to the previous tube as well as the subsequent tube according to Helmholtz's second theorem of vorticity (Kuethe and Schetzer, 1950). The schematic details this connection for the anterior part of these tubes, though they should also connect on the posterior ends of the tubes after being shed off the back of the fin. In between, the vortex tube remains roughly in the location where that portion of the tube was shed (Fig. 7B). Therefore, the downstream spacing of subsequent vortex tubes is mostly determined by the swimming speed of the fin. We will return to the issue of vortex spacing in further discussion below.

The wake structure remains consistent throughout all simulations, which range in Reynolds number over approximately one order of magnitude. Fig. 6 highlights the consistent vortex structures when viewed from below. Furthermore, Fig. 8 compares the midsagittal plane of the schematics with the DPIV data as well as all four simulations. In Fig. 8B–E, the line of blue vortices emanating from the fin along the direction of swimming is consistent with the lateral reorientation of the edge vortex tubes. Fig. 8E further shows the connection between the vorticity in the midsagittal plane and the 3D vortex structure. Whether robot scale, fish scale or without any body, the wake structure remains qualitatively similar.

The hydrodynamics of elongated ribbon fin propulsion is highly complex with significant 3D effects. Previously, we investigated the hydrodynamics of a non-translating undulating ribbon fin in stationary water (Shirgaonkar et al., 2008). This flow condition was used to explore the hydrodynamics of impulsive motion. Impulsive motion is common in knifefish as they are constantly changing their direction of motion during prey capture (MacIver et al., 2001).

The non-translating undulating fin produced both surge and heave forces, just as the average velocity of a freely swimming fin has both surge and heave components. We attributed the heave forces in the non-translating case to vortex tubes generated by the fin edge as it oscillates back and forth laterally. Similarly, Fig. 5 shows how the vorticity generated and shed by the fin edge creates a flow with large heave components. We attributed the surge forces in the non-translating case to the longitudinal jet created by the motion of the traveling wave. We inferred that the primary vortex rings encircling the jet were induced by the flow and were separate from the vortex tubes generated by the fin edge. However, in the current free-swimming simulations as well as the DPIV data, we see that the vortex tubes shed by the fin become distorted and actually wrap around the overall jet shown in Fig. 3. Consequently, the vortex rings and the fin edge vortex tubes that were previously thought to be separate phenomena associated with surge and heave forces, respectively, are likely different parts of the same vortex structures shed by the fin. A comparison of these vortex rings of the current FinSim and the non-translating fin from Shirgaonkar et al. (Shirgaonkar et al., 2008) is shown in Fig. 9.

St can be considered to indicate the ratio of local inertial forces to convective inertial forces, or in the case of a swimming fish, the ratio of local momentum input to the wake by the fin to the output momentum imparted on the fish by the wake, which moves the fish forward. Therefore, St can be considered a second measure of efficiency, where the fin aims to achieve the highest swimming velocity while minimizing lateral fin velocity. When St is considered as a type of efficiency, clearly a lower St implies better performance, with St=0 indicating perfect performance. The natural range of 0.2 to 0.4 for St indicates that there exists a practical lower limit of around 0.2 (gliding is a case where St=0, though the gliding velocity is not constant over time).

It should be noted that a better measure of efficiency for a cruising animal or robot would relate the total input power (P) to the system to the useful output power, such as propulsive efficiency (η_P=FU/P, where F is thrust or drag) or cost of transport (COT=P/U). For example, Anderson et al. found that η_P was maximal for St in the range of 0.25 to 0.40 for an oscillating foil with varying kinematics (Anderson et al., 1998). Moreover, Moored et al. found a local maximum of η_P when the driving frequency of an oscillating fin matched the hydrodynamic resonant frequency (Moored et al., 2012). However, these oscillating foils were not self-propelled. While we did not measure input power to the system, we would expect the power needed to drive the rays of either the robot or the fish to increase with ray velocity, which increases with both frequency and amplitude of oscillation, parameters in the numerator of St. Therefore, we expect COT to follow similar trends as St, and that the lowest St values observed are approaching the best propulsive efficiencies achievable.

From both of these relationships, it is clear that higher wave efficiencies and lower Strouhal numbers result in larger vortex spacing. Higher vortex spacing indicates that the fin is obtaining more impulse per vortex tube, which would indicate greater efficiency. It is hypothesized from the well-studied carangiform wake, however, that a vortex spacing resulting in a reverse von Kàrmàn street with associated linked small core vortex rings is the signature of highly propulsive thrust (Nauen and Lauder, 2002). These linked rings seem to limit the vortex spacing, which in turn limits St. For example, when St=0.2, S=2.5A, so for vortex spacing greater than this value the linked vortex pattern likely becomes unlinked and breaks down. Further discussion on predicting St from measured distances associated with the wake for ribbon fins as well as other forms of swimming is provided in a later section.

There have been a few studies relating the wake structures to St in both anguilliform and carangiform forms of swimming. A series of articles by Borazjani and Sotiropoulos show in simulation that as Reynolds number increases, St decreases for free-swimming carangiform and anguilliform swimming (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009). For low Reynolds numbers (Re=300 and 4000), St was always around 0.6 or above. The St fell within the observed range of 0.2 to 0.4 only for the inviscid simulations. For these low St scenarios, the wake was a single row wake structure consistent with the linked vortex rings associated with carangiform swimming. For higher St, the wake split into a double row wake structure. In self-propelled cases, the transition from single row to double row wake structure seems to occur as the vortex spacing decreases past a critical value, corresponding with increasing St. We expect that if the ribbon fin swam at higher St, a similar double row wake structure would develop. Other studies of heaving and pitching foils in an imposed flow show that changes in the effective angle of attack also can cause the transition from single row to double row wake structure (Blondeaux et al., 2005; Dewey et al., 2012).

For undulatory fins such as the knifefish fin, the vortex spacing also depends on wavelength and wave efficiency given by Eqn 5. This trend is visible in Fig. 8, where the spacing of the line of vortices in the midsagittal plane is determined by the wave efficiency, as the wavelength is scaled in each image to be equal. For a case with high wave efficiency, such as FinSim (Fig. 8E), the spacing is visibly larger than a case with low wave efficiency, such as FishSim1 (Fig. 8D). Clearly, if wavelength were further decreased, the vortex spacing would decrease, resulting in interference between adjacent vortex tubes. However, if wavelength was increased, the surge component of the thrust would vanish, as described in the following section.

The feature of the ribbon fin that differentiates it from most other forms of underwater propulsion is that the thrust produced has components in both the surge and heave axes. For that reason, there exists a consistent insertion angle of the fin to the body, so that the longitudinal axis of the body aligns more with the direction of thrust and the corresponding wake. The ribbon fin can even vary the proportion of the surge component of thrust to the heave component of thrust, though fish accomplish this mainly using two counter-propagating waves (Curet et al., 2011a; Sefati et al., 2013). Can the flow structure for a single traveling wave be broken into heave features and surge features? Consider the wake generated by an oscillating caudal fin of a carangiform swimmer such as a trout, as shown in Fig. 10A. As mentioned, the wake consists of a series of linked vortex rings, with an increased flow through those rings. There are two vortex rings per cycle, and each associated jet has a lateral component with alternating direction. The height of the ring is on the order of the height of the fin, while the axial spacing of the rings is determined mostly by the Strouhal number and the amplitude of oscillation.

The structures created by a trout surging can be related to the structures generated by a ribbon fin if it were to be actuated with a traveling wave of infinite wavelength. In this case, as shown in Fig. 10B, each ray of the fin is in phase with every other ray, so the fin oscillates similarly to the caudal fin of a trout. This oscillation is now centered around the heave axis of the ribbon fin instead of the surge axis of the caudal fin, and would tend to move the fish in the fin-fixed heave direction. Assuming an appreciable velocity could be reached so that St fell within the range of ~0.2 to 0.4, we predict that the wake would contain a similar trail of linked vortex rings, now with the width of the rings on the order of the length of the fin as shown in Fig. 10B.

Moving to the biological case of an undulating fin with a finite wavelength, as shown in Fig. 10C, the portions of the vortex tubes that bend up from the rostrocaudally stretched tube are no longer shed simultaneously, but over a period of time as the wave travels along the fin. This wave motion in the surge direction contributes momentum to the fluid along that axis so that a surge component of thrust emerges. The vortex rings are now realigned along a line at an angle ventral to the fin, approximately the same angle as the insertion angle of the fin to the body. The eel-like undulatory motion of the ribbon fin is responsible for this realignment of the wake structures and therefore the off-axis direction of thrust. However, the vortex ring pattern has more in common with the small core linked vortex rings found in the wake of carangiform swimmers than with the double row wake structures found in anguilliform locomotion (Tytell and Lauder, 2004; Kern and Koumoutsakos, 2006).

The wake of a fish can be considered the ‘footprint’ of the fish. We have described the footprint of an undulating fin and explored the similarities and differences compared with the footprints of other forms of underwater locomotion. The larger structure of the wake of a self-propelled fish can predict the swimming mode. First, a wake with double row vortex rings and associated lateral jets implies an anguilliform swimmer (Tytell and Lauder, 2004; Kern and Koumoutsakos, 2006). Second, a wake with a single row of linked vortex rings and associated downstream jet implies a carangiform swimmer (Nauen and Lauder, 2002). Third, a wake with a single row of diagonally distorted and linked vortex tubes implies an undulating ribbon fin swimmer. It is important to note that the transition from a single row wake to a bifurcated double row wake is due mostly to the differences in St exhibited by anguilliform and carangiform swimmers, and not simply to differences in kinematics (Tytell et al., 2010), leading us to predict that if a ribbon fin fish swam at a higher St, the wake would begin to bifurcate.

Furthermore, a closer examination of the vortex structure can produce an estimation of St. Using the equation for the vortex spacing along the swimming direction in terms of St (Eqn 6), we can predict St by measuring the distance between successive vortex tubes or rings as well as the width of the wake (which corresponds to A in Eqn 6). The lower the ratio of the wake width to the vortex spacing, the lower the St, indicating a fish that is able to swim further per ‘footprint’. These predictions can be extended to the bifurcating wake of an eel, where the distance S is the streamwise distance between two alternating vortex rings, which is indeed much smaller than that of a carangiform swimmer, hence the difference in St.

To study the flow structures of ribbon fin propulsion, we used a mechanical fin model that mimicked the undulatory kinematics of the ribbon fin. The mechanical fin was a slightly modified version of the knifefish-inspired ‘Ghostbot’ (Fig. 1B) used in a number of previous studies (Curet et al., 2011a; Curet et al., 2011b; Sefati et al., 2013). The fin remained 32.60 cm long (L_fin), but was modified to be 5.00 cm deep (h_fin), rather than the 3.37 cm in prior studies, to generate a higher amount of thrust. The Ghostbot is roughly three times the length of an adult black ghost knifefish like the one shown in Fig. 1A. The fin was actuated by 32 stainless steel rays, emulating the bony rays of the knifefish ribbon fin, each driven independently by a servo motor (Curet et al., 2011b). The fin membrane connecting these rays was lycra fabric, which had an elastic modulus of ≈0.2 MPa, similar to the collagenous membrane between fish fin rays (Lauder and Madden, 2006). The motors and other hardware were housed in a submersible cylindrical body (6.35 cm diameter) with rounded-conical caps at both ends. The total length of the robot was 60 cm.

The robot was submerged in a variable-speed flow tunnel and supported from above on frictionless air bearings (New Way Air Bearings, Aston, PA, USA), as shown in Fig. 2A. The working section of the tunnel was 100 cm long, 33 cm wide and 32 cm deep. The air bearings permit movement along the longitudinal direction of the flow tunnel; however, the platform can also be fixed to mechanical ground through a force transducer (LSB200, Futek, Irvine, CA, USA), so that longitudinal forces can be measured, similar to the methods used in Curet et al. (Curet et al., 2011b).

A DPIV system was built to illuminate 2D planes of fluid in the water tunnel. The DPIV system was modeled after that used in Schlicke et al. (Schlicke et al., 2007). Briefly, the beam from a 2 W laser (Verdi G2, Coherent Inc., Santa Clara, CA, USA) was directed to an oscillating galvanometer mirror system (6215HM40B, Cambridge Technology Inc., Bedford, MA, USA). This oscillating mirror was located at the focal point of a custom-designed parabolic mirror, so that successive instances of the beam reflected off the parabolic mirror remained parallel, creating a strobed laser-light sheet. This light sheet was then directed into the water, illuminating silver-coated glass spheres immersed in the fluid (44 μm beads, Potter Industries, Valley Forge, PA, USA). High-speed video (FastCam 1024P PCI, Photron, San Diego, CA, USA) of the illuminated plane was recorded at the same rate that the laser beam was scanned, so that each frame recorded one transversal of the laser beam. In all experiments, the image capture rate was 500 frames s⁻¹. The width of the laser sheet was ~30 cm.

DPIV data were analyzed using a commercial software package (DaVis, LaVision GMBH, Göttingen, Germany). Analysis consisted of a calibration using an evenly spaced grid, video data importation and cross-correlation for velocity vector calculation. Two passes were performed for each pair of images, decreasing from an interrogation window of 64×64 pixels with 50% overlap on the first pass to 32×32 pixels with 50% overlap on the second pass. Subsequent analysis of the velocity vectors was performed in MATLAB (The MathWorks, Natick, MA, USA), as described below.

While Curet et al. (Curet et al., 2011b) only measured free-swimming velocities in the surge direction, they noticed that the fin produces both surge and heave force components for single-traveling wave kinematics. Therefore, for a better measure of free-swimming performance, the fin must be pitched at the angle where the overall direction of thrust is aligned with the direction of flow. To estimate this angle, we performed a series of experiments over a range of pitch angles from 9 to 19 deg in which we expected to find the correct angle. For each angle, we fixed the robot to mechanical ground through a load cell to measure longitudinal force. Our setup did not allow for measurement of the vertical force. We then varied the speed of the flow around the swimming velocity expected for the nominal kinematic fin parameters while measuring force, where negative force indicated the robot swimming below its free-swimming speed and positive force indicated swimming above its free-swimming speed. The result was a linear relationship between flow speed and force with a zero crossing at the actual free-swimming speed where thrust and drag cancel. We also collected DPIV data from the midsagittal plane, which indicated the direction of the wake generated by the fin. We found that at a pitch angle of 13 deg, the swimming speed was maximized and the wake was aligned with the flow. Thus, for collection of all subsequent DPIV data, we pitched the fin at 13 deg relative to the direction of flow and set the flow speed to the measured free-swimming speed at that pitch angle. In these experiments, the deepest point of the fin was always twice the boundary layer thickness from the bottom of the tank as estimated from the DPIV data.

The DPIV planes chosen to investigate the wake generated by the fin are shown in Fig. 2B. The DPIV system is designed to allow for sagittal and horizontal planes parallel with the flow direction, and transverse planes perpendicular to the flow direction. Midsagittal plane (M) and horizontal plane (H) are composites of multiple DPIV data sets collected separately at different areas of the plane because of the limited width of the laser sheet (original cross-correlation on the images was performed separately as well). DPIV sets from a single plane were spatially and temporally aligned. Data were linearly interpolated onto a uniform grid for display, with a vector spacing greater than the outcome of the original cross-correlation of the DPIV data. A five-frame moving average was used to temporally filter the data, an important step for the transverse planes (T1–T5) because of the noise in the data due to the fact that particles generally only stayed within the laser light sheet for a small number of frames. Lastly, the velocity vector fields were phase averaged over five full cycles (T_cycle=f⁻¹=333 ms).

In order to relate the flow structures of our robotic model back to that of the actual fish, we used fluid dynamics simulations to step from a computational model based on the body shape (MacIver and Nelson, 2000), fin shape and fin kinematics (Ruiz-Torres et al., 2013) of A. albifrons to a computational model of the Ghostbot fin shape and kinematics. The first simulation, FishSim1, most closely matches the kinematics measured for the 14.9 cm s⁻¹ case from Ruiz-Torres et al. (Ruiz-Torres et al., 2013). A sinusoidal traveling wave was implemented on the measured fin morphology with the time-averaged kinematic parameters reported in Ruiz-Torres et al. (Ruiz-Torres et al., 2013) (see Table 1 for details). The amplitude profile of FishSim1 was based on the maximum excursion of each ray of the actual fish kinematics. In FishSim2, the morphology, frequency and number of undulations remained the same as FishSim1, and only the amplitude profile was changed to a constant 30 deg, as indicated by the dark gray shaded entries in Table 1. RoboSim maintained the same constant amplitude profile of FishSim2, but the robot body and fin morphology were modeled. Also, the frequency and number of undulations were adjusted to match the actual robot case, indicated by the light gray shaded entries in Table 1. However, the fin wave speeds (fλ) of the RoboSim and the FishSims were similar (48 and 40 cm s⁻¹, respectively). Lastly, we also include FinSim, which simulates an undulating rectangular fin on the scale of the knifefish fin without a body. The details of the computational method are described below.

Our computational method uses the constraint-based immersed body (cIB) method (Bhalla et al., 2013; Sharma and Patankar, 2005; Patankar et al., 2000; Patankar, 2001; Patankar and Sharma, 2005; Glowinski et al., 1999; Shirgaonkar et al., 2009; Curet et al., 2010), which is a variant of regular immersed body (IB) method of Peskin (Peskin, 2002). The cIB method solves the combined momentum equation of the fluid and immersed body (incompressible Navier–Stokes equation) while imposing the constraint of the prescribed motion of the immersed structure in a computationally efficient manner. Specifically, the combined momentum equation is solved in the computational domain D represented on a Cartesian grid in the Eulerian frame of reference. The immersed body is modeled as a Lagrangian collection of immersed particles that are free to move on the background Eulerian grid. The immersed particles do not require any connectivity information and thus no explicit meshing is needed in the body domain, as is generally needed in the finite element type (FEM) methods. The presence of the immersed body is modeled via an additional body force in the incompressible Navier–Stokes equation that is present only inside the volume region physically occupied by the immersed body at any given time-instant. The IB method and its variants, such as the cIB method, permit the usage of fast Cartesian grid solvers, which are otherwise difficult to use with the unstructured grids employed within the FEM framework.

The outlined numerical scheme is implemented in IBAMR (IBAMR, 2013), an open-source C++ library that provides support for fluid–structure interaction models that are based on the IB method. IBAMR uses SAMRAI (SAMRAI, 2007; Hornung and Kohn, 2002; Hornung et al., 2006), PETSc (Balay et al., 2009; Balay et al., 2008; Balay et al., 1997), hypre (HYPRE, 2011; Falgout and Yang, 2002) and libMesh (libMesh, 2013; Kirk et al., 2006) libraries, among others, for its functionality.

The simulations were carried out in a 3D rectangular computational domain D. In all simulations, zero-velocity boundary conditions were imposed on the domain boundaries, except for the simulations performed for the FinSim case. For the FinSim case, periodic boundary conditions were imposed on the domain in the rostrocaudal direction, and zero-velocity boundary conditions were used for the remaining boundaries. Details of the domain size, refinement ratio, number of grid levels and coarsest grid size, and vorticity threshold used for tagging high-gradient vortex regions are reported in Table 2.

The numerical simulations were carried at the Northwestern University Quest supercomputing facility. The cluster is made up of 502 nodes, each consisting of eight 64 bit, 2.26 GHz Intel Nehalem E5520 processors with 48 Gb of memory for each node.

We thank George Lauder for his helpful discussions concerning fish swimming mechanics and DPIV analysis. We also thank James Snyder for his help in programming the Ghostbot as well as Stuart Cameron for his advice on the design of the DPIV system.