Large-amplitude undulatory fish swimming: fluid mechanics coupled to internal mechanics

This paper will be not so much a complete review of the whole subject as a summary of the basic physical principles of undulatory swimming together with a more detailed survey of some recent research.

It is worth bearing in mind that Newton’s laws of motion apply to a swimming organism as a whole as well as to individual elements of it or of the water. Thus, the total force exerted by the water on the body must be balanced by the body’s inertia, i.e. rate of change of linear momentum, at every instant (neglecting buoyancy forces for now). Also, the total torque must balance the rate of change of angular momentum.

Here, we concentrate on fish (and cetaceans) that pass a wave of curvature backwards along the body (Fig. 1). The amplitude of the wave usually increases with distance from the nose. In all cases, the Reynolds number is large, and it is commonly assumed that it is therefore possible to consider the forces exerted by the water on the fish in two separate parts. There is the drag, D, exerted by viscous shear stresses in the boundary layers on the body (usually assumed not to separate, but there is an additional contribution to the drag from the pressure if the boundary layers do separate), which tends to slow the animal down. Then there is the thrust, T, generated by the body (and/or tail) motions, which pushes it along. The thrust comes from the reactive force, the difference in pressure between the two sides of the fish as it accelerates the fluid around it, and (sometimes) from the lift force on a lunate tail. At high Reynolds number, viscous forces make a negligible contribution to the thrust, except in so far that it is viscosity that forces flow to come tangentially off the sharp trailing edge of the caudal fin (the so-called Kutta condition) or, equivalently, that forces the pressure difference between the two sides to be zero at the trailing edge.

Lighthill’s theory was based on several simplifying assumptions: a body shape that varies slowly with s; small curvature of the centre surface; and neglect of any effect the vortex wake might have on the pressure distribution on the body (the fact that the flow comes smoothly off the trailing edge leads inevitably to a wake which rapidly distorts into a number of concentrated vortex structures; see below). These assumptions have been relaxed in more recent computational work by Cheng et al. (1991), although at the expense of assuming the fish to be an infinitely thin rectangular ‘waving plate’. Their model was designed to incorporate the effect of the wake on the distribution of pressure force P. They restricted attention to small amplitudes of undulation, with the consequence that the wake could be assumed to remain approximately planar. Their method was based on the recognition that a thin solid body can be replaced, for computational purposes, by an array of rectangular vortex ‘rings’ whose edges form a lattice that covers the surface of the body. At each time step, a row of vortex rings is shed from the trailing edge (caudal fin), and the strengths of these vortex rings are chosen so that the pressures on the two sides of the body, at the trailing edge, are the same; as stated above, this is equivalent to the Kutta condition, and means that P is zero at the trailing edge, not a maximum as in the Lighthill theory. After being shed, the wake vortices are carried downstream (relative to the waving plate) at the steady oncoming flow speed U₀. The strengths of the vortex rings that represent the body itself are chosen (at each time step) to make sure that the appropriate boundary conditions are satisfied on the body: i.e. that the component of fluid velocity perpendicular to the boundary at any boundary point is equal to the prescribed perpendicular component of the boundary velocity. This theory, like Lighthill’s, neglects viscous effects, assuming that the Reynolds number is large and that the boundary layer is thin and does not separate.

The two-dimensional version of this vortex-ring panel (or vortex lattice) method is particularly simple. Here, dorso-ventral variation is ignored, and the waving plate behaves as if it were infinitely deep. The vortex rings become line (or ‘point’) vortices, and the theory becomes equivalent to the analytical two-dimensional theory of Wu (1961); agreement with Wu’s predictions, e.g. of P(t), is a partial test of the computational method, a test that the method of Cheng et al. (1991) passed successfully.

A large-amplitude version of the vortex-lattice method has been developed, still for an infinitely thin rectangular ‘fish’, by Hill (1998). Some of its results will be outlined here. The computations are much more difficult in this case, because each vortex ring in the wake consists of a loop of vortex filaments which remain attached to the same elements of fluid all the time (when viscosity is not important), and the fluid elements have components of velocity perpendicular to the swimming direction because of the effect of the other vortex filaments. Thus, every vortex filament has to be tracked as it moves in the velocity field generated by all the others. The strengths of the vortex filaments in the wake remain constant as they move around. The vortex filaments on the body, in contrast, have unknown strengths, although their positions are known.

Nevertheless, the numerical method is feasible, and Hill has worked out the details in both the two-and three-dimensional cases, although at the time of writing he had not succeeded in implementing the recoil correction in three dimensions and the pressure difference P does not go precisely to zero at s=l (Fig. 2). The Lighthill theory for the same case gives results of similar magnitude which increase towards the tail. In two dimensions, good agreement is obtained with the results of both Cheng et al. (1991) and Wu (1961) at small amplitude. An example of the two-dimensional developing wake is shown in Fig. 3; it may be noted that the vortices representing the wake begin to cluster together and rotate about each other at 1–2 body lengths after being shed. This gradual rolling up of the trailing vortex sheet is in contrast to the shedding of already isolated vortices that is observed behind (three-dimensional) fish that swim with a rapid rotation of the caudal fin at the end of each tail-beat, which cannot be represented by equation 5 with the chosen parameter values (Müller et al., 1997; Videler et al., 1999). In three dimensions, good agreement is obtained with the semi-analytical theory of Chopra (1974) for a pitching and heaving flat plate. An example of the developing wake is shown in Fig. 4. We may note how the vortices roll up at the top (and bottom) edge(s) of the wake as well as across it. [For details of the numerical method and its testing, the reader is referred to Hill’s (1998) thesis].

Müller et al. (1997) and Videler et al. (1999) have recently published particle image velocimetry (PIV) measurements of the approximately two-dimensional velocity field in the plane of symmetry of two species of swimming fish, namely mullet (Chelon labrosus) and eel (Anguilla anguilla). A corresponding velocity field as computed by Hill’s model (for the eel) is shown in Fig. 5, where it can be seen that the sideways pushing and pulling of the fluid by the lateral undulations of the body generates apparent curved streamlines which, if they were not interrupted by the waving plate, would resemble the cross section of a concentrated vortex. Similar velocity fields are revealed by the pictures of Videler et al. (1999). These authors termed this flow structure the ‘body vortex’, which they stated would be shed, forming one of the concentrated wake vortices, when it reached the trailing edge. This interpretation is somewhat hard to follow, since it is not clear how the irrotational motion of Fig. 5 will be ‘shed’ as a rotational three-dimensional vortex. This fundamental question is one of great current interest.

A more complete three-dimensional ‘panel method’ simulation of fish swimming than that described above, including non-zero body thickness and one or two projecting fins, has been developed by Wolfgang et al. (1999). In their model, the fish body is replaced by a surface distribution of sources and dipoles, and only the thin fins are replaced by vortex lattices. As yet, few numerical details are available, and the recoil correction is not fully implemented, but if one takes the authors’ assurances about their code’s accuracy at face value, it is clear that this group has developed inviscid modelling of fish swimming further than anyone else. Like Müller et al. (1997), they devote considerable space to a discussion of ‘body vortices’ as a source of wake vortices, but the shedding process remains unclear.

All the above models of fish swimming have been essentially inviscid, except that the process of vortex shedding at the tail requires there to be viscosity. Hill (1998) is perhaps the first who has used the output of such a calculation as input to computing the velocity distribution in the boundary layers on the two sides of the fish and, hence, the viscous drag. As stated above, in the context of equation 2, that is a necessary prerequisite to calculating the time-dependent swimming speed, U(t), and checking whether the assumed average swimming speed, U₀, is indeed self-consistent. It is also necessary for a correct determination of internal mechanics at large amplitude (see below). So far, only the two-dimensional boundary-layer equations have been solved, revealing a moderate (30 %) mismatch between mean drag and mean thrust, indicating that the assumed value of U₀ needs to be modified, but there is nothing like the factor of four mismatch estimated by Lighthill (1971) on the basis of drag measurements on towed, inert fish. The limitation on the boundary-layer solver is that it breaks down as soon as reversed flow appears anywhere, a necessary (but not sufficient) precursor of boundary-layer separation. The calculations were therefore limited to moderate caudal fin amplitudes of less than 10 % of body length.

The ultimate in accurate computation of fish swimming should be the numerical solution of the full, three-dimensional Navier–Stokes equations, which would incorporate all fluid mechanical effects, given the fish kinematics as input. This has been tried by Liu et al. (1996), in a two-dimensional simulation of tadpole swimming; they used a new computational fluid dynamics (CFD) code, but did not give enough computational detail to make it easy to assess accuracy, especially at the high values used for the Reynolds number based on body length and swimming speed, up to 10⁵; Carling et al. (1998) are also developing such a code, but so far they too are restricted to two dimensions and Re up to 4000. The benefits of a Navier–Stokes solver will be that all fluid mechanical features will be included at once, so there will be no artificial distinction between (for example) thrust and drag forces, and the swimming speed U(t) will be computed as part of the solution. The disadvantage is the need for fine resolution of all the boundary layers, wakes and vortices, which requires huge computational resources at all but unrealistically small values of the Reynolds number. For now, an inviscid panel or vortex method plus boundary layer theory is likely to be the best approach to fish hydrodynamics.

An understanding of the hydrodynamics of fish swimming is not an end in itself for biologists, who are more interested in how the fish achieve their observed body motions and swimming speed, i.e. what their swimming muscles actually do. The rate at which a muscle fibre contracts depends upon its length (i.e. the current geometry), its state of activation, and the load against which it is contracting. In the case of fish muscle, the load comes primarily from the hydrodynamic forces, plus the body inertia and the stresses required to deform the body tissues. For undulatory swimming, analysis of this load is best achieved by treating the fish body as an active bending beam (Wu, 1961; Hess and Videler, 1984; Cheng et al., 1998). The equations of motion, relating the rates of change of longitudinal and transverse momentum to those components of force, and rate of change of angular momentum to the torque, are written down for a slice of the fish, of thickness δs, perpendicular to its spine (Fig. 6). These equations couple the external forces per unit length (pressure force per unit length, P, normal to the medial surface and viscous force per unit length, D_±, tangential to it) to the internal beam forces (F and G; see Fig. 6), the beam bending moment M, and the mass × acceleration of the slice itself. The resulting equations are as follows (see Hill, 1998):

The terms in these equations can, when multiplied by the slice length δs, be explained as follows, with reference to Fig. 6. In equation 9 (transverse forces), the first term represents the difference in the beam shear force F at the two ends of the slice, the second represents the force exerted towards the centre of curvature by the longitudinal tension G, the first term on the right-hand side is (transverse) body inertia, and P is the hydrodynamic pressure force. In equation 10 (longitudinal forces), the first term is the contribution from the shear force because the two end faces of the slice are not parallel, the second represents the difference in the longitudinal tension G at the two ends, the first term on the right-hand side is again body inertia, and D_± are the drag forces. In equation 11 (torque balance), the first term comes from the difference in bending moment, M, at the two ends of the slice, F represents the torque exerted by the shear force, and the third term represents the torque exerted by the difference in viscous drag on the two sides of the slice; the term on the right-hand side represents the rate of change of angular momentum of the slice. This term is new, in that it has not previously been included in analyses of fish, although it is well-known in the theory of dynamic bending beams (Wempner, 1973). It turns out to be quite small, even for large-amplitude motions. Also the (D₊− D₋ )l_z term is normally small for a thin fish with similar boundary layers on the two sides, and we henceforth neglect it.

Knowing the distribution of bending moment, M, along the fish and as it varies with time, it is but a small step to computing the bending moment M_m(s,t) actually generated by the swimming muscles themselves, as long as plausible values can be assigned to the passive elastic and viscoelastic properties of the body tissues (although very few data are available on such quantities). This has been done by Cheng et al. (1998), in an analysis of the small-amplitude steady-state swimming of saithe, for which morphological and kinematic data (i.e. measurements of m_b, I and h) were provided by Hess and Videler (1984). Fig. 7 shows the computed distribution along the fish of the amplitude and phase of the muscle bending moment, M_m, as well as of individual contributions to it. Two remarkable things can be seen. (i) The maximum amplitude of M_m occurs near the middle of the fish, where it is fattest, while there are significantly larger individual contributions close to the tail from both hydrodynamics and body inertia. These clearly compensate each other almost exactly. (ii) The wave of M_m travels down the fish much more rapidly than the wave of displacement; indeed, it is predicted to be almost a standing wave. This is qualitatively consistent with electromyographic measurements of the activation signal in saithe and other fishes (Wardle and Videler, 1993) and, more recently, in two species of tuna (Knower et al., 1999). Preliminary application of the large-amplitude theory (but in two dimensions and neglecting the drag terms) has been made by Hill (1998), with broadly similar results.

The above represents one approach to the investigation of how fish achieve their remarkable swimming performance: measure the kinematics (displacement wave, fin motions, etc); compute the hydrodynamics; and use the results to estimate the distribution of bending moment actually supplied by the locomotory muscles. Thus, one will be able to know (in principle) both the load experienced by the muscles and their rate of shortening and lengthening. Electromyography can provide data on the propagating wave of activation, so all the important factors in the muscle mechanics will be available for comparison with the results of experiments using isolated muscle preparations (Altringham et al., 1993); for a fuller discussion, the reader is referred to Wardle et al. (1995). The inevitable discrepancies will doubtless lead to many years of fruitful research.

Another approach is to work forwards from data and models of the central signal generator to calculate the activation wave. That can then be combined with information about the muscle properties obtained in the laboratory and with models of hydrodynamics and the passive internal properties to predict the displacement wave, h(s,t). This is the approach of Bowtell and Williams (1991, 1994) working on the lamprey. It is philosophically more satisfying, since it is what, in a sense, the fish actually does, but the calculation would have to be done iteratively, so the computational cost would become immense.

Finally, it should be remarked that steady-state undulatory swimming is a mode of locomotion that is employed only by some fish and only some of the time. The spectrum of behaviour is enormously wide: C-starts, S-turns, carangiform and lunate-tail propulsion (Lighthill, 1970; Weihs, 1972, 1973; Chopra, 1974; Wolfgang et al., 1999), balistiform and gymnotiform locomotion (Lighthill and Blake, 1990; Lighthill, 1990), propulsion using pectoral fins (Webb, 1973) and many more. The mechanics and control systems that lie behind them must be equally varied, and hydrodynamic theory, although richly informative, can provide only a small part of the desired understanding.

All the new models and results described in this chapter come from the PhD thesis of S.J.H. at the University of Leeds (Hill, 1998). We are grateful to the Wellcome Trust for supporting him and to the Engineering and Physical Science Research Council for the award to T.J.P. of a Senior Fellowship. We are also most grateful to Jianyu Cheng, John Altringham and (of course) Neill Alexander for the insight they helped us acquire. Finally, T.J.P. would like to acknowledge his profound debt to the late Sir James Lighthill who introduced him to biological fluid dynamics many years ago and whose early work on fish swimming was the inspiration for all that came since.