## 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 3-D 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 free-swimming dynamics from prescribed changes in body shape. We simulated cyclic swimming and a spontaneous C-start. A rigorous comparison with 2-D 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 non-linear 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 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^{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 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.

## 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 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.1*L*√(*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.

### Solution process to model free swimming

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.

### Computational fluid dynamic model

#### Governing equations

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.

*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 pseudo-compressibility coefficient λ. The term

**q**associated with pseudo time was designed for an inner-iteration 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: 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, boundary-fitted and cell-centred storage architecture, we rewrote Eqn 2 into a semi-discrete form: 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*is the volume of the cell (

_{ijk}*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

**q**

_{f}–

**q**) to the semi-discrete Eqn 3: 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 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.

### Computational swimming dynamic model

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.

**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

*y*-axis, α is its global angular acceleration about the

*y*-axis,

*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: where

*A*is the surface area of the

_{i,j}*i,j*th 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.

*N*_{XZ}was calculated as the sum of the cross-product of the positional vector and the force at each computational cell centre on the body surface, such that: where

**r**

*is the positional vector from the COM to the*

_{i,j}*i,j*th cell centre on the body surface and

**F**

*is the hydrodynamic force vector at the*

_{i,j}*i,j*th 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 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.

### Parameters to characterize flow conditions

*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 time-averaged speed from stage 3 of the C-start 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: 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

*P*

_{hydro}is the hydrodynamic power,

**F**

_{hydro,i,j}is the hydrodynamic force acting on the

*i,j*th surface element and

**u**

*is the velocity of the*

_{i,j}*i,j*th surface element.

**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.

**T**

***), drag (

**D**

***), net thrust (

**T**

_{net}

***), moment (

*N**) and power (

*P*

_{total}

***) as: where

**u**

_{COM}is the COM velocity, calculated by Eqn 5 and the second-order 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

*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.1*L* 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.

## 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±2*L* 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 tail-beat frequency (57±8 Hz, *N*=8) and amplitude (0.16±0.025*L, 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 (16*L* s^{−1}) than the observed swimming speed; it was at the lower boundary of the observed values of 18±2*L* 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 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.27*L versus* 0.16*L*), the snout amplitude was, by definition, wider (0*L versus* 0.04*L*), 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.

#### 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 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.3*L* long and 0.1*L* 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 3-D flow field was approximately 0.06*L* (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 1*L*. 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 1*L*, 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/6*T* to 3/6*T*). 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/6*T* to 5/6*T*). 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/6*T* and 1/6*T*). 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).

#### Vorticity distribution and wake topology

*Q*-surfaces. These surfaces helped to identify vortex rings in a wake. Fig. 12 shows the computed stable wake patterns as iso-surfaces of the second invariant of the velocity gradient tensor. The

*Q*-criterion was defined as: where Φ and Ω denote the symmetric and asymmetric parts of the velocity gradient, respectively, and ||•|| is the Euclidean matrix norm (Hunt et al., 1988). The iso-surfaces (all displayed iso-surfaces surfaces were generated after nine tail beat cycles, thus time was >9

*T*) 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/5*T* and rings B and C at 3/5*T*), 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/5*T*–3/5*T* and rings B and C at 3/5*T*–5/5*T*). 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.

### Body wave shape and swimming performance

*P*

_{R}) (Fig. 14).

*P*

_{R}was defined as the ratio of non-dimensionalized power to non-dimensionalized speed:

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 C-start

#### 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 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^{−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.0*L* s^{−1} 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.

#### Validation of the model: flow

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).

#### Pressure distribution and hydrodynamics of the flow during a C-start

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.

#### Vorticity distribution and wake topology during a C-start

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.

## 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 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.16*L* in the tethered fish instead of 0.09*L* 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 6*L* s^{−1}, a fraction of the observed swimming speed (18*L* s^{−1}).

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^{−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 Grant-in-Aid 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 Cheung-Kong Scholars Programme fellowship, and U.K.M. was funded by the Netherlands Organisation for Scientific Research (NWO-ALW grant 814.02.006), a JSPS fellowship and a National Science Foundation award (DBI-0821820).

## LIST OF SYMBOLS AND ABBREVIATIONS

**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