Body dynamics and hydrodynamics of swimming fish larvae: a computational study

Many aquatic animals, especially fishes, rely on body undulations to swim. The kinematics and hydrodynamics of this so-called ‘body-and-caudal-fin’ (BCF) swimming has been studied extensively. Many experimental studies have focused on cyclic behaviours, which is also called ‘steady’ swimming, in spite of the unsteady nature of the fluid flow [anguilliform swimmers (Wardle et al., 1995; Müller et al., 2001; Tytell and Lauder, 2004); carangiform swimmers (Blickhan et al., 1992, Müller et al., 1997; Wolfgang et al., 1999; Nauen and Lauder, 2002; Lauder and Drucker, 2002; Müller et al., 2002)]. Although cyclic swimming may be common during migrations in open water and in early larval stages (Hunter, 1972), most fish swimming behaviours are not steady. Fish usually do not maintain a constant swimming speed and direction for many minutes; rather, they turn regularly and rely on burst-and-coast swimming (Blake, 1983; Müller and van Leeuwen, 2004).

The kinematics of cyclic and non-cyclic BCF swimming, in particular C-starts, have been studied extensively in adult fish (for reviews, see Videler, 1993; Domenici and Blake, 1997). There is also a growing body of literature on the hydrodynamics of swimming (for reviews, see Lauder and Tytell, 2006; Lauder, 2011). However, hydrodynamic studies have focused on cyclic swimming more so than non-cyclic behaviours (Wolfgang et al., 1999; Müller et al, 2008; Tytell and Lauder, 2008; Gazzola et al., 2012). Many of the experimental hydrodynamic studies use two-dimensional (2-D) particle image velocimetry (PIV) to describe the flow patterns (e.g. Blickhan et al., 1992; Müller et al., 1997; Wolfgang et al., 1999; Nauen and Lauder, 2002; Lauder and Drucker, 2002; Tytell and Lauder, 2004) and to reconstruct the three-dimensional (3-D) structure from multiple 2-D slices (Tytell and Lauder, 2004; Müller et al., 2008). However, volumetric PIV data are also becoming available (Flammang et al., 2011).

Experimental studies of cyclic BCF swimming show that fish generate two main wave patterns during cyclic swimming. Anguilliform swimmers have been shown to generate a double vortex wake (Müller et al., 2001; Tytell and Lauder, 2004). Carangiform swimmers generate a double vortex wake or a vortex chain wake, depending on the Strouhal number (Müller et al., 2002). During C-starts, fish shed two vortex rings: one during the preparatory phase and the other during the propulsive phase (Müller et al., 2000; Müller et al., 2008; Tytell and Lauder, 2008).

To understand the relationship between body shape and body movements on the one hand and the resulting flow patterns and swimming performance on the other hand, we need computational fluid dynamics (CFD) to provide data on propulsive forces (Liu et al., 1996; Liu et al., 1997; Liu and Kawachi, 1999; Carling et al., 1998; Kern and Koumoutsakos, 2006). Several studies varied tail beat frequency to compute wake topology over a range of Strouhal numbers from 0 to 1.3 for carangiform and anguilliform swimmers (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009; Borazjani and Sotiropoulos, 2010; Reid et al., 2012). These studies predicted for both anguilliform and carangiform swimming that a double-row vortex street forms at high Strouhal numbers and a single-row vortex street forms at low Strouhal numbers; some of these studies tethered the fish (Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009) whereas others allowed one degree of freedom – swimming speed was not an input parameter to the model but was computed from the hydrodynamic forces generated by the fish (Borazjani and Sotiropoulos, 2010). So far, two free-swimming models have been developed that allow three degrees of freedom: a 2-D model (Reid et al., 2012) and a 3-D model (Kern and Koumoutsakos, 2006). These models both found similar relationships between wake shape and body wave kinematics. CFD also provides more detailed 3-D flow fields than experimental flow visualization (Liu et al., 1996; Liu et al., 1997; Liu and Kawachi, 1999; Liu, 2005; Kern and Koumoutsakos, 2006; Borazjani and Sotiropoulos, 2008; Borazjani and Sotiropoulos, 2009; Borazjani and Sotiropoulos, 2010; Katumata et al., 2009).

So far, computational studies of swimming hydrodynamics have focused mainly on cyclic swimming and adult fish; far fewer studies have addressed unsteady swimming behaviours in larvae (Katumata et al., 2009; Gen et al., 2011; Gazolla et al., 2012). Fish larvae swim in an intermediate flow regime, with Reynolds numbers (Re) ranging from 10¹ to 10³ (McHenry and Lauder, 2005). Computational studies show that swimming at these lower Reynolds numbers requires higher Strouhal numbers compared with adult fish (e.g. Borazjani and Sotiropoulos, 2008). Experimental studies show that larval swimming differs in both kinematics and flow from adult swimming (for a review, see Müller, 2008). Hydrodynamic modelling of larval swimming puts special demands on computational models: (1) larvae swim with a wide body wave amplitude and therefore affect a relatively large volume of fluid that needs to be mapped with a fine grid to resolve near-body flows; (2) non-cyclic swimming kinematics cannot be modelled readily by a simple wave function; and (3) the fish's centre-of-mass (COM) movements are the result of the interaction between the forces exerted by the water and those exerted by the fish body, requiring the computational model to couple the computation of fluid and body forces.

Here, we present a 3-D computational approach to model a fish larva swimming freely in a horizontal plane. This model is based on detailed experimental data from a C-start and cyclic swimming of zebrafish larvae (Müller et al., 2008), which provide both input and validation data for the model. This study presents a computational model with three features: (1) a chimera grid system that reduces the computational costs of a large body wave amplitude; (2) a ‘play-back’ kinematics module that transfers experimentally observed body shapes into the model to study non-cyclic swimming movements; and (3) a numerical coupling of fluid and body dynamics that models free swimming.

This study had two goals. First, by comparing the output of this computational model with experimental flow and kinematic data obtained for exact same body kinematics, this study provided a rigorous comparison of CFD and experiment. We characterized the mechanics of BCF swimming for both cyclic and non-cyclic swimming (C-start). To this end, we computed the time- and space-resolved pressure and shear stress distributions on the body surface, the fluid-dynamic power and efficiency, as well as the time course of the 3-D structure of the generated wake. Second, we explored the effect of body wave amplitude on performance during cyclic swimming in order to discover optimal swimming strategies.

This study modelled swimming bouts of zebrafish larvae [Danio rerio (Hamilton 1822)]. We focused on two age groups, 3 and 5 days post-fertilization (d.p.f.), with body lengths of 3.8 and 4.4 mm, respectively (Fig. 1). To describe the swimming motion, we defined a right-handed Cartesian coordinate system with the x- and y-axes in the horizontal plane (Fig. 2A). To model the fish's body shape, we traced 50 points along the dorsal and lateral outlines of the finfold and the body of five zebrafish larvae per age group (Fig. 1). We then chose the larva that was closest to its age-group average and applied a surface mesh of 33×45 cells, scaled by body length. The computational model had no pectoral fins, consistent with our observations that pectoral fins were usually adducted during fast and cyclic swimming (Müller and van Leeuwen, 2004).

The surface grid was extended by 20 grid layers to a distance one-third of a body length away from the fish, resulting in a 33×45×20 body-fitted grid (Fig. 2B, pink mesh). Large body deformations during C-starts can cause the body-fitted grid to overlap with itself. Such geometric conflicts were avoided as necessary by locally reducing the depth of the body-fitted grid. In order to resolve the boundary layer, the thickness of the grid cells at the surface was defined as at most 0.1L√(Re); cells can become thinner in areas of large body deformations.

To describe the flow in the far field, we surrounded the body-fitted grid with a global Cartesian grid (Fig. 2C, cyan mesh) of sufficient volume to ensure that the flow generated by the fish is not influenced by the boundary conditions at the outer boundaries of the virtual tank. The body-fitted and global grids were combined using a multi-blocked, overset-grid chimera scheme (Liu, 2009; Prewitt et al., 2000), which generated a virtual cubic water tank with a virtual swimming object submerged in the tank. Global grid cells whose centres fell within the body grid were emptied. Global cells (Fig. 3A) bordering on empty cells were labelled as inter-grid boundary cells (IGBCs); cells on the outer boundary of the body grid were also labelled IGBC (Fig. 3B). Because the body-fitted and global grids were incommensurate, we obtained the value for a given IGBC from one grid by interpolation from the other.

Fig. 4 shows step-by-step the processes coupling flow to body dynamics, which ensured free swimming – the flow forces generated by the fish drive the fish's translational movement and its orientation in space. To couple flow and body dynamics, we derived flow solutions in a body-fitted grid and a global grid. At each time step, first the flow fields were calculated for each grid (including pressure and shear forces at the body surface), the solutions from each grid were then interpolated from the inter-grid boundaries (grid assembly), and then the force distribution on the body was used to compute the translational and rotational acceleration of the fish body.

We used non-dimensional quantities throughout the computation, from which we later derived a set of dimensional quantities in order to compare the model predictions with experimental results.

We used an in-house CFD solver based on the multi-blocked and overset finite volume method coupled to the Navier–Stokes (NS) equations (Liu, 2009). This solver computes unsteady biological fluid dynamics in the intermediate and low Reynolds number regimes.

The model had two sets of boundary conditions. Besides the aforementioned boundary conditions required to combine the body-fitted and global grids, solutions to the NS equations require further boundary conditions at the solid walls of the undulating body and at the external boundary of the global grid. On the body surface, we used the non-slip condition to define the velocity components. When computing pressure, we accounted for dynamic effects due to the body's acceleration. For the background global grid, the boundary conditions for the velocity and the pressure were: (1) upstream (u,v,w)=0 while pressure p was set to zero; (2) downstream and at the external boundaries of the global grid a zero-gradient condition was enforced for both velocity and pressure, i.e. ∂(u, v, w, p)/∂n=0, where n is the unit outward normal vector at the external boundary.

The computational swimming dynamic (CSD) model was constructed by modelling the zebrafish as a deforming body with a time-varying COM and moment of inertia. Although the model can solve the six-degree-of-freedom motion of a free swimmer, we locked the following three degrees of freedom: (1) the vertical component of the swimming speed was assumed to be zero because the zebrafish larvae executed routine swimming mainly in a horizontal plane, and (2) pitch and (3) roll were assumed to be zero based on observations during routine swimming.

By coupling the body dynamic model with the NS solver, we directly solved the nonlinear equations of body motion numerically to obtain translational and angular acceleration variables (a_COM,X, a_COM,Y and α) at each time step. Then we used the second-order Runge–Kutta integration method to obtain the three time-varying velocities (u_COM,X, u_COM,Y and ω_XZ). We neglected the internal force during the fish's deformation and assumed that the fish's rotation was identical to the heading angle. We conserved the fish's angular momentum at each time step.

In summary, the model fish can translate and rotate freely within the allowed three degrees of freedom. The body-fitted grid system was updated at each time step to take into account both the swimming dynamics and the prescribed fish body deformation. The model then rotated the fish to its updated heading angle, and translated it to its updated COM position.

To obtain the amplitude envelope of the fish's body wave, we superimposed the midlines from nine tail beats of one cyclic swimming bout (Fig. 6A) so that they collapsed into one point at the snout tip (Fig. 6B). We further assumed that the anterior 20% of the body was too stiff to undergo bending. We sampled the amplitude envelope at 11 points at 0.1L intervals along the x-axis to obtain an ‘average amplitude envelope’. This average envelope was symmetric and slightly narrower than the maximal envelope (Fig. 6B).

The analytical-function approach is used commonly in computational modelling and is useful for parameter space mapping. However, more complex time-dependent functions are required to model burst-and-coast swimming or starts, when the characteristics of the body wave change from one tail beat to the next. To model C-starts, we developed the play-back module, which used the observed body waves as input for the CFD model. This method not only readily reflected changes in body wavelength and speed, but also facilitated a direct comparison between simulations and experimental results.

The first aim of this study was to characterize the mechanics of BCF swimming. To simulate cyclic swimming in larval fish, we developed a wave function to describe the body movements of the computational larva based on the actual movements of a 3 d.p.f. zebrafish larva (Müller et al., 2008) (Fig. 6A). In this observed bout, the larva swam at a mean speed of 18±2L s⁻¹ along an approximately straight path (change in direction is less than 15 deg over nine tail beats). The larva generated a body wave travelling down its body with an almost constant tail-beat frequency (57±8 Hz, N=8) and amplitude (0.16±0.025L, N=8, measured as lateral excursion of the tail tip from the mean path of motion). Table 1 lists the kinematic parameters used in the simulation, which were based on the analysis of nine consecutive tail beats.

The forward swimming speed of the model was slightly lower (16L s⁻¹) than the observed swimming speed; it was at the lower boundary of the observed values of 18±2L s⁻¹ (Fig. 7). In the simulation, the fish started from rest (initial U_COM=0) and reached the final speed after nine tail beat cycles (=18 tail beats). The discrepancy between the model's final speed and the experimentally observed speed might be caused by differences in the time history leading up to the experimentally observed swimming sequences, by slight discrepancies in morphology between the model and the real fish, or most likely, by the fact that the real fish changed its body wave shape from tail beat to tail beat, leading to small changes in swimming speed and swimming direction that did not occur in the model fish, which swam with the same symmetric tail beat amplitude envelope during the entire swimming bout. The model fish's average propulsive efficiency during cyclic swimming was 30.2% (averaged over a tail beat cycle; Table 2).

The free-swimming model revealed that the shape of the amplitude envelope was the result of fluid–body interaction. Fluid–body interaction caused the fish to side-slip and to yaw (which are excluded in a tethered computational model, explained in the  Appendix), resulting in an amplitude envelope whose shape depended on the frame of reference: although the tail amplitude was wider in the fish than in the earth frame of reference (0.27L versus 0.16L), the snout amplitude was, by definition, wider (0L versus 0.04L), given that the origin of the fish frame of reference was the snout.

Overall, the COM kinematics (translation and rotation) of the computational model closely resembled that of the actual fish larva. We therefore concluded that the model was sufficiently accurate to imitate the interaction between free-swimming dynamics and hydrodynamics during cyclic swimming.

Besides the body dynamics, the CFD model also accurately predicted the flow topology and time course. The symmetric body waves of the computational fish caused small differences with the experimental flow fields. However, CFD and experimental flow fields shared salient flow characteristics.

To validate our CFD data, we compared spatial and temporal flow patterns of the CFD-generated flows with the experimentally observed PIV patterns. The first step of the validation was the identification of the correct horizontal transect through the 3-D CFD flow field (Fig. 8). The following evidence suggested that the PIV flow field transected the larva ventrally near transect CFD-C. First, the distribution and size of the areas of elevated vorticity in CFD-C were very similar to the PIV results: the vorticity area (red area, labelled a) near the head extended as a narrow band from the snout to the peak of the wave crest, and the shed vortex pair (Fig. 8, blue area, labelled b) was ~0.3L long and 0.1L wide. Second, the magnitude of the vorticity peaks was similar in CFD and PIV [Fig. 8; vorticity peak in the wake (blue area, labelled b): 50 s⁻¹; vorticity peak near the head (red area, labelled a): −55 s⁻¹]. In conclusion, the CFD model generated horizontal transects that resembled the PIV data and helped to identify the location of the light sheet in the PIV experiments: the horizontal transect through the 3-D flow field was approximately 0.06L (0.34 mm) below the medio-frontal plane of the larva, only the yolk sac and the ventral finfold were situated in the light sheet.

The second step of the validation was a comparison of the temporal structure of the flow fields. We compared the PIV transects with CFD transects at nearly identical dorso-ventral positions (Fig. 9). Again, PIV and CFD shared salient flow characteristics. Clearly visible near the anterior body was an area of elevated vorticity in the boundary layer enveloping the head and yolk sac. This boundary layer separated at the eyes and remained close to the body until it was disturbed by the highly unsteady flow generated by tail beat. Also visible along the anterior half of the body was the drag wake that formed at the widest point of the larva at the protruding eyes and yolk sac. This drag wake extended posteriorly from the head and yolk sac and co-determined the flow along the posterior body. Along the posterior body, the model also correctly predicted the undulation-based vorticity regions. These regions formed halfway along the body, strengthened as they moved down the tail, and were finally shed at the tail tip and formed a distorted vortex ring, which appeared as vortex pairs in the horizontal cross-section of the flow field (medio-frontal plane) (Fig. 9). Each lateral tail beat produced one vortex pair in the horizontal cross-section. In conclusion, CFD accurately predicted the topology and time course of the flow patterns over a tail beat cycle.

Pressure distribution and hydrodynamics of the flow near the larva CFD can quantify flow characteristics that are difficult or impossible to measure, such as pressure and shear distributions on the fish body, which are important stepping stones towards understanding the hydromechanics underlying cyclic swimming. Fig. 10 shows the time-dependent relationship of body deformation, jet flow and surface pressure within one tail beat cycle.

We modelled an episode of cyclic swimming in which the body wave length was 1L. As the body wave travelled down the body, the deepening wave troughs caused negative-pressure regions and the growing crests caused positive-pressure regions near the inflexion points of the body wave (Fig. 10, emphasized by blue and red dashed lines indicating corresponding regions in the dorsal and lateral view). Note that the negative-pressure peaks were located more distally (at the finfold), while the positive-pressure peaks were located more centrally (halfway along the dorso-ventral extent on the body) (Fig. 10, dorso-ventral distributions are only visible in the lateral views, not the dorsal view). Negative- and positive-pressure regions always appeared in pairs near an inflexion point of the body wave, where the angle between the body surface and the forward direction of motion of the COM was maximal (Fig. 10, sketch). This position at the inflexion points of the body wave ensured that the force on the body due to the pressure differentials across the body had a substantial thrust component. Overall, the posterior body generated the highest force peaks and supplied most of the thrust, consistent with findings of other CFD models that evaluated pressure and shear stress in adult fish (Borazjani and Sotiropoulos, 2008; Reid et al., 2012).

With a body wave length of 1L, there were two pairs of high-pressure regions on the fish body at any given time (Fig. 10, sketch). The resultant force of the two pairs caused a slight yaw around the COM near the yolk sac – note that the position of the COM fluctuated slightly in the fish frame of reference due to the body undulations.

The pressure differentials across the body caused jet flows, most prominently towards the troughs of the body wave. The jets first formed just posterior of the yolk sac, pointing in a medio-rostral direction (Fig. 10, flows 1 to 3 at 1/6T to 3/6T). As the trough jets travelled posteriorly, they increased in speed and their flow direction changed gradually to a medio-posterior orientation (Fig. 10, flows 4 and 5 at 4/6T to 5/6T). By the time the trough jets reached the tail and were shed, strong jets also became visible at the crests (Fig. 10, flows 6 and 7 at 6/6T and 1/6T). These backward jet flows were a testament to the pressure differentials across the body that pushed the larva forward.

To show the close relationship between jet flows and vortex ring structures, we selected 12 vertical (Fig. 11, panels YZ1 to YZ12) and five horizontal cross-sections (Fig. 11, panels XY0 to XY–2) that illuminate the flow structures surrounding the swimming larva. Two undulation-based jet flows formed at the two crests and troughs (Fig. 11, jet flows 1 and 2). Both in horizontal and vertical cross-sections, the jet flows were sandwiched by vorticity regions of opposite rotation. This vorticity pattern indicated that two incomplete vortex rings (Fig. 11, rings Ra and Rb) had already formed near the larval body, with a jet flow (Fig. 11, golden arrows JF1 and JF2) located in the centre of the ring structure. The two vortex rings touched each other near the inflection point of the body wave (Fig. 11, panel YZ9). This close proximity might cause the flows of the two vortex rings to interact as indicated by the blue arrows in the central panel of Fig. 10. In the vertical cross-sections, the YZ vorticity showed that the two vortex rings were clearly separate entities at the point of contact of the two rings (Fig. 11, panel YZ9). However, no such separation between the vortex rings was visible in the XY plane (Fig. 11, panel XY0).

As the crests and troughs reached the tail tip, the jet flows and their vortex rings were shed in the wake, each ring now forming a complete loop. Newly shed vortex rings were prominent and initially remained closely together (Fig. 12, see rings A and B at 1/5T and rings B and C at 3/5T), but then gradually drifted apart as they moved away from the path of motion under their own momentum (Fig. 12, separation of rings A and B at 1/5T–3/5T and rings B and C at 3/5T–5/5T). The process of separation led to the patches of XY vorticity being stretched before the vortex rings separated completely. This process manifested as vortex rings that were initially stretched towards the path of motion rather than forming perfect tori (Fig. 12, left column), before vorticity dissipated sufficiently in these stretched regions to cause incomplete (U-shaped) tori in the more mature wake (Fig. 13).

In this cyclic swimming case, the wake comprised a staggered array of vortex rings along the mean path of motion. Unlike the single chain of vortex rings reported for adult carangiform swimmers (e.g. Blickhan et al., 1992), wake structures in this case comprised double-row vortex ring structures, which agreed roughly with the extrapolation from experimental results by Müller et al. (Müller et al., 2008). Vortex rings in the same row were parallel to each other, but they were not parallel to the mean path of motion or the rings of the contra-lateral row. The vortex rings continued to drift away from the mean path of motion under their own momentum as vorticity diffused and flow speed declined over time due to viscous effects.

Note that P_R only takes into account hydrodynamic power requirements but not muscle or metabolic terms. We did not observe an optimal amplitude within the examined range of amplitudes that would maximize any of the considered performance criteria. Swimming speed, power and power requirements per unit distance (P_R) were proportional to body wave amplitude. Of the four examined optimization parameters, only efficiency showed a strong slowing trend: decreasing the amplitude by 40% caused a 0.09 drop in efficiency while increasing the amplitude by 40%caused only a 0.01 increase in efficiency (Fig. 14D). We conclude, first, that the hydrodynamic efficiency of the experimentally observed amplitude was close to the predicted maximum efficiency attained at the highest wave amplitude. Second, fish can increase their swimming speed by increasing their body wave amplitude; our data suggest that fish larvae might be limited in their swimming speed more by physiological (e.g. muscle strains) and morphological (e.g. mechanical limits to body curvature) rather than hydrodynamic constraints. Third, increasing swimming speed through increasing body wave amplitude increased power and power consumption per unit distance. For fish larvae, swimming faster means proportional increases in hydrodynamic power requirements and transport costs. So, an increase in efficiency came at the price of increased power requirements. Fish larvae appeared to use a body wave amplitude during cyclic swimming that struck a balance between desired speed, efficiency and available muscle power rather than an amplitude that purely maximized hydrodynamic efficiency.

In the swimming event shown in Fig. 15, a 5 d.p.f. zebrafish larva spontaneously initiated a swimming bout with a turn (spontaneous C-start). During the first tail beat (preparatory stroke), the larva bent its entire body into a C-shape. During the subsequent tail beat (propulsive stroke), the larva accelerated forward. Following the C-start, the larva continued to swim for another five tail beats, while tail beat amplitude and swimming speed continually declined until the larva kept its body straight after the seventh tail beat and began to coast. The initial phase of this swimming event (0 to 70 ms) was published in Müller et al. (Müller et al., 2008).

To compute a start, we drove the fish model using the play-back module, which derived the fish's body undulations directly from the experimental data. The computational trajectory of the larva's COM closely resembled the experimental trajectory during the entire swimming bout (Fig. 16). The preparatory stroke (tail beat 1: 0–23 ms) resulted in only a slight lateral translation of the COM (Fig. 16A, Fig. 15; tail beat 1 from 0 to 23 ms). During the propulsive stroke (Fig. 16A, Fig. 15; tail beat 2 from 27 to 50 ms), the COM was accelerated forward and the fish reached maximum forward speed.

The largest differences in COM kinematics between the computational prediction and the experimental observation occurred during the preparatory stroke of the C-start: the computation predicted a 41 deg angle compared with 44 deg in the experiment; angular speed peaked at 2.1 deg ms⁻¹ after 9 ms in the experiment, compared with 2.8 deg ms⁻¹ after 13 ms in the computation (Fig. 17C,D). CFD and experimental observation agreed well regarding magnitude and time to maximum speed (7.4 versus 7.0L s⁻¹ after 39 ms; Fig. 17A). The propulsive stroke of the C-start (tail beat 2) produced not only the highest speed and power but also the highest efficiency (Fig. 17B) compared with the other tail beats of this swimming bout.

The computational model successfully predicted the flow patterns and vorticity levels generated during a C-start (Fig. 18). The preparatory stroke (10 ms) generated a strong jet flow towards the trough of the body wave plus two weaker jet flows adjacent to the head and the tail. The pressure differentials across the body that caused these jets also generated four patches of elevated vorticity along the body (Fig. 18, sketch, vortices 1–4), which later travelled down the body and were shed at the tail during stroke reversal. The good agreement between vorticity distributions and magnitudes of CFD and PIV flow fields suggested that the PIV transect was situated near the midline of the fish.

The computational model also correctly predicted the timing of the events, such as the splitting and shedding of vortex rings. Fig. 18 shows a transect through the medio-frontal plane of the wake. In this cross-sectional view, vortex rings appear as vortex pairs. Vortex pair 1–2 was shed during tail beat 1 at 30 ms, and another vortex pair 2b–3 was shed during tail beat 2 at 50 ms; vortex pair 3–4 was stretched, but not completely shed. These vortex pairs in fact form cross-sections through 3-D vortex rings (see below).

During the preparatory stroke, two regions with an elevated pressure across the body formed along the larva's body, one at the tail and one mid-body (Fig. 19, 10 ms, colour contour panels Lp and Rp). The pressure region at the tail produced a force that caused the fish to turn anti-clockwise. It was shed at the tail during the propulsive stroke 25 ms after the initiation of the start, forming a jet that died off quickly (Fig. 19, jet 1). The pressure region at the mid-body generated an opposite jet flow (Fig. 19, jet 2) that was most prominent in the trough of the C-bend. During the propulsive stroke, this pressure region travelled towards the tail, and the jet flow into the trough gradually changed direction, from being initially perpendicular to the body until it pointed mostly posteriorly and ran parallel to the body by the time the jet reached the tail. As the jet travelled down the body, the flow speeds in the jet continued to increase, accelerating the fish forward. The pressure forces of the preparatory stroke mainly caused an angular acceleration of the larva that changed the larva's heading. The pressure forces acting during the propulsive stroke mainly caused a translational acceleration of the larva's COM, but they also counteracted the anti-clockwise rotation initiated during the preparatory stroke, ending the turning motion and putting the COM on an approximately straight path of motion.

The larva generated not only pressures but also shear stresses (Fig. 19, colour contour panels Ls and Rs). The shear forces were lower than the pressure forces and mainly acted as drag. The largest shear forces arose at the head during the preparatory stroke, counteracting the yaw of the head, while the shear forces at the tail generated by the local jet flow were much lower and mainly acted as thrust. High shear forces also occurred along the finfold.

Overall, the larva generated mostly rotational forces during the preparatory stroke, causing a change of heading. During the consecutive propulsive stroke, the larva generated (counteracting) rotational and translational forces, which ensured that the larva proceeded along the new heading without further changes in heading.

Fig. 20 shows the 3-D vortex pattern as iso-surfaces at Q=0.02. Just like during cyclic swimming, we observed vortex rings with central jets. The preparatory stroke shed the tail jet and its surrounding vortex ring (Fig. 20, vortex ring A; cf. Fig. 18, vortex pair 1–2); the propulsive stroke shed the mid-body jet and vortex ring. Consistent with the angular and translational accelerations occurring during the generation of these vortex rings, the jet of the preparatory-stroke ring pointed slightly sideways and forward (Fig. 20); the jet of the propulsive-stroke ring pointed backward (Fig. 20). Both vortex rings quickly moved apart under their own momentum (Fig. 18, split of vorticity 2 and 2b).

Previous findings on adult fish that the preparatory stroke already causes a translation (e.g. Weihs, 1973; Webb, 1976; Domenici and Blake, 1991) could not be confirmed. Instead, the flow fields generated by the larva support the observed body and COM kinematics – the wake showed clear backward momentum only in the propulsive but not the preparatory stroke. Our findings, however, are consistent with those of another CFD study on larval fish (Gazzola et al., 2012), which computed similar wake and fluid force patterns based on the same experimental swimming sequence. That study also found that increasing the amplitude envelope increases swimming speed, consistent with our findings for cyclic swimming.

G.L. would like to specially thank his Japanese friends who helped and encouraged him through the harsh time of the Great East Japan Earthquake in 2011.

As our study of the effects of amplitude on swimming performance has shown, the magnitude of the body wave amplitude has a profound effect on swimming performance. However, many CFD models effectively tether the model swimmer, allowing the fish to only swim forward, but not to slip sideways (lateral displacement) and yaw (rotation). When imposing a body wave on the undulatory swimming using an analytical description of the body wave akin to Eqn 14 but not allowing side-slip and rotation, tethering causes CFD models to overestimate the amplitude of the tail beat in an earth frame of reference (supplementary material Fig. S1B). In our zebrafish larva, the lateral excursion of the tail in the earth frame of reference would be 0.16L in the tethered fish instead of 0.09L in the free-swimming fish. Our three-degree-of-freedom free-swimming CFD model allows lateral forces and moment to act on the fish, and we observe both side slip and yaw, which reduces the effective lateral excursions of the tail in an earth frame of reference. As shown in supplementary material Fig. S1, allowing slip results in a large reduction of propulsion (supplementary material Fig. S1C): the cyclic swim CFD speed is only 6L s⁻¹, a fraction of the observed swimming speed (18L s⁻¹).

Models that tether the fish can obtain more realistic results by using an amplitude-envelope function obtained from midlines in an earth frame of reference, in which the midlines arrange to form a pivot point near the rear of the head and show a strong head yaw. When using a free-swimming CFD model, which does not ignore the lateral forces and yaw, such an approach would introduce a severe error because the effects of hydrodynamic forces on lateral position and yaw would be included in the model twice, first embedded in the amplitude envelope equation, and then again in the flow solution. In fact, using such an amplitude envelope with an earth frame of reference leads to a severe underestimation of the tail beat amplitude and, consequently, a severe underestimation of swimming speed (supplementary material Fig. S1). To circumvent this problem, our model used an amplitude-envelope function that was derived from the fish midlines in a fish frame-of-reference with its point of origin at the snout tip. In this way, the model achieved a swimming speed of 16 m s⁻¹, a value very close to the experimentally observed value of 18 m s⁻¹.

 
• a(a_COM,X, a_COM,Y)

  translational acceleration

 
• a(x_axis)

  body wave amplitude-envelope function at x_axis

 
• BCF

  body and caudal fin

 
• COM

  centre of mass

 
• D

  drag

 
• D*

  dimensionless coefficient of drag

 
• d.p.f.

  days post-fertilization

 
• f

  frequency

 
• F(F_X, F_Y,, F_Z)

  force vector

 
• G(G_X, G_Y,, G_Z)

  force vector

 
• H(H_X, H_Y,,· H_Z)

  force vector

 
• IGBC

  inter-grid boundary cell

 
• I_XZ

  moment of inertia about the y-axis

 
• JF

  jet flow

 
• L

  body length

 
• M

  body mass of fish

 
• m_i,j,k

  mass of a body element

 
• n

  unit normal vector

 
• NS

  Navier–Stokes equation (symbols used in NS not included in list)

 
• N_XZ

  hydrodynamic moment about the y-axis

 
• p

  pressure

 
• P*_total

  dimensionless coefficient of power

 
• P_body

  pressure acting on the body

 
• P_hydro

  hydrodynamic power

 
• P_iner

  inertial power

 
• PIV

  particle image velocimetry

 
• P_total

  total power

 
• r

  positional vector from the COM to the centre of a grid cell

 
• Re

  Reynolds number

 
• Re_exp

  Reynolds number used in computation

 
• S

  body surface area

 
• S_body

  shear stress acting on body

 
• S_i,j

  area of surface element

 
• t

  time

 
• T

  thrust

 
• T*

  dimensionless coefficient of thrust

 
• T*_net

  dimensionless coefficient of net thrust

 
• T_net

  net thrust

 
• u

  velocity at the surface of a grid cell

 
• u_COM(u_COM,X, u_COM,Y, u_COM,Z)

  translational velocity vector of the COM

 
• V(u, v, w)

  fluid velocity vector

 
• VR

  vortex ring

 
• x_axis

  length along the body axis

 
• X, Y, Z

  world coordinate system

 
• (X_COM, Y_COM, Z_COM)

  position vector of the COM

 
• z(x_axis, t)

  propagating body wave function

 
• α

  angular acceleration

 
• η_pe

  hydrodynamic propulsive efficiency

 
• λ

  wavelength

 
• ν

  water viscosity

 
• ρ

  water density

 
• ω_XZ

  angular velocity about the y-axis

 
• ξ, ζ, η

  computational grid coordinate system