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 101 to 103 (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).

Fig. 1.

Zebrafish larvae and their corresponding body surface models. From top to bottom: larva at 3 days post-fertilization (d.p.f.), model larva at 3 d.p.f., larva at 5 d.p.f. and model larva at 5 d.p.f. Yolk sac and finfold (skin fold surrounding the posterior part of the body, within which the dorsal, anal and caudal fins develop) of the larva are marked in the upper image.

Fig. 1.

Zebrafish larvae and their corresponding body surface models. From top to bottom: larva at 3 days post-fertilization (d.p.f.), model larva at 3 d.p.f., larva at 5 d.p.f. and model larva at 5 d.p.f. Yolk sac and finfold (skin fold surrounding the posterior part of the body, within which the dorsal, anal and caudal fins develop) of the larva are marked in the upper image.

Fig. 2.

(A) Sketch of the fixed XYZ-frame and the degrees of freedom of the computational zebrafish. Free degrees of freedom: yaw, and forward and lateral translation; locked degrees of freedom: vertical translation, roll and pitch. (B,C) Computational grid systems: (B) body-fitted grid (33×45×20, pink); (C) global grid (43×123×73, cyan, the computational ‘water tank’) with a submerged ‘swimming’ object. Further explanation is given in the main text.

Fig. 2.

(A) Sketch of the fixed XYZ-frame and the degrees of freedom of the computational zebrafish. Free degrees of freedom: yaw, and forward and lateral translation; locked degrees of freedom: vertical translation, roll and pitch. (B,C) Computational grid systems: (B) body-fitted grid (33×45×20, pink); (C) global grid (43×123×73, cyan, the computational ‘water tank’) with a submerged ‘swimming’ object. Further explanation is given in the main text.

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

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

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.

Fig. 3.

Grid assembly. (A) In the global grid, cells whose centre falls inside the moving fish were emptied. Cells that border on empty cells were marked as inter-grid boundary cells (IGBCs). (B) In the body-fitted grid, all cells on the outer boundary of this grid were marked as IGBCs.

Fig. 3.

Grid assembly. (A) In the global grid, cells whose centre falls inside the moving fish were emptied. Cells that border on empty cells were marked as inter-grid boundary cells (IGBCs). (B) In the body-fitted grid, all cells on the outer boundary of this grid were marked as IGBCs.

Fig. 4.

Flow diagram of the solution procedure that occurred at each time step.

Fig. 4.

Flow diagram of the solution procedure that occurred at each time step.

The governing equations are the 3-D, incompressible and unsteady NS equations written in a strong conservative form for mass and momentum:
formula
(1)
where
formula
formula
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 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:
formula
(2)
where f=(F+Fv, G+Gv, H+Hv); S(t) denotes the surface of the control volume; n=(nX, nY, nZ) are the components of the unit outward normal vector corresponding to all the surfaces of a polyhedral cell; and ug 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:
formula
(3)
where
formula
for example,
formula
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 Vijk 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).
Fig. 5.

Sketch of the surface area vectors in ξ, ζ and η-directions, the local velocity of the moving cell surface (ug) and volume (Vi,j,k) of cells i, j and k.

Fig. 5.

Sketch of the surface area vectors in ξ, ζ and η-directions, the local velocity of the moving cell surface (ug) and volume (Vi,j,k) of cells i, j and k.

Table 1.

Swimming parameters characterizing cyclic swimming

Swimming parameters characterizing cyclic swimming
Swimming parameters characterizing cyclic swimming

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 σ(qfq) to the semi-discrete Eqn 3:
formula
(4)
where qf 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=qf, 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).
Fig. 6.

Kinematic input data for the cyclic swimming simulation. (A) Body axes for nine successive tail beats derived from a high-speed video of 1000 frames s−1 [cf. fig. 4 in Müller et al. (Müller et al., 2008)]. The only midlines shown are those at the moment when the tail tip reaches a lateral extreme. (B) Amplitude envelope of the body wave. Shown are all midlines obtained for this swimming bout (black lines), relative to the head of the fish. We sampled the envelope at 11 equidistant points (red diamonds) to obtain an average envelope (red lines) as input for the cyclic swimming computations.

Fig. 6.

Kinematic input data for the cyclic swimming simulation. (A) Body axes for nine successive tail beats derived from a high-speed video of 1000 frames s−1 [cf. fig. 4 in Müller et al. (Müller et al., 2008)]. The only midlines shown are those at the moment when the tail tip reaches a lateral extreme. (B) Amplitude envelope of the body wave. Shown are all midlines obtained for this swimming bout (black lines), relative to the head of the fish. We sampled the envelope at 11 equidistant points (red diamonds) to obtain an average envelope (red lines) as input for the cyclic swimming computations.

Fig. 7.

Computed speed of the centre of mass of the fish model for two amplitude envelopes (average envelope: red line; 1.1 × average envelope: purple line). The computation was started with the centre of mass at rest. Once a constant tail-beat average swimming speed was established after nine tail-beat cycles, the average computational speed equalled 16L s−1 (average envelope) and 17.5L s−1 (1.1 × average envelope), which is close to the experimentally observed speed of 18±2L s−1 (black line with grey boundaries) (Müller et al., 2008).

Fig. 7.

Computed speed of the centre of mass of the fish model for two amplitude envelopes (average envelope: red line; 1.1 × average envelope: purple line). The computation was started with the centre of mass at rest. Once a constant tail-beat average swimming speed was established after nine tail-beat cycles, the average computational speed equalled 16L s−1 (average envelope) and 17.5L s−1 (1.1 × average envelope), which is close to the experimentally observed speed of 18±2L s−1 (black line with grey boundaries) (Müller et al., 2008).

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.

Table 2.

Comparison of key parameters between experiment and CFD predictions

Comparison of key parameters between experiment and CFD predictions
Comparison of key parameters between experiment and CFD predictions

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.

Fig. 8.

Comparison of computational fluid dynamics (CFD) and experimentally observed vorticity patterns. A comparison of horizontal sections through the CFD vorticity patterns at four different dorso-ventral levels within the fish larva (A–D) (for exact dorso-ventral position see inset) shows that the experimentally obtained flow field (E) does not transect the larva in the medio-frontal plane, but slightly below.

Fig. 8.

Comparison of computational fluid dynamics (CFD) and experimentally observed vorticity patterns. A comparison of horizontal sections through the CFD vorticity patterns at four different dorso-ventral levels within the fish larva (A–D) (for exact dorso-ventral position see inset) shows that the experimentally obtained flow field (E) does not transect the larva in the medio-frontal plane, but slightly below.

The deformable-body 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:
formula
(5)
where F is the fluid force vector acting on the body COM, uCOM=(uCOM,X, uCOM,Y) is the translational velocity of the COM of the body, aCOM 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, IXZ is the yawing moment of inertia (which depends on the instantaneous body shape) and NXZ 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 (Fx,Fz) 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:
formula
(6)
where Ai,j is the surface area of the i,jth cell on the body surface, and Pbody,i,j and Sbody,i,j are pressure and shear stress vectors on this surface cell, respectively, which arose from the solution of NS equations.
Fig. 9.

Comparison of CFD and experimental flow patterns generated during cyclic swimming. The computational results are shown in the medio-frontal plane (left column) and 0.6L below the medio-frontal plane (centre column), which corresponds closely with the position of the experimental PIV plane (right column). Slight differences arose between CFD and PIV results because the actual fish made a weak turning manoeuvre whereas the modelled fish swam perfectly cyclically (symmetric amplitude envelope, identical tail beats). Furthermore, PIV is less able than CFD to capture strong gradients.

Fig. 9.

Comparison of CFD and experimental flow patterns generated during cyclic swimming. The computational results are shown in the medio-frontal plane (left column) and 0.6L below the medio-frontal plane (centre column), which corresponds closely with the position of the experimental PIV plane (right column). Slight differences arose between CFD and PIV results because the actual fish made a weak turning manoeuvre whereas the modelled fish swam perfectly cyclically (symmetric amplitude envelope, identical tail beats). Furthermore, PIV is less able than CFD to capture strong gradients.

The hydrodynamic moment NXZ 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:
formula
(7)
where ri,j is the positional vector from the COM to the i,jth cell centre on the body surface and Fi,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 (aCOM,X, aCOM,Y and α) at each time step. Then we used the second-order Runge–Kutta integration method to obtain the three time-varying velocities (uCOM,X, uCOM,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.

Fig. 10.

Temporal relationship between body deformation, jet flow and surface pressure within one swimming cycle (1/6 ~ 6/6T). Top panels: CFD results; bottom panel: sketch of salient features (S indicates a negative pressure region, and P a positive pressure region). The combination of these pressure regions with the local body orientations generates force vectors with a considerable thrust (i.e. forward) component.

Fig. 10.

Temporal relationship between body deformation, jet flow and surface pressure within one swimming cycle (1/6 ~ 6/6T). Top panels: CFD results; bottom panel: sketch of salient features (S indicates a negative pressure region, and P a positive pressure region). The combination of these pressure regions with the local body orientations generates force vectors with a considerable thrust (i.e. forward) component.

Fig. 11.

Three-dimensional structure of the body flows during cyclic swimming – spatial structure. Twelve vertical and five horizontal cross-sections through the fluid field around the larva are shown plus a sketch of the flow pattern at one moment during the cycle. In each cross-section, the positive direction of vorticity is anti-clockwise in the plane. Jet flows (JFs) on the body are labelled JF1 and JF2, incomplete vortex rings are labelled Ra and Rb. The faded vortex ring behind Rb in the sketch represents a vortex ring shed in the wake. Previously shed vortices are not shown.

Fig. 11.

Three-dimensional structure of the body flows during cyclic swimming – spatial structure. Twelve vertical and five horizontal cross-sections through the fluid field around the larva are shown plus a sketch of the flow pattern at one moment during the cycle. In each cross-section, the positive direction of vorticity is anti-clockwise in the plane. Jet flows (JFs) on the body are labelled JF1 and JF2, incomplete vortex rings are labelled Ra and Rb. The faded vortex ring behind Rb in the sketch represents a vortex ring shed in the wake. Previously shed vortices are not shown.

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 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:
formula
(8)
where UCOM is the final mean speed of the COM, Utail is the final mean speed at the tail tip and Reexp is the experimentally observed Reynolds number used in the CFD computations. We carried out all computations at Reexp=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:
formula
(9)
Hydrodynamic power was calculated as the sum of the dot products of the velocity and the hydrodynamic forces on the body, such that:
formula
(10)
where Phydro is the hydrodynamic power, Fhydro,i,j is the hydrodynamic force acting on the i,jth surface element and ui,j is the velocity of the i,jth surface element.
Inertial power was computed as:
formula
(11)
where Finer,i,j,k is the inertial force acting on the i,j,kth element of the body, abody,i,j,k and ubody,i,j,k are the acceleration and velocity of the i,j,kth element of the body and mi,j,k is the mass of the i,j,kth body element.
Fig. 12.

Three-dimensional structure of the wake during cyclic swimming – temporal pattern over one period T of cyclic swimming. Left: iso-surface of vorticity=40 s−1 and vectors of flow velocity; right: 2-D vorticity in medio-frontal plane (i.e. halfway along the height of the body) and sketch of the jet flows (orange arrows) and vortex rings (orange rings A, B and C; only three rings indicated).

Fig. 12.

Three-dimensional structure of the wake during cyclic swimming – temporal pattern over one period T of cyclic swimming. Left: iso-surface of vorticity=40 s−1 and vectors of flow velocity; right: 2-D vorticity in medio-frontal plane (i.e. halfway along the height of the body) and sketch of the jet flows (orange arrows) and vortex rings (orange rings A, B and C; only three rings indicated).

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 time-varying manner, as follows:
formula
(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 (Tnet*), moment (N*) and power (Ptotal*) as:
formula
(13)
where uCOM 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

To prescribe the body deformations of the zebrafish model, we developed two approaches: an analytical-function module and a play-back module. The analytical-function module described the fish's body wave using sinusoidal functions to define a propagating wave from the snout down to the tail:
formula
(14)
where t is time, xaxis is the longitudinal coordinate along the fish axis (xaxis=0 at the head and 1 at the tail tip), a(xaxis) is the amplitude envelope at xaxis, λ is the wavelength and z(xaxis, t) denotes the resultant lateral excursion at time t. The body wave amplitude envelope a(xaxis) 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 x-axis 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):
formula
(15)

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±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 tail-beat frequency (57±8 Hz, N=8) and amplitude (0.16±0.025L, N=8, measured as lateral excursion of the tail tip from the mean path of motion). Table 1 lists the kinematic parameters used in the simulation, which were based on the analysis of nine consecutive tail beats.

The forward swimming speed of the model was slightly lower (16L s−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 UCOM=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).

Fig. 13.

Overview of wake pattern during cyclic swimming (see supplementary material Movie 1). Iso-surfaces of Q=0.02 are coloured by fluid speed and are semi-transparent. All displayed iso-surfaces were generated after nine tail beat cycles from the start.

Fig. 13.

Overview of wake pattern during cyclic swimming (see supplementary material Movie 1). Iso-surfaces of Q=0.02 are coloured by fluid speed and are semi-transparent. All displayed iso-surfaces were generated after nine tail beat cycles from the start.

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

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

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.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 3-D flow field was approximately 0.06L (0.34 mm) below the medio-frontal plane of the larva, only the yolk sac and the ventral finfold were situated in the light sheet.

Fig. 14.

Predicted effect of body wave amplitude on five different performance indicators of cyclic swimming. The studied body wave amplitudes were 0.6, 0.8, 1.0, 1.2 and 1.4 times the experimentally observed value. (A) A family of swimming speed curves for five wave amplitudes. Higher amplitudes lead to higher average swimming speeds as well as larger cyclic variations in speed. (B–E) Mean speed during cyclic swimming (>9 tail beats after the initiation), mean non-dimensional power, mean efficiency and mean energy cost per unit distance (UDC, Eqn 17) as a function of normalized wave amplitude.

Fig. 14.

Predicted effect of body wave amplitude on five different performance indicators of cyclic swimming. The studied body wave amplitudes were 0.6, 0.8, 1.0, 1.2 and 1.4 times the experimentally observed value. (A) A family of swimming speed curves for five wave amplitudes. Higher amplitudes lead to higher average swimming speeds as well as larger cyclic variations in speed. (B–E) Mean speed during cyclic swimming (>9 tail beats after the initiation), mean non-dimensional power, mean efficiency and mean energy cost per unit distance (UDC, Eqn 17) as a function of normalized wave amplitude.

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.

Fig. 15.

Kinematic input data for C-start simulation. Body axes for seven successive tail beats during a swimming bout that was initiated by a C-start from rest (recorded at 1000 frames s−1). The only midlines shown are those at the moment when the tail tip reaches a lateral extreme (labelled by tail beat number; O and E denote original and end positions, respectively). The table provides the times for each of the selected body axes.

Fig. 15.

Kinematic input data for C-start simulation. Body axes for seven successive tail beats during a swimming bout that was initiated by a C-start from rest (recorded at 1000 frames s−1). The only midlines shown are those at the moment when the tail tip reaches a lateral extreme (labelled by tail beat number; O and E denote original and end positions, respectively). The table provides the times for each of the selected body axes.

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

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

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

Fig. 16.

Comparison of whole-body kinematics between experimental observation (red) and CFD (black). (A) Experimental (red) and computational (black) trajectories of the larva's centre of mass (COM) during the swimming bout of Fig. 15. (B) Comparison of outlines of the actual (grey) and the model fish (red) at the start, at the end of the preparatory stroke (first tail beat) and at the end of the bout.

Fig. 16.

Comparison of whole-body kinematics between experimental observation (red) and CFD (black). (A) Experimental (red) and computational (black) trajectories of the larva's centre of mass (COM) during the swimming bout of Fig. 15. (B) Comparison of outlines of the actual (grey) and the model fish (red) at the start, at the end of the preparatory stroke (first tail beat) and at the end of the bout.

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

To describe the wake topology, we used iso 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:
formula
(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 iso-surfaces (all displayed iso-surfaces 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.
Fig. 17.

Swimming performance parameter for a spontaneous C-start and consequent swimming bout. (A) Computational (red) and experimental (black) speed of the fish's centre of mass. (B) Fluid dynamic power and efficiency of the seven tail beats comprising the swimming bout. (C) Computational (red) and experimental (black) angular speed of the larva's centre of mass. (D) Computational (red) and experimental (black) heading angle of the larva.

Fig. 17.

Swimming performance parameter for a spontaneous C-start and consequent swimming bout. (A) Computational (red) and experimental (black) speed of the fish's centre of mass. (B) Fluid dynamic power and efficiency of the seven tail beats comprising the swimming bout. (C) Computational (red) and experimental (black) angular speed of the larva's centre of mass. (D) Computational (red) and experimental (black) heading angle of the larva.

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

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

Fig. 18.

PIV and CFD vorticity field (colour map) and flow velocity (black arrows) of preparatory and propulsive strokes of a C-start of a 5 d.p.f. zebrafish larva. Major flow features are displayed in the sketches at the left and right hand side. The CFD-predicted flow fields resemble the spatial and temporal patterns in the experimentally observed flow.

Fig. 18.

PIV and CFD vorticity field (colour map) and flow velocity (black arrows) of preparatory and propulsive strokes of a C-start of a 5 d.p.f. zebrafish larva. Major flow features are displayed in the sketches at the left and right hand side. The CFD-predicted flow fields resemble the spatial and temporal patterns in the experimentally observed flow.

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 (PR) (Fig. 14). PR was defined as the ratio of non-dimensionalized power to non-dimensionalized speed:
formula
(17)

Note that PR 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 (PR) 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.

Fig. 19.

Temporal relationship between body deformation, jet flow, surface pressure and shear flow within a C-start. Each panel shows pressure (Lp, Rp) and shear stresses (Ls, Rs) acting on the body for the left and right side of the body; in the lateral views, the larva's body is shown straight for clarity. Corresponding regions are indicated by stippled lines (red: positive forces; blue: negative forces). The C-start generated strong jets (black arrows). Shear stresses are much lower than pressure forces; note the scale bars for both at the left edge of the figure. Rostrad shear force vectors are defined positive while caudad vectors are negative. Vertical shear force components are not included in this figure.

Fig. 19.

Temporal relationship between body deformation, jet flow, surface pressure and shear flow within a C-start. Each panel shows pressure (Lp, Rp) and shear stresses (Ls, Rs) acting on the body for the left and right side of the body; in the lateral views, the larva's body is shown straight for clarity. Corresponding regions are indicated by stippled lines (red: positive forces; blue: negative forces). The C-start generated strong jets (black arrows). Shear stresses are much lower than pressure forces; note the scale bars for both at the left edge of the figure. Rostrad shear force vectors are defined positive while caudad vectors are negative. Vertical shear force components are not included in this figure.

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.0L 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.

Fig. 20.

Three-dimensional structure of the wake during a spontaneous C-start – temporal pattern (see supplementary material Movie 2). Left column: iso-surface of Q=0.02 and flux vectors 30, 40 and 50 ms after the initiation of the C-start, coloured by fluid speed. Right column: accordant 2-D vorticity in the medio-frontal plane crosses the middle height of the body, and sketch of jet flows (orange arrows) and vortex rings (orange rings).

Fig. 20.

Three-dimensional structure of the wake during a spontaneous C-start – temporal pattern (see supplementary material Movie 2). Left column: iso-surface of Q=0.02 and flux vectors 30, 40 and 50 ms after the initiation of the C-start, coloured by fluid speed. Right column: accordant 2-D vorticity in the medio-frontal plane crosses the middle height of the body, and sketch of jet flows (orange arrows) and vortex rings (orange rings).

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.16L in the tethered fish instead of 0.09L in the free-swimming fish. Our three-degree-of-freedom free-swimming CFD model allows lateral forces and moment to act on the fish, and we observe both side slip and yaw, which reduces the effective lateral excursions of the tail in an earth frame of reference. As shown in supplementary material Fig. S1, allowing slip results in a large reduction of propulsion (supplementary material Fig. S1C): the cyclic swim CFD speed is only 6L s−1, a fraction of the observed swimming speed (18L 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(aCOM,X, aCOM,Y)

    translational acceleration

  •  
  • a(xaxis)

    body wave amplitude-envelope function at xaxis

  •  
  • 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(FX, FY,, FZ)

    force vector

  •  
  • G(GX, GY,, GZ)

    force vector

  •  
  • H(HX, HY,HZ)

    force vector

  •  
  • IGBC

    inter-grid boundary cell

  •  
  • IXZ

    moment of inertia about the y-axis

  •  
  • JF

    jet flow

  •  
  • L

    body length

  •  
  • M

    body mass of fish

  •  
  • mi,j,k

    mass of a body element

  •  
  • n

    unit normal vector

  •  
  • NS

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

  •  
  • NXZ

    hydrodynamic moment about the y-axis

  •  
  • p

    pressure

  •  
  • P*total

    dimensionless coefficient of power

  •  
  • Pbody

    pressure acting on the body

  •  
  • Phydro

    hydrodynamic power

  •  
  • Piner

    inertial power

  •  
  • PIV

    particle image velocimetry

  •  
  • Ptotal

    total power

  •  
  • r

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

  •  
  • Re

    Reynolds number

  •  
  • Reexp

    Reynolds number used in computation

  •  
  • S

    body surface area

  •  
  • Sbody

    shear stress acting on body

  •  
  • Si,j

    area of surface element

  •  
  • t

    time

  •  
  • T

    thrust

  •  
  • T*

    dimensionless coefficient of thrust

  •  
  • T*net

    dimensionless coefficient of net thrust

  •  
  • Tnet

    net thrust

  •  
  • u

    velocity at the surface of a grid cell

  •  
  • uCOM(uCOM,X, uCOM,Y, uCOM,Z)

    translational velocity vector of the COM

  •  
  • V(u, v, w)

    fluid velocity vector

  •  
  • VR

    vortex ring

  •  
  • xaxis

    length along the body axis

  •  
  • X, Y, Z

    world coordinate system

  •  
  • (XCOM, YCOM, ZCOM)

    position vector of the COM

  •  
  • z(xaxis, 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

REFERENCES

Blake
R. W.
(
1983
).
Fish locomotion
.
Cambridge
:
Cambridge University Press
.
Blickhan
R.
,
Krick
C.
,
Zehren
D.
,
Nachtigall
W.
(
1992
).
Generation of a vortex chain in the wake of a subundulatory swimmer
.
Naturwissenschaften
79
,
220
-
221
.
Borazjani
I.
,
Sotiropoulos
F.
(
2008
).
Numerical investigation of the hydrodynamics of carangiform swimming in the transitional and inertial flow regimes
.
J. Exp. Biol.
211
,
1541
-
1558
.
Borazjani
I.
,
Sotiropoulos
F.
(
2009
).
Numerical investigation of the hydrodynamics of anguilliform swimming in the transitional and inertial flow regimes
.
J. Exp. Biol.
212
,
576
-
592
.
Borazjani
I.
,
Sotiropoulos
F.
(
2010
).
On the role of form and kinematics on the hydrodynamics of self-propelled body/caudal fin swimming
.
J. Exp. Biol.
213
,
89
-
107
.
Carling
J.
,
Williams
T. L.
,
Bowtell
G.
(
1998
).
Self-propelled anguilliform swimming: simultaneous solution of the two-dimensional navier-stokes equations and Newton's laws of motion
.
J. Exp. Biol.
201
,
3143
-
3166
.
Domenici
P.
,
Blake
R. W.
(
1991
).
The kinematics and performance of the escape response in the angelfish (Pterophyllum eimekei)
.
J. Exp. Biol.
156
,
187
-
205
.
Domenici
P.
,
Blake
R. W.
(
1997
).
The kinematics and performance of fish fast-start swimming
.
J. Exp. Biol.
200
,
1165
-
1178
.
Flammang
B. E.
,
Lauder
G. V.
,
Troolin
D. R.
,
Strand
T. E.
(
2011
).
Volumetric imaging of fish locomotion
.
Biol. Lett.
7
,
695
-
698
.
Gazolla
M.
,
Van Rees
W. M.
,
Koumoutsakos
P.
(
2012
).
C-start: optimal start of larval fish
.
J. Fluid Mech.
698
,
5
-
18
.
Gen
L.
,
Liu
H.
,
Müller
U. K.
,
van Leeuwen
J. L.
(
2011
).
Swimming hydrodynamics and maneuverability in C-start of zebrafish larvae: an integrated computational study
.
ASME-JSME-KSME 2011 Joint Fluids Engineering Conference
, Vol.
1
,
AJK2011-19020
, pp.
2059
-
2066
.
Hunt
J. C. R.
,
Wray
A. A.
,
Moin
P.
(
1988
).
Eddies, streams, and convergence zones in turbulent flows
. In
Studying Turbulence Using Numerical Simulation Databases, 2. Proceedings of the 1988 Summer Program
,
SEE N89-24538 18-34
, pp.
193
-
208
.
Hunter
J. R.
(
1972
).
Swimming and feeding behavior of larval anchovy Engraulis mordax
.
Fish Bull.
70
,
821
-
838
.
Katumata
Y.
,
Müller
U. K.
,
Liu
H.
(
2009
).
Computation of self-propelled swimming in larva fishes
.
J. Biomech. Eng.
4
,
54
-
66
.
Kern
S.
,
Koumoutsakos
P.
(
2006
).
Simulations of optimized anguilliform swimming
.
J. Exp. Biol.
209
,
4841
-
4857
.
Lauder
G. V.
(
2011
).
Swimming hydrodynamics: ten questions and the technical approaches needed to resolve them
.
Exp. Fluids
51
,
23
-
35
.
Lauder
G. V.
,
Drucker
E. G.
(
2002
).
Forces, fishes, and fluids: hydrodynamic mechanisms of aquatic locomotion
.
News Physiol. Sci.
17
,
235
-
240
.
Lauder
G. V.
,
Tytell
E. D.
(
2006
).
Hydrodynamics of undulatory propulsion
. In
Fish Biomechanics
(ed.
Shadwick
R. E.
,
Lauder
G. V.
), pp.
425
-
468
.
San Diego, CA
:
Academic Press
.
Liu
H.
(
2005
).
Simulation-based biological fluid dynamics in animal locomotion
.
Appl. Mech. Rev.
58
,
269
-
282
.
Liu
H.
(
2009
).
Integrated modeling of insect flight: from morphology, kinematics to aerodynamics
.
J. Comput. Phys.
228
,
439
-
459
.
Liu
H.
,
Kawachi
K.
(
1999
).
A numerical study of undulatory swimming
.
J. Comput. Phys.
155
,
223
-
247
.
Liu
H.
,
Wassersug
R.
,
Kawachi
K.
(
1996
).
A computational fluid dynamics study of tadpole swimming
.
J. Exp. Biol.
199
,
1245
-
1260
.
Liu
H.
,
Wassersug
R.
,
Kawachi
K.
(
1997
).
The three-dimensional hydrodynamics of tadpole locomotion
.
J. Exp. Biol.
200
,
2807
-
2819
.
McHenry
M. J.
,
Lauder
G. V.
(
2005
).
The mechanical scaling of coasting in zebrafish (Danio rerio)
.
J. Exp. Biol.
208
,
2289
-
2301
.
Müller
U. K.
(
2008
).
Swimming and muscle
. In
Fish Larval Physiology
(ed.
Finn
R. N.
,
Kapoor
B. G.
), pp.
523
-
549
.
Enfield, NH
:
Science Publishers
.
Müller
U. K.
,
van Leeuwen
J. L.
(
2004
).
Swimming of larval zebrafish: ontogeny of body waves and implications for locomotory development
.
J. Exp. Biol.
207
,
853
-
868
.
Müller
U. K.
,
van den Heuvel
B. L. E.
,
Stamhuis
E. J.
,
Videler
J. J.
(
1997
).
Fish foot prints: morphology and energetics of the wake behind a continuously swimming mullet (Chelon labrosus Risso)
.
J. Exp. Biol.
200
,
2893
-
2906
.
Müller
U. K.
,
Stamhuis
E. J.
,
Videler
J. J.
(
2000
).
Hydrodynamics of unsteady fish swimming and the effects of body size: comparing the flow fields of fish larvae and adults
.
J. Exp. Biol.
203
,
193
-
206
.
Müller
U. K.
,
Smit
J.
,
Stamhuis
E. J.
,
Videler
J. J.
(
2001
).
How the body contributes to the wake in undulatory fish swimming: flow fields of a swimming eel (Anguilla anguilla)
.
J. Exp. Biol.
204
,
2751
-
2762
.
Müller
U. K.
,
Stamhuis
E. J.
,
Videler
J. J.
(
2002
).
Riding the waves: the role of the body wave in undulatory fish swimming
.
Integr. Comp. Biol.
42
,
981
-
987
.
Müller
U. K.
,
van den Boogaart
J. G. M.
,
van Leeuwen
J. L.
(
2008
).
Flow patterns of larval fish: undulatory swimming in the intermediate flow regime
.
J. Exp. Biol.
211
,
196
-
205
.
Nauen
J. C.
,
Lauder
G. V.
(
2002
).
Hydrodynamics of caudal fin locomotion by chub mackerel, Scomber japonicus (Scombridae)
.
J. Exp. Biol.
205
,
1709
-
1724
.
Prewitt
N. C.
,
Belk
D. M.
,
Shyy
W.
(
2000
).
Parallel computing of overset grids for aerodynamic problems with moving objects
.
Prog. Aerosp. Sci.
36
,
117
-
172
.
Reid
D. A. P.
,
Hildenbrandt
H.
,
Padding
J. T.
,
Hemelrijk
C. K.
(
2012
).
Fluid dynamics of moving fish in a two-dimensional multiparticle collision dynamics model
.
Phys. Rev. E
85
,
021901
.
Tytell
E. D.
,
Lauder
G. V.
(
2004
).
The hydrodynamics of eel swimming: I. Wake structure
.
J. Exp. Biol.
207
,
1825
-
1841
.
Tytell
E. D.
,
Lauder
G. V.
(
2008
).
Hydrodynamics of the escape response in bluegill sunfish, Lepomis macrochirus
.
J. Exp. Biol.
211
,
3359
-
3369
.
Videler
J. J.
(
1993
).
Fish Swimming
.
London
:
Chapman & Hall
.
Wardle
C. S.
,
Videler
J. J.
,
Altringham
J. D.
(
1995
).
Tuning in to fish swimming waves: body form, swimming mode and muscle function
.
J. Exp. Biol.
198
,
1629
-
1636
.
Webb
P. W.
(
1976
).
The effect of size on the fast-start performance of rainbow trout Salmo cairdneri, and a consideration of piscivorous predator–prey interactions
.
J. Exp. Biol.
65
,
157
-
177
.
Weihs
D.
(
1973
).
The mechanism of rapid starting of slender fish
.
Biorheology
10
,
343
-
350
.
Wolfgang
M. J.
,
Anderson
J. M.
,
Grosenbaugh
M. A.
,
Yue
D. K. P.
,
Triantafyllou
M. S.
(
1999
).
Near-body flow dynamics in swimming fish
.
J. Exp. Biol.
202
,
2303
-
2327
.