SUMMARY
To understand the mechanics of fish swimming, we need to know the forces exerted by the fluid and how these forces affect the motion of the fish. To this end, we developed a 3D computational approach that integrates hydrodynamics and body dynamics. This study quantifies the flow around a swimming zebrafish (Danio rerio) larva. We used morphological and kinematics data from actual fish larvae aged 3 and 5 days post fertilization as input for a computational model that predicted freeswimming dynamics from prescribed changes in body shape. We simulated cyclic swimming and a spontaneous Cstart. A rigorous comparison with 2D particle image velocimetry and kinematics data revealed that the computational model accurately predicted the motion of the fish's centre of mass as well as the spatial and temporal characteristics of the flow. The distribution of pressure and shear forces along the body showed that thrust is mainly produced in the posterior half of the body. We also explored the effect of the body wave amplitude on swimming performance by considering wave amplitudes that were up to 40% larger or smaller than the experimentally observed value. Increasing the body wave amplitude increased forward swimming speed from 7 to 21 body lengths per second, which is consistent with experimental observations. The model also predicted a nonlinear increase in propulsive efficiency from 0.22 to 0.32 while the required mechanical power quadrupled. The efficiency increase was only minor for wave amplitudes above the experimental reference value, whereas the cost of transport rose significantly.
INTRODUCTION
Many aquatic animals, especially fishes, rely on body undulations to swim. The kinematics and hydrodynamics of this socalled ‘bodyandcaudalfin’ (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 burstandcoast swimming (Blake, 1983; Müller and van Leeuwen, 2004).
The kinematics of cyclic and noncyclic BCF swimming, in particular Cstarts, 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 noncyclic 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 twodimensional (2D) 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 threedimensional (3D) structure from multiple 2D 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 Cstarts, 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 doublerow vortex street forms at high Strouhal numbers and a singlerow 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 freeswimming models have been developed that allow three degrees of freedom: a 2D model (Reid et al., 2012) and a 3D model (Kern and Koumoutsakos, 2006). These models both found similar relationships between wake shape and body wave kinematics. CFD also provides more detailed 3D 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^{1} to 10^{3} (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 nearbody flows; (2) noncyclic swimming kinematics cannot be modelled readily by a simple wave function; and (3) the fish's centreofmass (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 3D computational approach to model a fish larva swimming freely in a horizontal plane. This model is based on detailed experimental data from a Cstart 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 ‘playback’ kinematics module that transfers experimentally observed body shapes into the model to study noncyclic 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 noncyclic swimming (Cstart). To this end, we computed the time and spaceresolved pressure and shear stress distributions on the body surface, the fluiddynamic power and efficiency, as well as the time course of the 3D 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.
MATERIALS AND METHODS
Geometric model and grid system
This study modelled swimming bouts of zebrafish larvae [Danio rerio (Hamilton 1822)]. We focused on two age groups, 3 and 5 days postfertilization (d.p.f.), with body lengths of 3.8 and 4.4 mm, respectively (Fig. 1). To describe the swimming motion, we defined a righthanded Cartesian coordinate system with the x and yaxes 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 agegroup 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 onethird of a body length away from the fish, resulting in a 33×45×20 bodyfitted grid (Fig. 2B, pink mesh). Large body deformations during Cstarts can cause the bodyfitted grid to overlap with itself. Such geometric conflicts were avoided as necessary by locally reducing the depth of the bodyfitted 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 bodyfitted 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 bodyfitted and global grids were combined using a multiblocked, oversetgrid 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 intergrid boundary cells (IGBCs); cells on the outer boundary of the body grid were also labelled IGBC (Fig. 3B). Because the bodyfitted and global grids were incommensurate, we obtained the value for a given IGBC from one grid by interpolation from the other.
Solution process to model free swimming
Fig. 4 shows stepbystep 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 bodyfitted 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 intergrid 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 nondimensional quantities throughout the computation, from which we later derived a set of dimensional quantities in order to compare the model predictions with experimental results.
Computational fluid dynamic model
Governing equations
We used an inhouse CFD solver based on the multiblocked 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 governing equations are the 3D, incompressible and unsteady NS equations written in a strong conservative form for mass and momentum: (1) where where u, v and w are the velocity components in the Cartesian coordinate system, p is pressure, t is physical time and Re is the Reynolds number. The equations also contain the pseudo time parameter τ and pseudocompressibility coefficient λ. The term q associated with pseudo time was designed for an inneriteration at each time step, and vanished when the divergence of velocity approached zero so as to satisfy the equation of continuity. By introducing the generalized Reynolds transport theorem and by employing the Gauss integration theorem, Eqn 1 was transformed in a general curvilinear coordinate system as: (2) where f=(F+F_{v}, G+G_{v}, H+H_{v}); S(t) denotes the surface of the control volume; n=(n_{X}, n_{Y}, n_{Z}) are the components of the unit outward normal vector corresponding to all the surfaces of a polyhedral cell; and u_{g} is the local velocity of the moving cell surface. For a structured, boundaryfitted and cellcentred storage architecture, we rewrote Eqn 2 into a semidiscrete form: (3) where for example, The parameters and refer to the two opposite surfaces i of a particular computational cell as indicated by the subscript and . The parameters , , and were derived in a similar manner and refer each to two opposing cell surfaces j and k. The term V_{ijk} is the volume of the cell (i, j, k). The superscript ξ denotes the ξdirection perpendicular to two particular cell surfaces (indexed by i, whereas the ζ and ηdirections are indexed by j and k, respectively; see Fig. 5). A more detailed description can be found in Liu and Kawachi (Liu and Kawachi, 1999) and Liu (Liu, 2009).
Fortified solution algorithm
The geometries of the body grid and the global grid were incommensurate. It was therefore necessary interpolate values in two overlap regions between the grids (IGBCs) and these interpolated values needed to be fortified to enhance stability of the algorithm. The solution algorithm, which is based on a fortified NS approach, was derived by adding forcing terms σ(q_{f}–q) to the semidiscrete Eqn 3: (4) where q_{f} is derived from the interpolation (flow solution of the other overlapping grid in the last time step) in the IGBCs, and set to zero in the rest of the computational domain. In the IGBCs, the switching parameter σ was given a large positive value such that the NS equation algorithm was effectively turned off and reduced simply to q=q_{f}, so the solution was fortified there. In the rest of the computational domain, we defined σ=0, which reduced the algorithm to the standard solution. More details and further derivations from Eqn 4 can be found in Liu and Kawachi (Liu and Kawachi, 1999) and Liu (Liu, 2009).
Boundary conditions
The model had two sets of boundary conditions. Besides the aforementioned boundary conditions required to combine the bodyfitted 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 nonslip 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 zerogradient 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.
Computational swimming dynamic model
The computational swimming dynamic (CSD) model was constructed by modelling the zebrafish as a deforming body with a timevarying COM and moment of inertia. Although the model can solve the sixdegreeoffreedom 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.
The deformablebody dynamic model was constructed based on the Newton–Euler equations of six degrees of freedom of body motion. By locking three degrees of freedom, we reduced these equations to a set of three coupled nonlinear ordinary differential equations: (5) where F is the fluid force vector acting on the body COM, u_{COM}=(u_{COM,X}, u_{COM,Y}) is the translational velocity of the COM of the body, a_{COM} is the translational acceleration of the COM, ω_{XZ} is its global angular velocity about the yaxis, α is its global angular acceleration about the yaxis, M is the body mass, I_{XZ} is the yawing moment of inertia (which depends on the instantaneous body shape) and N_{XZ} is the hydrodynamic moment about the COM. Body mass was assumed constant but the yawing moment of inertia varied with time and was updated at each time step. The earth frame of reference was used for all dynamic computations to avoid any additional terms that the body deformations might otherwise introduce into the equations of motion. After the NS equations were solved, the two force components (F_{x},F_{z}) exerted on the body surface in the global system were evaluated by a summation of the force vectors of all the cells at the body surface due to the local pressures and shear stresses: (6) where A_{i,j} is the surface area of the i,jth cell on the body surface, and P_{body,i,j} and S_{body,i,j} are pressure and shear stress vectors on this surface cell, respectively, which arose from the solution of NS equations.
The hydrodynamic moment N_{XZ} was calculated as the sum of the crossproduct of the positional vector and the force at each computational cell centre on the body surface, such that: (7) where r_{i,j} is the positional vector from the COM to the i,jth cell centre on the body surface and F_{i,j} is the hydrodynamic force vector at the i,jth cell centre on the body surface.
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 secondorder Runge–Kutta integration method to obtain the three timevarying 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 bodyfitted 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.
Parameters to characterize flow conditions
The ratio of inertial to viscous forces is expressed by the Reynolds number (Re=UL/v), where U is the speed of the fish relative to the surrounding water, L is the body length of the fish and ν is the kinematic viscosity of water. To derive a meaningful Reynolds number for our computations from the experimental data, we used a timeaveraged speed from stage 3 of the Cstart when the fish began to swim cyclically and speed plateaued. We used two speed averages, one for the COM and one for the tail tip, the body part that usually moves the fastest in an earth frame of reference, to compute an average Reynolds number as follows: (8) where U_{COM} is the final mean speed of the COM, U_{tail} is the final mean speed at the tail tip and Re_{exp} is the experimentally observed Reynolds number used in the CFD computations. We carried out all computations at Re_{exp}=550 (Table 1).
Evaluation of thrust, drag, power and propulsive efficiency
The power required for cyclic swimming, i.e. the work done per unit time, can be calculated as the sum of the work done by the inertial and hydrodynamic forces on the body surface: (9)
Hydrodynamic power was calculated as the sum of the dot products of the velocity and the hydrodynamic forces on the body, such that: (10) where P_{hydro} is the hydrodynamic power, F_{hydro,i,j} is the hydrodynamic force acting on the i,jth surface element and u_{i,j} is the velocity of the i,jth surface element.
Inertial power was computed as: (11) where F_{iner,i,j,k} is the inertial force acting on the i,j,kth element of the body, a_{body,i,j,k} and u_{body,i,j,k} are the acceleration and velocity of the i,j,kth element of the body and m_{i,j,k} is the mass of the i,j,kth body element.
To evaluate the instantaneous mechanical performance of the undulating swimmer, we defined the mechanical propulsive efficiency, i.e. the ratio of the effective work rate and the total power in a timevarying manner, as follows: (12) where T is the thrust vector. Thrust (T) and drag (D) forces were defined relative to the direction of the path of motion as the sum of the positive and negative components of the hydrodynamic forces acting on each surface element.
All forces, moments and powers were further defined as dimensionless coefficients of thrust (T*), drag (D*), net thrust (T_{net}*), moment (N*) and power (P_{total}*) as: (13) where u_{COM} is the COM velocity, calculated by Eqn 5 and the secondorder Runge–Kutta integration method, ρ is the density of water, S is the surface area of the body and L is body length.
Kinematic modules to simulate swimming movements
To prescribe the body deformations of the zebrafish model, we developed two approaches: an analyticalfunction module and a playback module. The analyticalfunction module described the fish's body wave using sinusoidal functions to define a propagating wave from the snout down to the tail: (14) where t is time, x_{axis} is the longitudinal coordinate along the fish axis (x_{axis}=0 at the head and 1 at the tail tip), a(x_{axis}) is the amplitude envelope at x_{axis}, λ is the wavelength and z(x_{axis}, t) denotes the resultant lateral excursion at time t. The body wave amplitude envelope a(x_{axis}) was derived from digitised midlines of swimming fish (Fig. 6A). The amplitude function is expressed in the fish frame of reference (Fig. 6B) to remove all external kinematics from the COM and body kinematics.
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 xaxis to obtain an ‘average amplitude envelope’. This average envelope was symmetric and slightly narrower than the maximal envelope (Fig. 6B).
We applied a correction factor to Eqn 14 to ensure that the body length of the fish remained constant, such that (in normalized body coordinates): (15)
The analyticalfunction approach is used commonly in computational modelling and is useful for parameter space mapping. However, more complex timedependent functions are required to model burstandcoast swimming or starts, when the characteristics of the body wave change from one tail beat to the next. To model Cstarts, we developed the playback 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.
RESULTS AND DISCUSSION
Cyclic swimming
Validation of the model: kinematics
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^{−1} 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 tailbeat 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^{−1}) than the observed swimming speed; it was at the lower boundary of the observed values of 18±2L s^{−1} (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 freeswimming model revealed that the shape of the amplitude envelope was the result of fluid–body interaction. Fluid–body interaction caused the fish to sideslip 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 freeswimming dynamics and hydrodynamics during cyclic swimming.
Validation of the model: flow
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 CFDgenerated flows with the experimentally observed PIV patterns. The first step of the validation was the identification of the correct horizontal transect through the 3D CFD flow field (Fig. 8). The following evidence suggested that the PIV flow field transected the larva ventrally near transect CFDC. First, the distribution and size of the areas of elevated vorticity in CFDC 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^{−1}; vorticity peak near the head (red area, labelled a): −55 s^{−1}]. 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 3D flow field was approximately 0.06L (0.34 mm) below the mediofrontal 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 dorsoventral 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 codetermined the flow along the posterior body. Along the posterior body, the model also correctly predicted the undulationbased 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 crosssection of the flow field (mediofrontal plane) (Fig. 9). Each lateral tail beat produced one vortex pair in the horizontal crosssection. 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 timedependent 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 negativepressure regions and the growing crests caused positivepressure 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 negativepressure peaks were located more distally (at the finfold), while the positivepressure peaks were located more centrally (halfway along the dorsoventral extent on the body) (Fig. 10, dorsoventral distributions are only visible in the lateral views, not the dorsal view). Negative and positivepressure 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 highpressure 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 mediorostral 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 medioposterior 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 crosssections (Fig. 11, panels XY0 to XY–2) that illuminate the flow structures surrounding the swimming larva. Two undulationbased jet flows formed at the two crests and troughs (Fig. 11, jet flows 1 and 2). Both in horizontal and vertical crosssections, 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 crosssections, 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).
Vorticity distribution and wake topology
To describe the wake topology, we used iso Qsurfaces. These surfaces helped to identify vortex rings in a wake. Fig. 12 shows the computed stable wake patterns as isosurfaces of the second invariant of the velocity gradient tensor. The Qcriterion was defined as: (16) where Φ and Ω denote the symmetric and asymmetric parts of the velocity gradient, respectively, and • is the Euclidean matrix norm (Hunt et al., 1988). The isosurfaces (all displayed isosurfaces surfaces were generated after nine tail beat cycles, thus time was >9T) were coloured by contours of flow velocity, which made the higher flow speeds of the central jets stand out at the centre of each vortex ring.
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 (Ushaped) 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 doublerow 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 contralateral 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.
Body wave shape and swimming performance
To study the effect of body wave shape on swimming performance, we varied the width of the amplitude envelope between 0.6 and 1.4 times the experimentally observed value during cyclic swimming (Fig. 14). All studied cases reached their final swimming speed after approximately 15 tail beats (Fig. 14A). Our model predicted that increasing the body wave amplitude causes a monotonic increase in swimming speed, power, hydrodynamic efficiency and power requirement per unit distance (P_{R}) (Fig. 14). P_{R} was defined as the ratio of nondimensionalized power to nondimensionalized speed: (17)
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.
Spontaneous Cstart
Validation of the model: kinematics
In the swimming event shown in Fig. 15, a 5 d.p.f. zebrafish larva spontaneously initiated a swimming bout with a turn (spontaneous Cstart). During the first tail beat (preparatory stroke), the larva bent its entire body into a Cshape. During the subsequent tail beat (propulsive stroke), the larva accelerated forward. Following the Cstart, 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 playback 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 Cstart: the computation predicted a 41 deg angle compared with 44 deg in the experiment; angular speed peaked at 2.1 deg ms^{−1} after 9 ms in the experiment, compared with 2.8 deg ms^{−1} 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^{−1} after 39 ms; Fig. 17A). The propulsive stroke of the Cstart (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.
Validation of the model: flow
The computational model successfully predicted the flow patterns and vorticity levels generated during a Cstart (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 mediofrontal plane of the wake. In this crosssectional 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 crosssections through 3D vortex rings (see below).
Pressure distribution and hydrodynamics of the flow during a Cstart
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 midbody (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 anticlockwise. 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 midbody generated an opposite jet flow (Fig. 19, jet 2) that was most prominent in the trough of the Cbend. 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 anticlockwise 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.
Vorticity distribution and wake topology during a Cstart
Fig. 20 shows the 3D vortex pattern as isosurfaces 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 midbody jet and vortex ring. Consistent with the angular and translational accelerations occurring during the generation of these vortex rings, the jet of the preparatorystroke ring pointed slightly sideways and forward (Fig. 20); the jet of the propulsivestroke 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.
ACKNOWLEDGEMENTS
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.
APPENDIX
Effects of tethering on swimming performance and flow
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 sideslip 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 freeswimming fish. Our threedegreeoffreedom freeswimming 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^{−1}, a fraction of the observed swimming speed (18L s^{−1}).
Models that tether the fish can obtain more realistic results by using an amplitudeenvelope 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 freeswimming 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 amplitudeenvelope function that was derived from the fish midlines in a fish frameofreference with its point of origin at the snout tip. In this way, the model achieved a swimming speed of 16 m s^{−1}, a value very close to the experimentally observed value of 18 m s^{−1}.
FOOTNOTES

FUNDING
This work was partly supported by the Precursory Research for Embryonic Science and Technology, Japan Science and Technology Agency program and a GrantinAid for Scientific Research (nos 18656056 and 18100002) from the Japan Society for the Promotion of Science (JSPS). G.L. was funded by a Japanese Government (MEXT) Scholarship, H.L. was funded by a CheungKong Scholars Programme fellowship, and U.K.M. was funded by the Netherlands Organisation for Scientific Research (NWOALW grant 814.02.006), a JSPS fellowship and a National Science Foundation award (DBI0821820).

Supplementary material available online at http://jeb.biologists.org/cgi/content/full/215/22/4015/DC1
LIST OF SYMBOLS AND ABBREVIATIONS
 a(a_{COM,X}, a_{COM,Y})
 translational acceleration
 a(x_{axis})
 body wave amplitudeenvelope function at x_{axis}
 BCF
 body and caudal fin
 COM
 centre of mass
 D
 drag
 D*
 dimensionless coefficient of drag
 d.p.f.
 days postfertilization
 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
 intergrid boundary cell
 I_{XZ}
 moment of inertia about the yaxis
 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 yaxis
 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 yaxis
 ξ, ζ, η
 computational grid coordinate system
 © 2012.