Hydrodynamics of the bluegill sunfish C-start escape response: three-dimensional simulations and comparison with experimental data

The body motion of a bluegill sunfish during a fast start is quite different than that during steady swimming. During a fast start (also called C-start), the sunfish bends its body into a C shape, which is followed by one or more alternating tail beats (Domenici and Blake, 1997; Tytell and Lauder, 2008). This behavior has been divided into two major phases: Stage 1, the initial C-bend; and Stage 2, the stroke in which the body bends out the C shape (Weihs, 1973). Such fast start behavior has direct consequences for fitness, as it is used to flee predators (Walker et al., 2005) or to capture prey (Canfield and Rose, 1993; Harper and Blake, 1991; Webb, 1984). Therefore, the fast start behavior has been studied extensively for different species from different perspectives. For example, numerous studies deal with muscle activity (Ellerby and Altringham, 2001; Jayne and Lauder, 1993; Tytell and Lauder, 2002; Westneat et al., 1998), neural control (Eaton et al., 2001; Fetcho, 1991; Koyama et al., 2011; Tytell and Lauder, 2002; Zottoli et al., 1995) and kinematics (reviewed by Domenici and Blake, 1997).

However, there are relatively few studies dedicated to understanding the hydrodynamics of the C-start, even though this is crucial for evaluating the C-start performance of different species with different body and/or fin shapes and kinematics. Throughout the evolutionary history of fishes, and, by extension, all vertebrates, there has likely been a strong selective pressure for fast start performance, because fishes use this behavior to escape from predators (Walker et al., 2005). To understand the evolution of body and fin morphology across all vertebrate species, one must examine how body form and kinematics affect the underlying mechanisms of hydrodynamic force and power generation.

To the best of our knowledge, however, only one group has directly studied the hydrodynamics of true C-start maneuvers. Tytell and Lauder studied the hydrodynamics of bluegill sunfish during the C-start using particle image velocimetry (PIV) on the mid-depth plane of the fish (Tytell and Lauder, 2008). They found that the hydrodynamics is dominated by three distinct jets induced by the body motion, which help the fish initiate and complete the turn. By calculating the momentum of each jet from the measured two-dimensional (2-D) velocity fields via a control-volume approach and summing the momentum of the three jets, Tytell and Lauder (Tytell and Lauder, 2008) provided a 2-D estimate of the force acting on the fish during the two stages of the C-start maneuver.

Several other studies have examined the hydrodynamics of turning maneuvers. The neural control of these behaviors is very different from that of a C-start, but the kinematics and hydrodynamics are similar, albeit slower. Wolfgang et al. performed pioneering simulations of turning of a giant danio using a potential flow model (Wolfgang et al., 1999). They prescribed the body kinematics from experimental observations (Anderson, 1996) and compared the computed flow field against PIV measurements. Their model yielded a simplified description of the complex three-dimensional (3-D) structure of the near-body flow but captured the large 2-D flow features on the mid-depth plane with good accuracy. Epps and Techet (Epps and Techet, 2007) employed PIV and observed two distinct vortices forming during a turn of a giant danio. They further calculated the hydrodynamic impulse acting on the fish by assuming that the wake consists of two axisymmetric vortex rings. More recently, Conte et al. tested a simplified mechanical fish that can accelerate in manner similar to a pike and found two vortex rings shed near the lateral directions when the acceleration was at its peak (Conte et al., 2010).

Although these recent studies have provided the first insights into the hydrodynamics of the C-start, the 3-D structure of the flow during this behavior has not been explored and remains poorly understood. Yet, the 3-D flow field is crucial for: (1) accurately calculating forces during turning and (2) understanding the 3-D wake structure, which could determine how well predators or conspecifics can sense the wake (Hanke and Bleckmann, 2004). In this work, we integrate the experimental work of Tytell and Lauder (Tytell and Lauder, 2008) with a high-resolution 3-D computational fluid dynamics model to study, for the first time, the 3-D hydrodynamics of the C-start of a fish. We use the experiments to reconstruct the fish body motion and C-start kinematics, which are prescribed as input to the 3-D numerical model of Borazjani and Sotiropoulos (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009; Borazjani and Sotiropoulos, 2010), to obtain the 3-D velocity and pressure fields throughout the C-start maneuver. We carry out simulations both for finite Reynolds number (Re; viscous flow simulations at Re=4000) and infinite Re (inviscid flow simulations) to compare the forces and wake structure produced in each flow regime. We compare the calculated instantaneous 2-D flow fields with the PIV measurements and show that they are in very good agreement with each other throughout the entire C-start maneuver. Further, we elucidate for the first time the 3-D structure of the vortices shed by the turning fish and calculate the forces imparted by the flow on the fish during Stages 1 and 2 of the turning process. We calculate the forces using both the full 3-D flow field as well as the 2-D flow field at the mid-body horizontal plane using the same method as Tytell and Lauder (Tytell and Lauder, 2008). This dual approach enables us to validate the numerical simulations by comparing the measured and computed 2-D force estimates and explore the adequacy of using 2-D measurements to estimate the force during C-start. To our knowledge, our work is the first to employ 3-D computational simulations informed by and validated with high-resolution laboratory experiments with live fish to illuminate the 3-D flow patterns during any animal behavior.

The paper is organized as follows. First we describe the numerical method and the method we employed to reconstruct the C-start kinematics from experimental measurements. Subsequently, we present our results and compare them to the experimental measurements. Finally, we discuss our results as they relate to previous studies and underscore the insights gained with integrating numerical simulations with experiments. We conclude by discussing the effect of viscosity and Reynolds number on wake structure and forces during the C-start.

We employed the sharp-interface immersed boundary method (Gilmanov and Sotiropoulos, 2005) to solve the incompressible, unsteady, 3-D Navier–Stokes equations in a background Cartesian grid within which the fish body is immersed. The method has been described in detail in our previous publications (Borazjani et al., 2008; Ge and Sotiropoulos, 2007; Gilmanov and Sotiropoulos, 2005) and only a very brief description of the technique is given herein. We employed an unstructured, triangular mesh to discretize and track the position of the 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 nodes were reconstructed by quadratic interpolation along the local normal to the boundary. The reconstruction method has been shown to be second-order accurate (Borazjani et al., 2008; Gilmanov and Sotiropoulos, 2005). The immersed boundary nodes at each time step were recognized using an efficient ray-tracing algorithm (Borazjani et al., 2008).

The method has been validated extensively for flows with moving boundaries (Borazjani et al., 2008; Gilmanov and Sotiropoulos, 2005) and has also been applied to simulate steady swimming of tethered (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009) and self-propelled (Borazjani and Sotiropoulos, 2010) carangiform and anguilliform swimmers as well as the wake structure of anatomically realistic copepods (Borazjani et al., 2010). As in our previous works (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009; Borazjani and Sotiropoulos, 2010), we employed an efficient fractional step method to advance the flow equations in time (Ge and Sotiropoulos, 2007). The Poisson equation was solved with FGMRES (Saad, 2003), with multigrid as a preconditioner using parallel libraries of PETSc (Balay et al., 2011).

We obtained the kinematics of the C-start from experimental observations that captured the C-start motion of bluegill sunfish, Lepomis macrochirus, using high-speed cameras [see Tytell and Lauder (Tytell and Lauder, 2008) for details of kinematic and experimental hydrodynamic methods]. High-speed videos of bluegill C-starts and velocity vector fields of C-start flow patterns from that paper were used here to generate the input for computations and for comparison with the flow field output respectively. The kinematics data extracted from high-speed videos of bluegill C-starts consist of the instantaneous coordinates of 20 evenly spaced points (x_i and z_i) along the midline of the fish for 97 time instants, each 1 ms apart (Δt_exp=1 ms), during the entire turning maneuver [x_i=x_i(t_n) and z_i=z_i(t_n), where i=1,20 and t_n=nΔt_exp for n=1 to 97, as shown by blue circles in Fig. 1]. Because the spatial and temporal resolution of our simulations are much finer than that of the experiments (see below for more details), interpolation is required to reconstruct kinematics suitable to be used as input to the numerical model. We interpolated the experimental data in both time and space using cubic spline interpolation (Jacobs and Lott, 1989) to reconstruct the midline of the fish at each time step of the simulation as shown in Fig. 1. The cubic spline interpolation is implemented in our computer code based on the CMath library (Jacobs and Lott, 1989).

The procedure we employed to reconstruct the kinematics at any time instant t is summarized as follows. In general, and because the time step used in the numerical simulations is much smaller than Δt_exp, the time instant t was located within two experimental time instants where data were available (t_n<t<t_n+1). Therefore, we first preserved the spatial resolution of the experiments (20 points along the fish midline) and interpolated in time using cubic splines to obtain the coordinates of these 20 points at time t [x_i(t) and z_i(t)] shown by red circles in Fig. 1. Subsequently, we employed these 20 data points (red circles) to interpolate in space, also using cubic splines, to reconstruct the instantaneous shape of the fish midline at the desired numerical resolution. It should be noted that the midline obtained by the above procedure is smooth and does not have discontinuities at the data points, as shown by the black line that represents the reconstructed fine spatial and temporal resolution midline at time t passing through the red circles in Fig. 1.

The fish body geometry was re-constructed from a computed tomography scan of a bluegill sunfish similar to that used in the experiments of Tytell and Lauder (Tytell and Lauder, 2008). To facilitate the validation of our simulations with the experiments, we used the same bluegill sunfish geometry and reconstructed the kinematics used as input for the simulations from the specific experiment from which the PIV data that we use for comparison were obtained.

The forces on the fish body were calculated using two approaches: (1) the full 3-D approach, which requires knowledge of the instantaneous 3-D velocity and pressure fields all around the fish and can only be employed in the numerical simulations; and (2) an approximate 2-D approach, which uses instantaneous 2-D velocity data on a plane and can thus be employed both in the PIV experiments and the numerical simulations.

The accuracy of our numerical procedure for calculating the pressure and viscous forces for moving boundary problems has been demonstrated in Borazjani and Sotiropoulos (Borazjani and Sotiropoulos, 2008). Borazjani and Sotiropoulos (Borazjani and Sotiropoulos, 2008) simulated the flow induced by an axially vibrating cylinder and compared the results of their simulations with benchmark computational data (Dutsch et al., 1998). Excellent agreement was reported both for the total force and its two components.

To analyze the 3-D flow field, we employed fluid mechanics quantities that are typically used to quantify the kinematics of a fluid element, which in general can be expressed as the superposition of a rigid body translation, a rigid body rotation and a deformation. These three components of motion are quantified by the fluid velocity, vorticity and deformation (strain rate) fields, respectively (Panton, 1996).

Fig. 4 compares the PIV measurements (left column), inviscid simulation (middle) and viscous simulation (right) results at several time instants during the C-start on the midplane of the sunfish. The measurements and simulations in this figure are compared in terms of out-of-plane vorticity contours and in-plane velocity vectors. In what follows, we first briefly discuss the flow features observed in the experiments (for details, see Tytell and Lauder, 2008) and then compare the experimental results with the numerical simulations.

As the fish in the experiment (see left column in Fig. 4) begins to bend its body at the start of the turning maneuver in the middle of Stage 1 (t=21 ms), the tail arches to the right (positive x direction), generating a jet flow (Jet 1) and a clearly defined core of positive vorticity at its tip. The formation of Jet 1 is accompanied by the formation of Jet 2, which, as shown in Fig. 4B, transports fluid in the opposite direction and towards the fish body as required for mass conservation to be satisfied.

At t=41 ms in Stage 2 (Fig. 4C), the fish body has turned 90 deg with respect to its original orientation. The tail is now moving in the negative x direction to align with the rest of the body and as it does so it sheds a clockwise-rotating eddy (negative vorticity) that acts to strengthen both Jets 1 and 2. As seen in Fig. 4C, Jet 1 is now driven by the common flow between the two tail-shed vortices and has become wider and intensified as evidenced by the higher horizontal velocity component in this region. Jet 2 has also intensified and the increased momentum in the negative x direction provides the thrust needed to subsequently propel the fish forward along the positive x direction. This is clearly seen in Fig. 4D, where Jet 2 has fully formed into a strong propulsive jet along the x direction. A new feature of the flow that has emerged at t=41 ms is the formation of a third jet, Jet 3, which is located at the right side of the fish body (looking forward along the body) and is directed towards the body along the normal to the body direction. The formation of this jet is driven by the combined motion of the head and the body and is weaker overall than the other two jets.

During the final stage of the C-bend (see Fig. 4D), the tail has completed its counterclockwise rotation, moving out of the way of Jet 2 and allowing it to develop into a strong propulsive jet along the x direction. As a result, a strong shear layer develops at the interface between Jets 1 and 2, which, as seen in Fig. 4D, becomes unstable and breaks down into small-scale vortical structures. Jet 3 has moved towards the posterior of the body, where it is strengthened by counterclockwise (positive) vorticity entrained from Jet 2 as its lower shear-layer rolls up around the tail.

Both the inviscid and viscous simulations yield flow fields that essentially capture all of the major features of the flow observed in the experiments. In particular, the sequence of initiation of the various jets relative to the motion of the body, the spatial structure of the velocity field in the vicinity of the body, and the ultimate emergence of the characteristic triple-jet structure at the completion of the C-bend are reproduced with good accuracy by the numerical simulations. By definition, the inviscid simulations do not yield a viscous layer along the body because of the imposed slip-wall boundary condition, but they do capture the vortices generated around the body because of the acceleration imparted on the flow by the rapid body motion. The viscous simulation, in contrast, yields clearly defined viscous boundary layers around the body. The simulated boundary layer, however, appears to be much too thick compared with the measurements, which show little near-body generated vorticity. For the most part, this discrepancy should be attributed to the significantly lower Reynolds number in the simulations (4000 vs 275,000 in the experiment), but is also partly due to insufficient near-body resolution in the PIV experiments. The lack of near-body resolution in the experiments is also evident by the large velocities near the anterior of the fish body observed in the simulations. Such high velocities are not observed in the experiments, presumably because there are no measurement points sufficiently close to the body to resolve them.

Both the viscous and inviscid simulations capture the shear layer that develops during the end of the C-bend between Jets 1 and 2 (see Fig. 4D). In the viscous simulation, the shear layer breaks down into two smaller vortices, but in both the experiment and the invsicid simulations, the shear layer instability results in many small-scale vortical structures. This is to be expected because the shear-layer instability is a Reynolds-number-dependent phenomenon, which is suppressed by the relatively low Reynolds number used in the viscous simulations.

One persistent discrepancy between the experiments and the simulations is the angle of Jet 3 relative to the body. In the experiments, this jet is oriented nearly perpendicular to the body whereas in the simulations it is inclined at and angle of approximately 45 deg. The fact that both the inviscid and viscous simulations yield essentially the same orientation for Jet 3 suggests that this discrepancy is not a Reynolds number effect but should most likely be attributed to the fish not being completely in the horizontal xz midplane, i.e. the fish movement up and down or tilting relative to the horizontal midplane, which has not been accounted for in the reconstructed kinematics. Tytell and Lauder attempted to analyze only sequences in which the movements were primarily in the horizontal plane, but the behavior almost always included some rolling or vertical motion (Tytell and Lauder, 2008). Because the flow during C-start is highly three-dimensional (as will be shown in the next section), the flow in planes above or below the midplane can have higher out-of-plane but lower in-plane velocity. Also, the motion of the fish along the upward, downward or tilting directions changes the direction that momentum is transferred to the fluid by the fish body. All these factors can collectively change the velocity field near the body of the fish.

Finally, another discrepancy between the experiments and the simulations is with regard to the rate of decay of the positive vorticity eddy (VT1) shed by the tail at the very start of the C-bend (see Fig. 4B). At t=41 ms (Fig. 4C), VT1 is present in both simulations and its strength, as indicated by the magnitude of positive vorticity, is in good agreement with that calculated during the experiments. At the end of Stage 2 (Fig. 4D), however, it breaks down into two eddies with diminished strength in the simulated velocity fields, whereas it stays as a single strong vortex in the measurements. The reasons for this discrepancy are not entirely clear, especially because grid refinement tests have shown that the grid resolution is not an issue. As the initial acceleration of the tail will determine the strength of this eddy, we could, as suggested before, reasonably speculate that the fish was not completely horizontal during the C-start and might have moved up and down or tilted relative to the horizontal midplane. In our simulations, however, the fish has been assumed to stay on the same horizontal plane throughout the C-start maneuver. Also, the position of the tip of the tail is difficult to extract precisely from the videos, and some digitization error may contribute to the discrepancy.

The 2-D flow fields shown in Fig. 4 primarily illustrate the role of the tail in generating vorticity, but cannot show the vortices created by the dorsal and anal fins. To visualize these vortices, we plotted the velocity vectors and out-of-plane vorticity contours on a horizontal plane passing through the dorsal fin (Fig. 5). Note that for clarity we have zoomed in the area near the dorsal fin. It can be observed that at the beginning of Stage 1 a counterclockwise vortex (VF1) is created by the motion of the fin in the positive x direction (Fig. 5A). However, during Stage 1, as the fish body bends into a C shape, the fin’s motion changes direction and starts moving in the negative x direction. Consequently, a new vortex (VF2 in Fig. 5B) is created, rotating in the opposite direction of VF1. The common flow of these two vortices creates a jet flow in the positive x direction, similar to Jet 1 (Fig. 5B). As the fins continue to move along the negative x direction while the fish bends out of the C shape in Stage 2, new clockwise-rotating eddies are created (denoted as VF3 and VF4 in Fig. 5C). At the end of Stage 2, the vortices shed from the fins are breaking down into smaller vortices and the motion of the fin in the positive z direction now creates counterclockwise vortices (Fig. 5D).

The 2-D vorticity field observed in Figs 4 and 5 and discussed in the previous section comprises the footprints of an intricate 3-D wake, which could not be elucidated by the experiments. In this section we analyze the numerical simulations to present a complete description of the 3-D structure of the wake and establish linkages between the various jets identified in the previous section and coherent vortical structures. Fig. 6 shows several snapshots of the 3-D wake visualized using the iso-surfaces of the q criteria (Hunt et al., 1988) at the same time instants as those shown in Figs 4 and 5. The wake structure during the complete C-start maneuver is shown in supplementary material Movies 1 and 2. We begin by discussing the results of the inviscid simulations (Fig. 6, left column), followed by comparison and discussion of the inviscid and viscous simulation results.

At the beginning of Stage 1 (Fig. 6A), and as the body moves in the positive x direction, a tube-like vortical structure develops at the posterior of the body that outlines the fins and the tail (see Fig. 6A). This vortical structure consists of vortex tubes created by the fins (denoted as F1 in Fig. 6A) and the tail (denoted as T1 in Fig. 6A), which at this instant are connected together to form a single continuous starting vortex loop. Juxtaposing now Fig. 6A and Fig. 4A, it is evident that VT1 in Fig. 4A marks the intersection of the tail-attached vortical structure seen in Fig. 5A with the fish midplane. Similarly, juxtaposing Fig. 6A and Fig. 5A, we note that vortex VF1 is the footprint of vortex tube F1 on the plane passing through the dorsal fin plane.

As the tail continues its clockwise rotation and starts bending to form the C shape (Fig. 6B), the T1 and F1 vortex tubes are seen to deform along with the moving body. Each of the two F1 tubes detaches from its respective fin and assumes a step-like shape, connecting with the T1 tube, which still remains wrapped around the tail, along the right side of the tail via a flattened vortex sheet. As seen in Fig. 6B, yet another pair of vortex tubes (denoted by F2 in Fig. 6B) starts to develop at the edges of the dorsal and anal fins because of the motion of the fish body and the fins in the negative x direction. The sense of rotation in vortex loop F2 is opposite to that of the starting F1 loop, as demonstrated by their respective footprints, VF2 and VF1, on the dorsal fin plane shown in Fig. 5B. The sense of rotation of the F2 structures creates a flow in the negative x direction, which gives rise to the previously discussed Jet 2 in Fig. 4B.

At Stage 2, when the tail starts to move in the opposite direction to get out of the C shape, it creates a new tail-generated vortical loop (denoted as T2 in Fig. 6C), whose footprint is VT3 in Fig. 4C. This new tail-generated vortex loop connects to the vortical structure T1 (created by the previous motion of the tail and evolved into an approximately square-shaped vortex ring) via the two S-shaped tubes F1, which were generated earlier by the fins.

It can also be observed in Fig. 4C that a blob of negative vorticity, denoted as VT2, is located next to VT3. From this 2-D plot it appears as though VT3 has broken off from VT2 (Fig. 4C), which in 3-D terms corresponds to the pinching off of vortex loop T1 from vortex loop T2 to form the previously discussed square vortex ring. The complex transition from one vortex topology to another is illustrated in the schematic shown in Fig. 7 and involves pinch-off of the initial tail vortex loop T1 to form what we now refer to in Fig. 6C as vortex ring T1 and vortex loop T2.

Another new feature of the flow that is evident in Fig. 6C is the growth of vortex tube F2 and the generation of new similar vortex tubes from each fin, denoted as F3 and F4. All these fin-generated tubes start from the edge of each fin and connect with vortical structures generated along the tail. The footprints of F3 and F4 in the 2-D plane shown in Fig. 5C are eddies VF3 and VF4, respectively. The sense of rotation of these structures is similar to that of F2 and, as a result, they induce flow that acts to strengthen Jet 2.

The sequence of images shown in Fig. 6A–C shows that both the tail and the fins contribute significantly to the structure of the wake during the first two stages of the C-start and that vortex pinch-off and reconnection phenomena play an important role in the formation of the leading vortex ring T1 shown in Fig. 6C. As the tail continues its motion in the negative x direction and then in the negative z direction to get out of the C-bend at the end of Stage 2 (Fig. 6D), a similar sequence of events, involving formation of tail and fin generated vortex loops followed by vortex pinch-off and reconnection, gives rise to the intricate web of vortical structures observed in Fig. 6D, which consist of a series of vortex rings interlinked with vortex loops.

It is worth noting that both the inviscid and viscous calculations have very similar vortex topologies at all stages of the process. The only difference between the two flow fields is that for the latter case the vortical structures are smoother and somewhat more diffused and therefore easier to discern. Other than this anticipated Reynolds number effect, however, both simulations yield essentially the same 3-D wake dynamics.

To further highlight the complexity and three-dimensionality of the wake, we visualized the flow at the end of Stage 2 from different views and with different values of the quantity q (Fig. 8). Note that based on the definition of q (see Eqn 8), higher values of q indicate more intense vortical structures. Fig. 8A shows a side view of the fish with isosurfaces at q=1, whereas Fig. 8B shows the same view with q=20, highlighting only the strong vortices. In particular, the vortex loop T1 is visible in Fig. 8A at a low q value, but is not seen in Fig. 8B at the high q value, which indicates that its strength is lower relative to the strength of loops T2 and T3, which are present in both panels. The lower strength of T1 is in agreement with the breakup and diminished strength of its 2-D footprint, VT1 (Fig. 4D). It can also be observed that two additional vortex loops (denoted by F5 and F6 in Fig. 8B) are created by the anal and dorsal fins at the end of Stage 2. Similar to the F1 and F2 loops shown in Fig. 6B, the sense of rotation of these loops is opposite of each other as they are created by the motion of the fins in opposite directions, as demonstrated by their 2-D footprints (Fig. 5D).

We plotted the calculated time series of the x- and z-components of the force exerted by the flow on the sunfish during the C-start obtained from the inviscid and viscous simulations (Fig. 9). The force in this figure has been calculated using the 3-D approach, i.e. by integrating directly on the fish body the calculated instantaneous pressure and shear stress fields (see Eqn 2). Both the inviscid and viscous simulations exhibited the same overall trends, with the former yielding a somewhat more spiky distribution and larger magnitude extremes than the latter. This trend is consistent with the previously presented visualization of the respective flow fields, in which the inviscid simulations were characterized by a more complex overall wake with smaller-scale flow structures than their viscous counterpart. These differences notwithstanding, it is evident from Fig. 9 that both simulations yield force components with very similar characteristics. An important common feature is that, regardless of whether viscous effects are taken into account, the variation of the calculated forces is not monotonic and smooth, but instead is characterized by large-amplitude fluctuations occurring over a time scale of the order of a few milliseconds. During Stage 1, these force fluctuations are larger in the x-component of the force whereas during Stage 2 they are far more pronounced in the z-component of the force. This trend is obviously related to the rapidly changing orientation of the fish body, which is mainly oriented in the z-direction during Stage 1 and the x-direction during Stage 2, and the fact that the primary force component during each stage is oriented along the direction normal to the body.

We also calculated the force of each jet from the numerical simulations by applying the same method used to estimate the force from the PIV measurements of Tytell and Lauder (Tytell and Lauder, 2008). This is possible because the 2-D flow field in the midplane of the fish, where PIV measurements were obtained, can be readily extracted from our 3-D simulations. The computed momenta of the three jets obtained from the simulated instantaneous flow fields at the fish midplane were compared against the momenta calculated from the PIV measurements in Fig. 10A. It can be observed that the jet momenta from the numerical simulations are in excellent overall agreement with those obtained from the PIV measurements, especially when one takes into consideration the considerable uncertainty inherent in comparing numerically simulated flow fields with PIV data obtained from one realization of the C-start maneuver of a live fish. This level of agreement is of course not surprising considering that we have already established, through the comparisons shown in Fig. 4, that the simulated 2-D flow fields capture essentially all of the flow features observed in the PIV measurements. This level of agreement validates the numerically calculated force records obtained from the fully 3-D approach shown in Fig. 9 and points to the conclusion that the results of 2-D measurements should be treated with caution when used to derive even qualitative conclusions about the forces exerted on the fish body during complex and highly dynamic maneuvers such as the C-start.

One difference between the computed and measured results is the large-amplitude force fluctuations, such as those observed in our simulations (Fig. 10B), which were not observed to as large a degree in the experimental measurements of Tytell and Lauder (Tytell and Lauder, 2008). Recall, however, that the forces in the experiment were not measured directly; rather, they were calculated from PIV measurements using a 2-D control-volume approach (Tytell and Lauder, 2008). Smoothing due to the finite spatial and temporal resolution may limit the ability to detect these force fluctuations experimentally. Nevertheless, experimental force estimates do show fluctuations with a similarly short time course as the computed forces, though of smaller amplitude (compare dashed and solid blue lines in Fig. 10B).

The results shown in Figs 9 and 10 collectively suggest that the high-amplitude and frequency fluctuations of the force computed via the 3-D approach should be the result of a real physical process rather than a numerical artifact. To investigate the origin of these fluctuations, we plotted the calculated pressure field and fish midline velocity vectors at few instants in time (A–D) corresponding to selected pairs of peaks and troughs in the calculated forces shown in Fig. 9 (Fig. 11). The calculated pressure fields shown in Fig. 11 show that the large fluctuations of the force during Stages 1 and 2 are mainly related to the temporal fluctuations of the pressure field. For instance, it can be observed in Fig. 9, that at time instant A, the x-component of the force has a large positive peak whereas time instant B corresponds to a large negative peak. The visualized pressure field at time instant A reveals the presence of a large pocket of high pressure on the right side of the fish body (looking towards the head along the midline) and an equally large pocket of negative pressure on the left side of the body (Fig. 11A). At time instant B, however, the relative location of the pockets of high and low pressure has flipped, with the former now located on the left side of the fish near the head and the latter located on the right side of the fish near the tail (Fig. 11B). These large temporal fluctuations in pressure are entirely consistent with the large fluctuations in the x-component of the force between time instants A and B as shown in Fig. 9, which, in response to the pressure field, should fluctuate from a positive x-component force peak in Fig. 11A to a negative force peak in Fig. 11B. Similar trends are also observed between time instants C and D (Fig. 11C,D) and time instants E and F (Fig. 11E,F), thus pointing to the conclusion that the observed fluctuations of the computed force are entirely physical and are induced by the large temporal variation of the pressure field.

We have carried out the first 3-D simulation of the C-start of a bluegill sunfish with the kinematics prescribed from the experimental observations. The simulated instantaneous flow fields are compared against the PIV measurements in the midplane of the same fish from which we prescribed the kinematics. The level of agreement between the simulated and measured flow field is excellent, and is quite remarkable considering the uncertainties involved in simulating biological behaviors such as the C-start. One such uncertainty is the assumption that the swimmer stays in a horizontal plane and does not move upward or downward during the C-start. We estimate from our measurements that the sunfish has moved approximately 4 mm (0.04L) upwards during the C-start, but this displacement was not taken into account in our simulations. Another source of uncertainty can be attributed to the discrepancy in the Reynolds number of the simulations and the experiments, which is discussed in more detail in the next section.

The above limitations notwithstanding, our simulations have captured all of the flow features observed in the PIV measurements. The simulations have captured the three distinct jets observed in the experiments as well as the vortices shed and the breakdown of vortices at the end of Stage 2. The simulated flow fields exhibit general similarities with the slow turn of a giant danio studied experimentally by Epps and Techet (Epps and Techet, 2007). The vortices observed at the end of the giant danio turn [see fig. 5 of Epps and Techet (Epps and Techet, 2007)] are very much consistent with the vortices we found for the sunfish (Fig. 4D). Nevertheless, Epps and Techet (Epps and Techet, 2007) assumed that two vortex rings comprised the 3-D structure (denoted by Γ₁ and Γ₂ in fig. 5 of their paper) and ignored the middle vortex (VT3 in Fig. 4D), which we now know is the footprint of multiple connected vortex loops. Müller et al. (Müller et al., 2008) studied the rapid maneuver of a larval fish and observed similar vortices to those we observed in sunfish during both Stage 1 and Stage 2 of the maneuver.

Our numerical simulations provide a complete description of the 3-D wake structure responsible for creating the various jets and vortices observed in the PIV experiments. We found that a vortex loop (T1) is generated by the tail at the start of Stage 1 when the fish bends its body into a C shape. The flow at the center of this loop is the Jet 1 observed in the experiments. Through pinch-off, this loop is transformed at the start of Stage 2, when the tail starts moving in the opposite direction, into a square-shaped vortex ring interlinked with vortical structures generated by the fins and the tail. The footprints of this vortex ring are two patches of vorticity (vortices 1 and 2 in Fig. 4), which are clearly visible in the PIV measurements and simulations on the fish midplane. At Stage 1, another loop is created by the body and anal and/or dorsal fins (F2), and the induced velocity of this loop gives rise to Jet 2. At Stage 2, the movement of the tail to get out of the C-bend and the following tail beat adds to the complexity of the wake by creating a series of intertwined vortex loops and rings that induce flow that add momentum to Jet 2. The 2-D footprints of these connected loops have the same sign vorticity, and become unstable and break into smaller-scale vortices of the same sign.

The forces based on both the 3-D flow field (Fig. 9) and the 2-D midbody plane (Fig. 10) exhibited large-amplitude fluctuations over time scales spanning a few milliseconds. We argue that the experimentally documented rapid changes in the acceleration of the fish midline, which give rise to the large pressure fluctuations on the fish body, are physical and not the result of under-resolved experimental measurements. Although experimental measurements of insufficient temporal resolution could indeed yield spurious discontinuities in the second-order temporal derivative of the measured fish midline position, the fish kinematics used in the simulation were reconstructed from digitized experimental images taken at 1 ms intervals using cubic spline interpolation (see Fig. 1). Note that the complete C-start takes place within 60 ms, during which the orientation (position) changes completely. Therefore, it is reasonable to anticipate that the second derivative of position (acceleration) changes rapidly within a few milliseconds and that the 1 ms time interval of the experimental data is adequate to resolve these changes. Our results, therefore, make a strong case that the calculated large fluctuations in the force might indeed be necessary to change the position and orientation of the fish in such a short time.

Similar fluctuations were also observed in the simulations of Wolfgang et al. (Wolfgang et al., 1999) (see fig. 14 of their paper). However, the force fluctuations in their simulations (Wolfgang et al., 1999) were not as large and as fast as the fluctuations we observed here, probably because they simulated a slow turn that took approximately 700 ms whereas our turn only lasted approximately 60 ms. Two-dimensional experimental estimates of force also had rapid fluctuations, but not nearly at the same amplitude as the computed forces. The reduction in amplitude may be due to the finite temporal and spatial (2-D vs 3-D) resolution of the experimental measurements. As direct measurement of the forces in experiments with live fish is challenging, if not impossible, our work points to the necessity for the approach we adopted herein, i.e. coupling high-resolution numerical simulations with experiments, as the best path forward for understanding the dynamics of aquatic swimming, especially in complex situations such as the C-start maneuver.

We start by comparing the viscous and inviscid simulations in Fig. 4 (middle and right columns). A clear difference is the presence of a thick boundary layer next to the fish body in the viscous simulations. This is a viscous effect because as the Reynolds number is increased, the boundary layer becomes thinner and approaches zero in the inviscid limit. Another important difference between the two simulations is the fact that during the end of the C-bend (Fig. 4D) the inviscid simulations capture the experimentally observed breakdown of the shear layer between Jets 1 and 2, yielding a more complex and smaller-scale vorticity field compared with that obtained by the viscous simulations. This breakdown is a Reynolds number effect, because the shear layer becomes more susceptible to the onset of instabilities as the Reynolds number is increased. The Reynolds number used in the viscous simulation (Re=4000) is sufficiently low for viscous forces to stabilize the flow and suppress the experimentally observed shear-layer breakdown.

We observed similar patterns in the visualized 3-D structures in Fig. 6. Comparing the inviscid and viscous simulations, it is evident that the main 3-D features and vortex loops are present in both simulations. However, the viscous and inviscid simulations are different in the details of the 3-D wake in two ways. First, in the inviscid simulations, the vortex loops break up into smaller loops, whereas in the viscous simulations, the vortex structures are much smoother as the higher viscosity stabilizes the flow and prohibits the breakdown of vortices. Second, some 3-D structures are created near the head and around the body of the sunfish in the viscous simulation because of the thick boundary layer; these 3-D structures do not exist in the inviscid simulations.

These differences between viscous and inviscid simulations notwithstanding, it is evident that both simulations capture the key features of the flow during the entire C-start maneuver. This leads to the important conclusion that the vorticity dynamics of the flow during the C-start is driven primarily by the rapid acceleration of the fish body. Provided that the Reynolds number is sufficiently large, as the Re=4000 value used in the present simulations appears to be, the importance of viscous effects is secondary and restricted mainly to the near-body boundary layer and the late breakdown of the shear layers. Note that in the simulations we have lowered the Reynolds number by increasing the viscosity of the fluid. The forces were similar in both environments (inviscid vs viscous), mainly because the forces caused by the rapid acceleration of the body (added mass) dominate the viscous forces during the C-start. The Reynolds number of approximately 4000 is still in the transitional regime where the inertial forces dominate the viscous forces. The effectiveness of such kinematics in the viscous regime (e.g. for fish larvae) has yet to be tested.

Escape responses to flee predators are crucial for fitness. Consequently, studies of fast starts have been pivotal in evolutionary and ecological studies of predatory–prey interactions (Bergstrom, 2002; Domenici, 2001; Ghalambor et al., 2004; Hale, 1996; Harper and Blake, 1988; Harper and Blake, 1991). In fact, fish with slower or less-efficient fast starts tend to get captured more frequently than fish with a higher escape performance (Walker et al., 2005). The performance of the C-start is directly related to the hydrodynamics and how power is transferred from muscle to the fluid. There is still ongoing debate regarding how much of the power and/or force is generated in each stage of the C-start (Tytell and Lauder, 2008; Wakeling, 2006; Weihs, 1973). It has been argued, for instance, that Stage 1 is preparatory and does not contribute much to total thrust (Weihs, 1973). This assertion, however, has been countered by others who showed that thrust generation during Stage 1 is based on theoretical calculations (e.g. Wakeling, 2006). Such division of power and force generation between the two stages has consequences for both neural control and performance (Tytell and Lauder, 2008). More specifically, if substantial thrust is produced during Stage 1, then the Mauthner circuit that controls Stage 1 also has a direct effect on the overall escape performance (Tytell and Lauder, 2008). In this case, the whole body, which contributes to the C-bend, is crucial for force production. However, if Stage 2 is dominant, then the Mauthner circuit acts more like a trigger for a behavior in which other circuits may have a greater impact on performance (Tytell and Lauder, 2008). In this case, because Stage 2 involves more caudal fin movement than body movement (Domenici and Blake, 1997), the caudal fin would be more important for force output than the rest of the body (Tytell and Lauder, 2008).

From Fig. 9 it is evident that significant force is generated during Stage 1 of the C-start. Consequently, according to our results and in agreement with those of Tytell and Lauder (Tytell and Lauder, 2008), the Mauthner circuit has a crucial role in the performance of the C-start. Furthermore, we demonstrated that the forces based on inviscid and viscous simulations were quite similar (Fig. 9), which shows that the acceleration (added mass) forces are dominant during the C-start maneuver. The added mass forces are proportional to the acceleration of the body (Fig. 11), i.e. the whole body contributes to the force generation during the C-start. Based on this argument, we hypothesize that the tail and side fins contribute less than the main body to the force generation. In the future, we will test this hypothesis by comparing the forces generated by the C-start of a tailless and a finless virtual sunfish with those of the one tested and validated here against the experimental data. Such a virtual numerical experiment is similar to Webb’s experiments comparing the escape performance of an unmodified body shape in trout with the performance of fish on which he had amputated the dorsal and anal fins (Webb, 1977). Webb stated that his data were too variable to formulate conclusions on the effects of median fin amputation and he could not detect a significant effect of the dorsal and anal fins (supporting our hypothesis), despite the theory suggesting that the increase in depth of the fish (by fins) should increase the C-start performance (Weihs, 1973). Moreover, the anal and dorsal fins are controlled actively by the muscles (Jayne et al., 1996; Lauder et al., 2007; Standen and Lauder, 2005). During the C-start, these fins are rapidly erected to increase their surface area (Jayne et al., 1996; Tytell and Lauder, 2008). In this work we did not change the shape or motion of the fins during the C-start and the motion of the fins were assumed to be similar to the fish midline. The effect of active control of the fins during the C-start can be pursued as part of future investigations.

Another important feature of the C-start is the wake signature that it produces, which makes the fish conspicuous to predators. Based on our results (Figs 4, 5, 6), the 3-D wake structure is contained in a cube of size L×0.5L×L surrounding the fish in the x, y and z directions, respectively. This is a small volume for such a quick and high-velocity maneuver, considering the fact that the wake signature will not reach a predator at a distance approximately 2L from the fish until the fish is far enough away to avoid the predator. Additionally, predators are known to track fish using their wake flow patterns produced by steady swimming or low-frequency flapping (Dehnhardt et al., 2001; Hanke and Bleckmann, 2004; Wieskotten et al., 2010). Based on the data presented here, such tracking might prove difficult by a predator based on the C-start wake flows of the prey. At least three major jet flow directions are generated by a C-start with a complex 3-D pattern. If a predator in an environment where visual input is minimal can compute the direction of the escaping fish from the C-start wake flows, then that predator could in theory continue to pursue the prey. However, it remains to be determined whether a predator can track prey based on escape wake flow patterns, and this represents an intriguing area for future study. We can also compare the size of the wake signature of different species with different body forms and kinematics in future studies.

Finally, the effects of body shape on the hydrodynamics and the C-start performance of different fish is not completely understood. A previous study of Webb can guide our future numerical and experimental investigations (Webb, 1978). He compared seven different species of fish and found that among these species the bluegill sunfish had the best performance. He hypothesized that the body form of the bluegill that is dorso-ventrally flattened provides the best compromise for C-starts (Webb, 1978). Having demonstrated our capability to simulate numerically C-starts and validated our simulations against experiments, this enables us in the future to design and carry out virtual numerical experiments that can fully settle the above issues. Such numerical experiments comprise a powerful tool for biological research as they can be carried out under fully controlled conditions, are inherently free of extraneous effects and can easily disentangle the effects of different parameters (see Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009; Borazjani and Sotiropoulos, 2010).