Self-propelled anguilliform swimming: simultaneous solution of the two-dimensional navier–stokes equations and Newton’s laws of motion

It is well known that an anguilliform swimmer propels itself forward by propagating waves of curvature backwards along its body. As a consequence of the interaction between the body and the surrounding water, forces due to fluid pressure and to drag act on the body of the creature. At any instant, an integration of these forces around the body will, in general, give rise to either an acceleration or a deceleration of the creature, not only in the forward direction but also laterally and in the rotational sense. After accelerating from a resting position, in the absence of external forces, the creature will tend towards an asymptotic state. In this state, the mean velocity in the forward direction, over a cycle of body movement, will be constant whereas the mean lateral and mean rotational velocities will be zero. Under these circumstances, the mean force due to pressure will be equal to the mean drag. Underlying this, however, a cyclic variation in all variables, whether forces or velocities in forward, lateral or rotational directions, will remain an ever-present feature of the creature’s behaviour. The aim of the present work is to investigate the behaviour of this system by means of a two-dimensional mathematical model of the fluid mechanics and associated dynamics of a time-dependent moving boundary responding to forces generated by its own movement.

Most previous studies of the hydrodynamics of aquatic locomotion have relied upon the assumption that water is an inviscid fluid. The pioneering work of Lighthill (e.g. Lighthill, 1960, 1969, 1970), the studies by Wu (e.g. Wu, 1971) and the work of Weihs (1972), Webb (1975), Hess and Videler (1984) and Cheng et al. (1991) are testimony to this approach. More recently, others have considered it necessary to quantify the viscous effects in the water in order to specify more accurately the drag forces acting on the body of the fish. To do this requires the techniques of computational fluid dynamics, and studies by Carling et al. (1994), Williams et al. (1995) and Liu et al. (1996, 1997) are now beginning to show how these techniques can be used.

In addition to the inviscid assumption, most studies make the further assumption that some point on the creature’s body, say the nose, remains stationary whilst a uniform upstream flow is prescribed. This could represent the behaviour of a fish in a flow tank or it might approximate swimming at constant speed in an open domain. However, this assumption does not allow the velocity to change in response to the cyclic variation in forward force that emerges from the hydrodynamical calculation; for example, Hess and Videler (1984) state that ‘the calculated thrust fluctuated during each half-period between zero and approximately twice its average value’. To be able to understand the biological significance of the fluid dynamics it is necessary to take into account the cyclic variation in both force and velocity. Even with viscous effects taken into account, the imposition of a specified forward speed may give rise to an asymptotic state in which the mean forward force due to pressure, as computed by the model, cannot be expected to balance the mean drag. Thus, in previous work in two dimensions, both Williams et al. (1995) and Liu et al. (1996) have found the calculated mean thrust to be greater than the mean drag. The three-dimensional study by Liu et al. (1997) also shows this effect, although the magnitude of the net force is smaller than with their two-dimensional work. The appearance of net thrust in all these studies indicates that the results presented, such as the creature’s position, the flow field and the force distribution, could only come about by restraining the creature’s movement by extraneous force. In other words, those model creatures were not swimming freely but were effectively tethered and being prevented from accelerating away from their specified positions. Any forces calculated in these studies could differ markedly from the forces developed by a freely swimming animal.

In the present study, the model creature is not subject to these constraints and simply finds its own path according to the forces acting upon it. It is self-propelled.

The model creature described here has a predetermined time-dependent body shape in the form of a backward-travelling wave with increasing amplitude from head to tail. Both the position and movement of this body are calculated from the forces produced by its interaction with the water. It is the computational equivalent of a clockwork fish in which the mechanically driven changes in the body shape propel the creature through the water.

The computational model may be considered partly as a calculation to determine the motion of the creature’s centre of mass and partly as a computation of the fluid dynamics. This is related to the recoil correction in the work of Lighthill (1960, 1970), although the present approach extends that basic idea. In practice, the complete model is a single computation iterating to find a balance between an acceleration required at any instant to achieve a new position and the force on the body that then arises from the surrounding water in response to the new shape and position. In the following description of this system, the determination of the changing shape and position of the creature’s body will be considered first and then the computation of the hydrodynamics. These two aspects of the creature’s behaviour will then be brought together by considering the dynamical interaction between the body and the surrounding water.

For ease of computation, all the variables and constants in the system have been converted to dimensionless quantities, by dividing by characteristic values. All distances have been divided by a characteristic length, d₀, which has been defined as the amplitude of the travelling wave at the tail. Time has been divided by a characteristic time, t₀, which is the locomotor cycle duration. (A list of symbols and their definitions is provided at the end of Appendix B.)

A distinction is made between the shape of the creature’s body and its position and orientation relative to the enclosure in which swimming takes place. Two coordinate systems are used, as shown in Figs 1 and 2. One of these (x_s,y_s) describes the body shape, independently of the enclosure. The second coordinate system (x,y), fixed in the enclosure, enables the position of the creature’s centre of mass to be defined as (x_m,y_m). The final link between these two systems is the angle θ_m, shown in Fig. 1B, which defines the orientation of the x_s axis relative to the x direction in the fixed coordinates. All three quantities, x_m, y_m and θ_m, vary with time. The central feature in the present work is that, although the creature’s shape (x_s,y_s) is a prescribed function of time, its position is not predetermined. Instead, the creature’s position (x_m,y_m) and its orientation, θ_m, arise as part of the computation of the dynamical system.

A numerical integration of equation 2 using an increment in s of Δs can then be used to determine x_s, thus completing the specification of the body centre-line coordinates.

Equation 1 gives an idealised description of the changing body shape of an anguilliform swimmer (see also Hess, 1983). The ordinate y_s(s,t) represents a backward-travelling wave with amplitude increasing linearly from b/(1+b) at the nose to 1 at the tail, as s increases from zero to l. This gives a wave that travels down the body once in each locomotor cycle; in other words, the wavelength is 1 body length. The parameter b is included to allow an amplitude other than zero at the nose of the creature. The value used in the solution illustrated here is b=0.25, giving a maximum amplitude at the nose of 0.2 compared with the amplitude at the tail. The length of the creature is 8 cm (representing a small eel, as in Gray, 1933), and the amplitude of the travelling wave at the tail is 1 cm. Dimensionless length, l, is therefore 8.

Equation 3 gives a body outline which tapers from width w₀ at the nose of the creature to w₁ at the tail and is parallel to the body centre-line at the nose and the tail. The values used in the present study are w₀=0.64 and w₁=0.16, giving a creature with a body width at the nose equal to 8 % of body length tapering to a tail width of 2 %.

Equations 1 and 3 characterise the body of the creature, with equation 1 defining the centre-line and w_s in equation 3 expressing body width perpendicular to this centre-line. Starting with w₀ at the nose, the outline shape has been constructed by forming quadrilateral segments such that, for each increment Δs along the centre-line, the area of each segment is constrained to be precisely w_sΔs, where w_s is the width at the mid-segment position. In addition, within each segment, the areas on either side of the centre-line are equal. The body shape and the body segments shown in Fig. 1B correspond to the parameter values given above. The segment area constraint conserves the total area of the creature’s body as its shape changes with time, as required for the satisfactory computation of two-dimensional, incompressible fluid flow. Effectively, this conserves the creature’s mass for a two-dimensional system in which the aggregate density of the body tissues is assumed to be constant.

Having defined the creature’s shape, the position of the centre of mass in the local (x_s,y_s) coordinate system is computed. The local coordinates defining shape are then related to the global coordinates, which will eventually determine the creature’s position in its enclosure. The domain of computation has been chosen to be a rectangular tank 25 body lengths by 4 body lengths, as shown in Fig. 1A. The global coordinates (x,y) are measured from the lower left-hand corner of the tank, as viewed from above, giving overall dimensions of x=200 and y=32, in dimensionless units. An impression of the position of the body of the fish can be seen in Fig. 1A in the centre of the tank at a point part-way through the overall journey of the creature. To arrive at this position, the animal starts from rest at a point approximately 4 body lengths from the right-hand end of the tank, in the mid-line. The travelling wave then begins to develop (described in detail below) and, in response, the creature propels itself forward, close to the mid-line, until it reaches the position shown.

The coordinates of the centre of mass (x_m,y_m) are shown in Fig. 1B. The remaining variable in Fig. 1B, the angle θ_m, enables the relationship to be established between the global coordinates (x,y) and the local coordinates (x_s,y_s) which define the body shape as given in equation 1. The line at angle θ_m shown in Fig. 1B passes through the centre of mass. This line is parallel to the axis of generation of the body shape, that is, the x_s axis. The x_s axis itself, as indicated in Fig. 2, does not necessarily pass through the centre of mass.

The motion of the creature’s body has now been specified in terms of three unknown quantities x_m, y_m and θ_m, all three of which are functions of time and must emerge as part of the computation of the dynamical system. Knowing x_m, y_m and θ_m at any instant and given the prescribed shape change enables a complete determination of velocity components around the outline of the creature’s body. This velocity distribution is required in order to provide boundary conditions for the computation of the fluid dynamics.

where μ is the dynamic viscosity of water. For the common anguilliform fish, the Reynolds number based on the above definition lies in the range 50–5000. Although the Reynolds number is high, viscous effects cannot be ignored because in all cases there will be regions of relative stagnation. For example, these will occur at vortex centres, between vortices rotating in the same sense and in the far field. In addition, the viscous driven flow around the body is of paramount importance in determining the drag forces as the creature moves through the water.

The boundary condition for velocity at the body of the creature is the so-called no-slip condition in which the fluid velocity components (u, v) at the boundary are set to the velocity components of the body itself and are consistent with the current shape and position of the body as described above. The boundary conditions remote from the body assume the creature to be immersed in a fixed tank and so, at the boundaries of the tank, u=v=0.

Equations 4 and 5 must now be solved numerically, together with the continuity constraint (equation 8), to yield a velocity field (u, v) and a pressure field, p. On a rectangular grid, this solution process is relatively straightforward. Here, whatever grid is used for numerical solution, there must be some means for tracking the changing shape and position of the creature as it moves through the tank. There are advantages in choosing a grid which adapts to follow the motion of the creature. This enables pressure and velocity fields around the body to be readily available for use in the calculation of the body dynamics. Such a grid must be non-rectangular and must change with time as the creature moves through its cycle and makes progress on its journey through the tank. An example of such a grid is shown in Fig. 3. On a non-rectangular grid, the numerical solution becomes more involved and there is essentially a choice between the finite element method (see Zienkiewicz, 1977), the finite volume method (Patankar and Spalding, 1972) and a finite difference approach used in conjunction with an appropriate coordinate transformation. Here, transformation methods are preferred; such methods work well when applied to rapidly moving boundary problems, such as the present swimming problem.

The governing equations 4, 5 and 8 are expressed in terms of the global Cartesian coordinates (x,y). These equations must be transformed to new coordinates ⁠, where these new independent variables follow the grid. The time variable, t, remains unchanged under the transformation and so ⁠. Accordingly, a set of local transformations is established in the given moving curvilinear mesh, such as that shown in Fig. 3. The details of the transformation are given in Appendix A. The grid is characterised by the angles and ⁠, as shown in Fig. 4A, where and are defined in Appendix A. As can be seen in Figs 3 and 4, the grid used here has been chosen to be orthogonal. The reasons for making this choice are explained in detail in Appendix A. Briefly, such a choice gives computational advantages by allowing the derivatives of and to represent curvature in the computational domain. In addition, new velocity components (U, V) can be defined to be locally perpendicular to one another and to lie in the local directions, respectively. The computational variables corresponding to U, V and pressure p are defined at discrete points, as shown in Fig. 4B, and so for cell i,j the computational counterparts to the components of velocity are U_i,j and V_i,j, respectively, and the corresponding pressure is p_i,j. Arranging the computational variables in this way simplifies the computation. The procedure for satisfying the continuity constraint, together with calculating the pressure, can then be established without computational complications.

where the transformed variables ⁠, the operators and the quantities S₁, S₂, S₄ and S₅ are defined in Appendix A (see equations A57 and A58). A brief explanation of the numerical solution of equations 10 and 11 is also given in Appendix A. This is accompanied by results for a test case in which, for a representative geometry, comparison can be made with an analytical solution. For a grid comparable in resolution to the coarsest part of the full grid shown in Fig. 3, the degree of agreement between the numerical and analytical solutions is good.

Each segment of the body is subjected to forces due to fluid pressure and to drag. Fig. 5 shows these forces for a segment i along the body. The position of segment i is defined by polar coordinates (r(i),ϕ(i)), where these are expressed relative to the centre of mass and the current inclination of the x_s axis. This axis is shown in Fig. 2 and its inclination, θ_m, is indicated in Fig. 1B. So, in Fig. 5A it can be seen that the angle α(i) is defined as:

In equation 12, θ_m is a function of time only, whereas ϕ(i) is a function of both time and position along the body. Furthermore, θ_m is an unknown function whereas ϕ(i) is a known function, calculated from the given body shape (equation 1).

Fig. 5B shows the total force and moment acting on the centre of mass and also the local forces acting on a segment. The local force f_p(i) due to pressure acts normal to the segment–fluid interface in the direction of of the local coordinate transformation. The drag force f_d(i) acts tangentially, in the direction of x̄. Local forces due to pressure and drag are taken into account at the nose and the tail of the creature as well as along both sides. Both local pressure and local drag forces are readily available from the computation of the fluid mechanics. Their determination is straightforward given the dependent variables p and U associated with the fluid cell adjacent to each segment interface. For example, knowing the velocity U in the x̄ direction along the interface allows the velocity gradient normal to the interface, ⁠, to be determined. This velocity gradient is then used to calculate the local drag. The grid resolution around the body in the transformed domain has been chosen to ensure that the calculation of local pressure force and drag is sufficiently accurate. The degree of accuracy which might be expected in a typical part of the full computation is demonstrated in the test case in Appendix A.

Having found the local forces acting on each body segment, these may be resolved in the global x and y directions and then summed over all body segments to give the total force. This may be regarded as acting on the centre of mass. These components are shown as F_f(t) and F_l(t) in Fig. 5B, expressing, respectively, the instantaneous force in the forward direction propelling the creature along and the lateral force that arises as a consequence of the cyclic movement of the body in the water. To complete the specification of the forces about the centre of mass, the total moment, M(t), of all the local forces must also be calculated. Again, this is straightforward given f_p(i) and f_d(i) and knowing the geometry.

where Δt is the time step size; it should be noted that a(i) does not change with time, only with segment position, as discussed above.

The right-hand sides of equations 19 and 20 are equal to the rate of change in the position of the centre of mass. Hence, the forward velocity can also be written as ⁠, and the lateral velocity, v_l, as ⁠.

Equations 22 and 23 are seen now to govern the dynamics of the creature’s body in the forward and lateral directions expressed as time-dependent variables associated solely with the movement of the creature’s centre of mass. During the iterative process used to solve these equations, taking equation 22 as an example, a small change in F_f(t) is fed back to the kinematics by means of a single, corresponding change in the velocity v_f(t). This is explained in more detail in Appendix B.

During this part of the computation, it is important to keep the newly established positions of the body segments, i.e. the positions at time t and the new position of the centre of mass, constrained such that equations 19 and 20 are satisfied precisely throughout the iteration. In order to do this, the new angular position of the x_s axis, θ_m(t), is required. Analysis of the rotation must therefore be taken further.

Equation 28 is now the rotational counterpart to equations 22 and 23, thus completing the specification of the dynamics of the body.

The physical basis for the interaction between the body and the surrounding water has been set out above, in essence in equations 22, 23 and 28. Within each time step, an iteration has to take place in order to satisfy these equations. Controlling this iteration is the key to obtaining satisfactory computational results. The numerical procedure must be stable and has to produce results of an acceptable accuracy in a reasonable number of iterative steps. The main obstacle to achieving this is that a small departure from the correct position of the creature’s body, as may well occur during any iterative process, produces a relatively large corrective force arising from the fluid dynamics. Without due care, this process can easily become numerically unstable.

Appendix B gives details of the numerical procedure used here to update the position of the creature’s centre of mass. In essence, each velocity component, v_f(t), v_l(t) and v_r(t), needs to be corrected cautiously at each iterative step. The first approximation to the velocities then determines a new position of the centre of mass and a new orientation of the creature’s body. At this position, new forces are computed arising from the fluid dynamics. These forces produce new accelerations which in turn give rise to adjusted velocities. The formation of new velocity approximations, each by a small incremental change from the previous value, then allows the next iteration to begin. This continues until all three components of force, F_f(t), F_l(t) and M(t), have converged. The forces must be used to test convergence because it is the forces, related to the second derivative of positions, that are varying most rapidly during the iteration. The mathematical detail of the iterative process is described in Appendix B.

In the numerical calculation, it is necessary to specify the initial conditions, the simplest being a creature at rest in still water. From this state, some form of starting manoeuvre must be used to develop the travelling wave on the body. Two possibilities come to mind, which resemble the commonly observed forms of behaviour during fast starts (Webb, 1976, 1978). The first is to increase the speed of the travelling wave gradually from zero to its final value. This produces a start from an initial ‘S-shape’. The second is to scale the amplitude of the body shape y_s, again increasing gradually from zero to the distribution defined in equation 1. This second manoeuvre is the one used here. The combined effect of the increasing amplitude and the travelling wave produces an apparent ‘C-shape’.

Mathematically, the important feature of the given variation in β is that it ensures a continuity in both body shape and the rate of change of shape as the creature settles into its rhythmic behaviour at t=1.

Using the computational model described above, results are presented here for a creature 8 cm long, swimming with a cycle time t₀ of 1.2 s. The travelling wave has an amplitude at the tail d₀ of 1 cm and a wavelength of 1 body length. These values approximate the young eel in the work of Gray (1933). Other geometric parameters defining the body shape are as described above. The Reynolds number based on d₀ as characteristic length and d₀/t₀ as characteristic speed is Re=83.

At the start of the computation, the creature is positioned with its nose at the point x=155.5, y=16.0 and lies with the centre-line of its body along the mid-line of the tank (see Fig. 1A for tank dimensions). The body curvature and backward-travelling wave are then allowed to build during the first cycle of activity as described above. At the end of the first cycle, the creature’s body shape and its change with time are fully established. This established pattern of changing body shape then continues for a further 19 cycles of activity, as the creature accelerates and swims through the tank. In the computation, each cycle is subdivided into 32 time steps, to provide sufficient resolution for the rapidly changing forces which act on the body.

The results are presented as follows. Figs 6–9 show the structure of the fluid flow at various stages of the creature’s journey. The kinematics and dynamics of the creature’s centre of mass are shown in Figs 10–13. Fig. 14 shows the amplitude envelope and the nose and tail trajectories that result from the calculations.

Streamlines (contours of the instantaneous stream function) can be a useful way of representing a two-dimensional velocity field. In the present study, the stream function distribution has been recovered from the velocity field, after the computation has been completed, by a process of integration. Fig. 6A shows the streamline pattern at the end of cycle 5 of the creature’s journey. The general direction of the flow ahead and to the side of the nose of the creature is indicated by the arrows. Red streamlines represent flow circulating in a clockwise direction, and green streamlines represent anticlockwise flow. At the end of cycle 5, the nose of the creature has reached a position at which x=138.0 representing a total journey of 2.19 body lengths from the starting point. Fig. 6 shows the progress of the creature at intervals of three cycles. In Fig. 6B, at the end of cycle 8, the nose is at x=121.3, representing a further overall progress of 2.09 body lengths during the three cycles from cycle 5 to cycle 8. In moving from cycle 8 to cycle 11 and then to cycle 14 (Fig. 6B–D), the nose of the creature makes progress in two approximately equal steps of 2.16 body lengths. This suggests that after an initial acceleration over the first eight cycles an asymptotic state is being approached. The mean velocity of the creature in this state, based on the movement of the nose, is 0.60 body lengths s⁻¹.

In each panel in Fig. 6, the leading or outer streamlines correspond to a stream function value of +1.5 for anticlockwise flow and −1.5 for clockwise. Since the stream function is zero along the side walls of the tank, this means that the volume flow of water returning past the creature between each outer streamline and the corresponding side wall is the same in all the plots. From the outer streamline working inwards, the increment in stream function from contour to contour is 1.0 in all plots, indicating equal increments in volume flow between adjacent streamlines. It follows, therefore, that the closer the streamlines are to one another the faster will be the local fluid velocity. Interpreting Fig. 6 using this information enables an overall impression to be gained of the two-dimensional water flow in an enclosed tank through which an anguilliform creature is accelerating.

Regions of self-contained circulation can be seen in Fig. 6. To depict these regions more clearly, Fig. 7 shows the ‘separating’ streamlines for the contour plots at cycles 8, 11 and 14. These are superimposed for comparison and correspond to Fig. 6B,C,D, respectively. Separating streamlines are drawn by first identifying the stream function values corresponding to saddle points on the stream function surface. These saddle points represent instantaneous stagnation points in the fluid flow, which occur between vortex structures rotating in the same direction. This gives rise to flows locally which are in opposite directions. Viscous interaction between these flows then gives rise to a point between the flows at which the fluid velocity vector is zero. These stagnation points appear in Fig. 7 as small open circles, three in each of the three plots. The cyclic repetition of the flow structure in the neighbourhood of the body is apparent (compare b, c and d in Fig. 7). Although the physical positions of these separate structures are similar, the strength of circulation within corresponding parts need not necessarily be the same.

Figs 8 and 9 show the progress of the creature through cycle 15. Fig. 8 shows the full streamline pattern, continuing from Fig. 6D, with Fig. 8B showing streamlines and body shape half-way through the cycle. The flow structure can be seen to be an approximate reflection in the mid-line of that in Fig. 8A. Fig. 8C then depicts the completion of the body shape cycle for cycle 15. Fig. 9 shows the cyclic behaviour in more detail, on an expanded scale and including more time steps during cycle 15. The stream function increment used to produce the plots in Fig. 9 has been reduced to 0.5, half that in Fig. 8, in order to highlight the local flow. The major circulating flow structures developing around the body can be seen. Following their formation, these structures are first displaced backwards slightly, dictated by the backward-travelling wave, and then tend to remain stationary for a short time as the creature travels forwards. Shortly after this, the effect of viscosity tends to break down the rotating structure and the fluid becomes part of the larger-scale circulation behind the body (see Figs 8A,B, 9A–E). In general, the panels in Fig. 9 show three circulating fluid structures near the body, which were depicted explicitly in Fig. 7. The first circulation, on one side of the body near the nose (on the right side, e.g. in Fig. 9A), is seen to grow and gain strength, enhanced by the viscous boundary layer along the side of the body. There is a second circulation on the opposite side of the body at approximately the mid-length. A third circulation, near the tail, is seen on the same side of the body as the first. The growth and dissolution of these structures in cyclic repetition is the key to understanding the behaviour of the flow. The growth is caused by the motion of the body and the viscous effects between the body and the water. The circulation within these structures is maintained as the creature passes by. Finally, the local circulation is diminished by viscous dissipation as the creature moves away.

In all the stream function plots, but most clearly in Fig. 9, a region behind the body can be seen in which no streamlines are plotted. The demarcation of this region can be seen to be the outer streamlines. The volume flow between these streamlines is therefore 3.0 (since the stream function values are ±1.5). From the physical distance separating these streamlines, it is therefore possible to estimate the mean velocity in the wake following behind the creature. Over cycle 15, this wake velocity is in the range 0.25–0.5 body lengths s⁻¹ in the direction of forward motion. Since the mean forward speed of the creature is 0.60 body lengths s⁻¹, the relative velocity of the wake, as seen by the creature, is in the range −0.35 to −0.1 body lengths s⁻¹, depending on the point in the cycle. Thus, a recession of the wake is apparent to the creature. This relative speed would be expected to change, however, depending on the dimensions of the tank, particularly the distance of the side walls from the creature.

Fig. 10A shows the trajectory of the centre of mass of the creature during cycle 15; the forward and lateral components of the progression are shown in Fig. 10B. A cyclic lateral movement can be clearly seen as the creature travels from point a to point i. There is also a small cyclic change in the rate of forward progression, which appears as a barely discernible variation in the gradient of the upper curve of Fig. 10B. The second derivative of each curve in Fig. 10B, which gives the respective component of acceleration, varies far more noticeably over a cycle of body movement. This can be seen from the force (acceleration) vectors in Fig. 9.

Fig. 11A shows the velocity components associated with the movement of the centre of mass in the forward and lateral directions over 20 cycles of body movement. The full-cycle running mean velocity is also shown in each case. For each cycle of body movement, two cycles of oscillation in forward velocity occur, one for each leading edge of the backward-travelling wave. From the mean value of the forward velocity trace, an initial phase of rapid acceleration can be seen. During approximately one cycle of activity from a point approximately half-way through the first cycle, the mean forward velocity increases from near zero to half its asymptotic value. After a further three cycles of body movement, the mean velocity has attained 90 % of its final value. Finally, at around cycle 9, the asymptotic state appears to have been reached. Here, the amplitude of the oscillations in forward velocity represents a departure of ±9 % from the mean. The creature then continues, largely undisturbed, until approximately cycle 15. From this time onwards, a slight deceleration is apparent, and by cycle 20 a reduction in the mean velocity of approximately 3 % has occurred.

For the oscillations in forward velocity after cycle 15, it can be seen that the two peaks per cycle become unequal. This behaviour arises because, although the specified changing body shape exhibits left–right symmetry, the creature is not swimming precisely along the mid-line of the tank (see Fig. 10A). This is because the starting manoeuvre produces an initial trajectory of the centre of mass in the form of a spiral, ending with a sideways displacement as shown in Fig. 10A. All starting manoeuvres are necessarily asymmetric and so some displacement of the asymptotic trajectory is to be expected. The surrounding fluid domain is therefore not symmetrical and this gives rise to asymmetric forces. The mathematical creature, as defined here, has no mechanism for countering this asymmetry by corrective action. Furthermore, it is probable that the tendency away from symmetry becomes more pronounced as the creature approaches the end wall of the enclosure. At the end of cycle 20, the position of the creature’s centre of mass is approximately 7 body lengths from the end wall, well within the domain of influence in which such effects might be expected.

The asymptotic mean velocity of the centre of mass, indicated in Fig. 11A, is 0.60 body lengths s⁻¹. This forward speed corresponds to a stride length of 0.72 body lengths cycle⁻¹. The velocity of the backward-travelling wave is also indicated in Fig. 11A. The ratio of asymptotic mean forward velocity to backward-travelling wave velocity is 0.77 for this case.

The lateral velocity trace is also shown in Fig. 11A. In this case, only one cycle of oscillation in velocity occurs for each cycle of body movement. The amplitude of the oscillations in lateral velocity is ±10 % of the asymptotic mean forward velocity. The running mean lateral velocity is close to zero, but swings from slightly negative to slightly positive values during the overall journey.

The propulsive forward force and the drag on the body of the animal, which give rise to its velocity, are shown in Fig. 11B. After 1.5 cycles from the start, the drag reaches a maximum and thereafter tends to decrease as the animal and the fluid flow develop a cyclic pattern. The dashed line indicates the full-cycle running mean in the forward pressure force. It can be seen that when the mean velocity is at or near its asymptotic value the mean pressure force tends to balance the drag over a cycle of activity. Fig. 11B shows that after the creature has settled into its asymptotic state the maximum excursion of the forward pressure force is approximately eight times the mean drag. The acceleration and deceleration corresponding to these positive and negative net forces are reflected in the corresponding increase and decrease in forward velocity in Fig. 11A. A change in the pattern of the forward pressure force can be seen clearly from cycle 17 onwards. The two peaks in pressure force that occur per cycle become unequal, corresponding to the behaviour in the forward velocity. As argued above, such asymmetric forces are to be expected when the creature is swimming in a domain without precise lateral symmetry.

The lateral forces are shown in Fig. 11C. Overall, there is one cycle of oscillation in the lateral pressure force per cycle of body movement, although near to each maximum a small secondary peak occurs. A similar behaviour can be seen around the minimum. The maximum excursion in the lateral pressure force is approximately 0.55 times that for the forward pressure force (Fig. 11B).

Fig. 12A shows the angular position of the body axis, θ_m, plotted against time. In the asymptotic state, the amplitude of the oscillations in θ_m is approximately ±0.1 rad. The effect of the cyclic variation in θ_m can be seen in Fig. 2, which shows how the orientation of the x_s axis varies during cycle 15. The total angular velocity, v_r, is shown in Fig. 12B, where v_r from equation 27 includes both a contribution due to the rate of change in θ_m and also that due to changing body shape. This total angular velocity arises from the moment about the centre of mass (Fig. 12C).

Examining more closely the variation of the forces in the asymptotic state, Fig. 13A,B shows the forward and lateral components of the net force during cycle 15. Points a–i in Fig. 13A indicate the forces acting on the creature at times in the cycle corresponding to the body positions shown in Fig. 9A–I, respectively. The force vectors in Fig. 9 correspond to the components of force shown in Fig. 13A,B. During the second half of the cycle, the maximum net propulsive force occurs when the creature is approximately half-way between the positions in Fig. 9G and Fig. 9H. There is a significant lateral force acting at the same time, however, producing an aggregate acceleration of the centre of mass not only in the forward direction but also to the creature’s left. At around this point in the cycle, the creature is farthest to the right in its trajectory (see Fig. 10A). The lateral force is therefore seen to be restoring the overall motion of the creature towards the mean direction of travel.

The creature’s position is also determined partly by the moment of forces about the centre of mass. Fig. 13C shows the variation of this moment during cycle 15. It can be seen that between points g and h the moment is positive, or anticlockwise, and from Fig. 2 the rotation of the x_s axis is also anticlockwise. The combined effect of this rotation and the lateral progression of the centre of mass is to produce a relatively large lateral movement of the creature’s nose, as can be seen in Fig. 2.

The rotational effect occurs along the whole body, as illustrated in Fig. 14, which shows the body centre-line amplitude envelope and the trajectories of the nose and the tail. This shows both the prescribed linear envelope of the backward-travelling wave and the resulting non-linear envelope found in the computed behaviour. The difference between the two envelopes is caused by the rotation of the body about the centre of mass in combination with its lateral movement.

The trajectory of the tail exhibits the expected figure-of-eight pattern, with the solid triangle indicating the direction of travel around the trajectory. The nose trajectory also exhibits a figure-of-eight form. The direction of travel indicated for the nose (open triangle) is shown at the same point in the cycle as for the tail. Thus, it can be seen that, for a small part of the cycle, the lateral direction of travel of the tail is opposite to that of the nose. This is brought about by the effect of rotation, since for the wavelength of 1.0 used here, the trajectories of the nose and tail would otherwise cycle in phase.

The results presented here illustrate the behaviour of a dynamical system representing the motion of an anguilliform swimmer. The present model differs from previous computational models in that the forward swimming speed of the creature has not been specified beforehand but has arisen from the computation. The mean forward speed attained is that at which the forward pressure force balances the backward drag, over each cycle of activity. In contrast, in previous studies in which viscosity was included in the computation (Williams et al. 1995; Liu et al. 1996, 1997), the calculated forward force due to pressure was greater than the calculated backward drag. Hence, those creatures were effectively prevented from swimming forward.

One consequence of this artificial imbalance may be in the resulting vortex structure in the wake. In the previous two-dimensional studies, vortices shed by the creature remained as clearly identified structures in the wake (see Fig. 8 in Williams et al. 1995; Fig. 11 in Liu et al. 1996), whereas in the present study the vortices join the overall circulation soon after being shed from the creature (Figs 6, 9). As those creatures were effectively restrained, much of their body motion was probably used to build up the strength of the generated vortices. Those circumstances perhaps resemble a swimmer anchored in a stream while at the same time engaging in vigorous paddling of the feet.

Another assumption of previous work has been that the forward speed of the creature is constant. At first, from observing either real or calculated swimming, this may seem reasonable. For example, the forward progression shown in Fig. 10B does not appear to differ significantly from a straight line. However, the time derivative of this progression, the forward velocity, does undergo significant oscillation, as can be seen in Fig. 11A. This oscillation in velocity reflects the fact that the net forward force on the body fluctuates between positive and negative values, and the excursions are even larger than the oscillations in lateral force. Since it is the forces rather than the velocities that are important to the understanding of the biomechanics of swimming, oscillations in all forces must be allowed for in the mathematical formulation. Models that assume a constant forward speed cannot be used to reach reliable conclusions about the development of forces during swimming.

The importance of rotational forces in the behaviour of the model creature is illustrated by the fact that a linearly changing amplitude with respect to the axis of generation, x_s, produces a non-linear profile with respect to the axis of forward progression (see Fig. 14). In addition, the imposed non-zero amplitude at the nose has been enhanced considerably by rotation of the body about its centre of mass. Published accounts of eel kinematics do show a significant non-zero movement of the nose during swimming (Grillner and Kashin, 1976; Hess, 1983). However, there is no obvious way to disaggregate lateral and rotational movements in kinematic data, since there is no straightforward way to discover a unique equivalent to the axis of generation, x_s. Hence, choosing the value for nose amplitude in the model is not straightforward. The choice of wavelength raises further complications. Published values of wavelength for anguilliform swimmers (Gray, 1933; Grillner and Kashin, 1976; Hess, 1983) are in the range 0.79–0.67 body lengths. This suggests that between 1.2 and 1.5 waves are on the body at any time. However, near the head and the tail, eel body-shape data are difficult to interpret. It is not clear to what extent the lateral movement of the centre of mass, the creature’s rotation and the increasing wave amplitude contribute to the final amplitude profile along the body. Over the first 20 % of the body length, where the lateral movement is small, the observed amplitude envelope lies in the range from 0.1 to 0.2 times the amplitude at the tail. In addition, the body curvature is smaller in both the head and the tail regions. For these reasons, the published pictures of anguilliform swimmers (e.g. Gray, 1933) tend to show approximately one wave on the creature’s body at any time. Thus, the effective wavelength in anguilliform swimming may be larger than previously supposed. As a consequence of the interdependence of these geometric parameters, the following configuration was chosen for the present study: a wavelength of 1.0 body length and a linearly increasing amplitude envelope with a value at the nose of 0.2 times the value at the tail. In future work, the effects of using different nose amplitudes and wavelengths will be investigated. The results of such studies should produce an amplitude envelope that more closely matches the published profiles.

Forces and movements were originally calculated in the lateral and forward directions only, on the assumption that moments of force and rotation about the centre of mass would have only minor effects. However, the resulting forward pressure forces in that system were significantly greater than in the system described here, and the asymptotic mean forward velocity reached a value some 40 % greater than the asymptotic velocity found here. As Fig. 11A indicates, such a value would then be slightly greater (by some 8 %) than the velocity of the backward-travelling wave (J. Carling, T. L. Williams and G. Bowtell, unpublished data). This is not physically realistic. Imposing zero rotation on the creature was equivalent to providing an externally applied moment of force with the required direction and magnitude to prevent the animal rotating about its centre of mass. The achievement of a forward velocity greater than that of the backward-travelling wave indicates that the creature was able to gain net forward force from this implicit externally applied moment. Similar arguments apply to the imposition of zero acceleration in the forward direction, thus emphasising that it is ill-advised to assume a constant forward velocity in modelling swimming fish.

The work described here is based on two-dimensional fluid flow in the plane of the anguilliform motion. The results may therefore be interpreted as representing the behaviour of a swimming sheet extending in directions perpendicular to the plane of swimming. The results might also be thought to approximate the behaviour of an anguilliform swimmer in shallow water within an enclosure. Clearly, it is desirable to develop a three-dimensional version of the model, and work is currently in progress to achieve this. Until such work is completed, it is difficult to judge whether the effect of removing the two-dimensional restriction will substantially alter the results, particularly the magnitude of the forces acting on the creature. Nevertheless, whatever differences the three-dimensional results may produce, the present study is an important step since it represents the first dynamic model in which the position and the movement of the creature are not prescribed; it is a self-propelled swimming shape.

The idealised shape used to represent the body of an anguilliform swimmer exhibits a square-cornered nose and tail. This shape, not an obvious choice biologically, facilitates the particular solution procedure used here to solve the flow equations. It is fair to question whether the decapitation affects the general form of the flow and the forces generated. Close to the rectangular corners, it is clear that the local fluid velocity and pressure will depend on the local shape. However, as Fig. 9 shows, the main rotating fluid structures are created and sustained not only by displacement as the creature moves through the water but also by the action of the viscous layer along the side of the creature’s body. The general shape and magnitude of these vortex structures would probably not be affected to any substantial degree by changing the shape of the nose or the tail. The presence of sharp corners at the nose might enhance the strength of the circulation near the nose, but further along the body it is unlikely greatly to affect the overall behaviour. A greater influence is likely to be the effect of the lateral movement of the nose and its cyclic acceleration and deceleration. Taking proper account of the movement of the nose, as has been done here, may be the more important factor, rather than the precise detail of its shape. Nevertheless, to test this contention, it would be useful to experiment with a variety of nose shapes, and such a study is a longer-term aim.

The 8 cm anguilliform shape used here reaches an asymptotic mean forward speed of 0.72 body lengths cycle⁻¹. Although this is within the range obtained from experimental work on non-anguilliform swimmers, it is larger than that found by Gray (1933) for a 7 cm eel (Anguilla vulgaris) and a 13 cm butterfish (Centronotus gunnelus), which were swimming at mean forward speeds of approximately 0.4 and 0.5 body lengths cycle⁻¹, respectively. Hess (1983) calculated a value of 0.55 body lengths cycle⁻¹ for a 14 cm eel (species unidentified). Upon comparing the ratio of the forward velocity with the velocity of the backward-travelling wave, a similar pattern is apparent. Here, this ratio is 0.77, which is somewhat larger than the values of 0.65 and 0.67 found by Gray (1933) for the eel and butterfish, respectively. The value obtained by Hess (1983) was 0.70. Direct quantitative comparisons must be made with caution, however, since the present two-dimensional model may well give rise to higher swimming speeds than a more realistic three-dimensional equivalent.

The results have shown (see Fig. 11B) that the mean drag reaches a maximum during the early part of the creature’s journey and then reduces towards an asymptotic value as the cyclic behaviour of the flow becomes fully established. This comes about because, in the early stages, before the vortex structure around the body has developed its final cyclic pattern, there tends to be a larger difference between the speed of the creature and that of the slowly rotating fluid adjacent to the body. This difference gives rise to a steeper velocity gradient at the creature’s body and hence a greater drag force. As the speed of the creature increases, however, so too does the rotational speed within each vortex. The relative velocity between the creature and the surrounding water then tends to be less than it was in the early stages, and the smaller velocity gradient accounts for the resulting reduction in the mean drag. In this sense, this creature is gaining some small benefit from the inevitable presence of the vortex structures.

The present results clearly support the view that an anguilliform swimmer can use most of the full length of its body to generate vortex structures. In contrast, the generation of vortices in carangiform swimming is more closely associated with the movement of the tail and that part of the body near the tail. For example, Müller et al. (1997) have published results of experiments using particle image velocimetry to show the tail-dominated vortex wake behind a mullet (Chelon labrosus Risso). As Müller et al. (1997) point out, however, Rosen (1959) observed attached vortices on the body. Thus, just as there may not be a rigid demarcation into distinct modes of anguilliform and carangiform swimming (Breder, 1926; Webb, 1975), so there may be a spectrum of mechanisms of vortex generation, with both viscous and displacement mechanisms always present to some degree but with one or the other being predominant in each individual case.

The value of Reynolds number (Re) used here to non-dimensionalise the computation is based on d₀, the amplitude of the travelling wave as it reaches the tail, and d₀/t₀, where t₀ is the cycle duration. For the present case, this gives a value of 83. As far as the numerical work is concerned, the definition of Re does not affect the solution, since these same values of length and velocity are subsequently used to assign dimensions to the computed values. However, another common use of Reynolds number is to classify problems in fluid dynamics with respect to the relative importance of viscous and inertial forces. For this use, a standard length and a standard velocity must be chosen that somehow characterise the flow. It has been common practice to classify the swimming behaviour of fish by using a value of Re based on body length and forward speed. In the present study, having found the asymptotic mean forward speed, it becomes possible to calculate Re using these customary parameters. Using the given body length and the calculated forward speed, the value of Re would be 3840. It is worth considering which definition of Re more closely characterises the flow. The answer may be that neither gives a fair reflection of the flow regime. As discussed above, a typical ratio of pressure force to drag in this model is 8 (see Fig. 11B). This ratio is effectively an alternative way of defining Re, based on forces due to inertia and to viscosity. In the present case, therefore, the customary definition of Re is approximately 500 times as large as this ratio of forces. The results of this study therefore suggest that viscous effects may be far more important than previously supposed.

The remaining major simplifying assumption in the present work is that the body shape (as opposed to its position) is prescribed. Further development is therefore required before the model can be used to investigate the mechanisms underlying the curious relative timing of muscle stimulation and body curvature in fish (Williams et al. 1989; van Leeuwen et al. 1990; Wardle and Videler, 1993). Such relative timing has important implications for the physiology and energetics of fish muscle (Altringham and Johnston, 1990; Curtin and Woledge, 1993; Rome et al. 1993). The next step in the present approach will be to incorporate into the system a model of the dynamic behaviour of an anguilliform body (Bowtell and Williams, 1991, 1994). Instead of a predetermined body shape, the input to the combined model will then be the known patterns of muscle activation generated by the nervous system. An investigation using this more comprehensive model should increase our understanding of the dynamics of anguilliform swimming and should also provide some insight into why the particular relative timing between muscle activation and muscle shortening has evolved.

In describing the coordinate transformations used here, it is helpful first to consider the general non-orthogonal case, in which the angles and shown in Fig. 4A do not add to zero at the local origin, O. Reducing the general case to an orthogonal system can then follow by making the appropriate restrictions.

Having derived the transformed equations, a brief explanation of the numerical solution procedure is given, together with a test case to show that in a simplified geometry the computed results compare favourably with an analytical solution.

where is the inclination of the grid line. The angle is measured from the global x axis and is positive in the anticlockwise direction, whereas the angle is measured from the global y axis and is positive clockwise. Introducing and in this manner is effectively a separation in which only and only. In the numerical implementation, for a given grid at any time ⁠, a path length along x̄, say Δ x̄, may be obtained from equations A19–A20 by a numerical integration of with respect to x̄ from one grid point to its neighbour. Similarly, may be obtained from equations A21–A22 by an integration of along ⁠.

where the Cartesian translation has been introduced to represent the rate of change in the position of local origin, O, at time ⁠. As the grid point positions change with time, so do the values of and ⁠.

The transformation of first derivatives, already defined in equations A27–A29, is now simplified slightly by setting c=1 to give the orthogonal case.

Each component is in the same direction as the component in equations A27–A29; Ū is the velocity component perpendicular to and is the component perpendicular to x̄.

Equations A53 and A54 are now used to form derivatives with respect to x̄, whereas equations A55 and A56 are used for derivatives. For an orthogonal system, there is no need to distinguish between component U and Ū nor between V and ⁠. Henceforward, the dependent variables U and V will be used, where these are defined such that they are always perpendicular to one another with U in the direction of along x̄ and V in the direction of along ⁠.

In equations A57 and A58, the terms S₁, S₂, S₄ and S₅ all depend on derivatives of angular position. It is these derivatives alone (rather than f and g) that are used to specify the geometry of the curvilinear domain. Once an orthogonal grid has been established, the derivatives of angular position are implicit in the formulation. The use of an orthogonal system gives rise to a succinct formulation (equations A57 and A58) which in turn enables an accurate numerical approximation of the governing equations and the fluid flow.

The numerical variables U_i,j, V_i,j and p_i,j associated with the transformed equations A57 and A58 are depicted in Fig. 4B–D. Fig. 4C shows coordinates centred on U_i,j and applies to equation A57; Fig. 4D, centred on V_i,j, applies to equation A58. Fig. 4B shows the locations of U_i,j, V_i,j and p_i,j centred on a cell within the grid. This serves to illustrate how the continuity constraint in equation 8 is satisfied. Rather than use the transformed version of equation 8, the volume flow into and out of each cell is computed instead, as a divergence. Driving this divergence to zero to satisfy continuity can then follow the solution of equations A57 and A58 as part of the iteration at each time step. A well-established technique for solving such a system is the SIMPLE algorithm (Patankar and Spalding, 1972). That method is used here, with a modification to make the convergence more rapid. It can be shown analytically that this modified SIMPLE method is consistent, in the same sense as for the standard method, in that when the method converges it correctly produces a solution (u,v,p) that satisfies the Navier–Stokes equations. Thus, the numerical solution will be the same whether the standard or the modified SIMPLE method is used. The advantage of using the modified version is that the solution can be obtained to the same specified accuracy as the standard SIMPLE method but with an economy of computational effort.

To assess the degree of accuracy to be expected in a typical application of the present method, a test case has been considered in polar coordinates (r,α). The coordinates are shown in Fig. 15A, together with the associated velocity components. The boundary, r=r_b, has been set to rotate about the origin r=0 at a constant angular speed, with tangential velocity V_α=V_b at r=r_b. This represents, in an approximate sense, the geometry near the tail of the creature as shown in Fig. 3. A part of the full grid in Fig. 3 has been selected and is shown in magnified form in Fig. 15B. This indicates the sense in which the resolution of the grid used in the test case, as shown in Fig. 15C, is comparable to the grid in the full computation. The test case curvature and radial grid size are roughly the same as their counterparts in the full computation. The idealised problem provides a reasonable basis for assessing the accuracy of the numerical approach.

where an arbitrary constant has been set to zero. The profiles for and from equations A82 and A83, respectively, are plotted in Fig. 15D, together with the numerical solution for the grid in Fig. 15C. The degree of agreement between the numerical solution (symbols) and the analytical solution (solid lines) is good. For example, at point 1 ⁠, where the error is greatest, is in error by 1.94 %. The degree of agreement at the other points is considerably better than this.

In general, the numerical solution procedure conserves volume flow in the domain of computation. In this test case, the implication of this is that the pressure gradient, ∂p/∂α, is calculated such that the numerical integration of from to is zero (i.e. Q=0). Because the numerical integration is not precisely the same as the analytical integration, a discrepancy arises between the numerical value of ∂p/∂α and the analytical value corresponding to ⁠. In this example, the value of ∂p/∂α from the numerical solution which is comparable to was found to be 3.9772 at all radial positions. In other words, conserving volume flow leads to an error of 0.57 % in the calculation of the pressure gradient ∂p/∂α which is maintaining the computed velocity profile. Once again, the degree of agreement is good.

The numerical solution in the full computation can be expected to approximate its analytical formulation to a level of accuracy comparable to that found in the test case. In particular, good accuracy can be obtained with a grid which has relatively coarse spatial resolution. In this respect, there is an advantage in having an analytical formulation (such as the present orthogonal system) which represents curvature explicitly when approximating the governing equations in a curvilinear domain.

with v_f(t) positive in the negative x_m(t) direction. The angular velocity, v_r(t), defined in equation 27, has been separated such that only part of the rotational movement, that due to ⁠, remains to be determined.

where v_f(t)^(k) represents the forward velocity between time t−Δt and t at iteration k. The ordinate, y_m(t)^(k), can also be evaluated in this manner, by rearranging equation A85. An expression for the orientation of the body axis, θ_m(t)^(k), can be found by rearranging equation 27, bearing in mind that the summation, which depends on body shape alone and not orientation, will be known. At this stage, for given velocity components [v_f(t)^(k), v_l(t)^(k), v_r(t)^(k) ], the new values of x_m(t)^(k), y_m(t)^(k) and θ_m(t)^(k) are calculated from equations such as equation A86 together with the other constraints of the body geometry. This will completely determine the shape and position of a new body outline.

where ω is a weighting factor commonly chosen to be between 0.5 and 1.0. In this case, numerical experiment has shown that the optimum choice for ω is 0.75. Below this value, the calculation becomes numerically unstable, whereas increasing ω from 0.75 towards 1.0 reduces the accuracy of the numerical approximation.

To assist in understanding the relationship between the positions, velocities and forces in equations A86 and A87, it is helpful to interpret the numerical scheme in the following way. In equation A86, the velocity, v_f(t)^(k), may be regarded more appropriately to represent the value of velocity half-way between the body positions at time t−Δt and t. In the language of numerical analysis, equation A86 could then be regarded as a central difference approximation at time t−Δt/2, which is roughly one order more accurate than it otherwise would be. By the same argument, a weighted force, such as F_f(t) in equation A87, will tend increasingly to provide a better approximation as the value of ω is reduced from 1.0 towards 0.5. However, in common with many iterative schemes, there is a lower limit for stability. Here, as has been stated above, this limit is ω=0.75, corresponding to a time t−Δt/4, precisely half-way between the time (t) at which the new body position is being calculated and the effective time (t−Δt/2) at which the velocity may be thought to occur.

where another weighting parameter, γ_f, has been introduced. This weighting parameter allows a cautious increment in velocity towards the final converged state. In general, the parameter γ may take on different values for the forward, lateral and rotational components of the motion. Extensive numerical investigation has shown that, for optimum convergence, the value of γ should be set to γ_f=0.25, γ_l=0.05 and γ_r=0.05 for the three directions.

With the updated velocity in equation A89, it is now possible to return to equation A86, calculate the new position of the centre of mass and begin the calculation again. When the iterative process has converged (typically in 10 iterations), the three velocities in equation A89 will all have the same values (to within a specified tolerance) and so it can be seen that the numerical solution is independent of the particular value chosen for γ. There is therefore a difference between the role of ω in equation A87 and that of γ in equation A89. Whereas γ controls numerical stability alone and does not affect accuracy, the choice of ω influences both the stability and the accuracy of the numerical approximation.

 
• A

  total mass of body, equation 21

 
• a(i)

  mass of body segment i

 
• b

  amplitude parameter

 
• c

  trigonometric Jacobian, equation A30

 
• D/Dt

  total derivative with respect to t, equation 6

 
• 

  total derivative with respect to ⁠, equations A44, A59 and A69

 
• d₀

  amplitude of the travelling wave at the tail (characteristic reference distance)

 
• F_f(t)

  force acting on the centre of mass in the forward direction, equations 13 and 22

 
• 

  weighted average force, equation A87

 
• F_l(t)

  force acting on the centre of mass in the lateral direction, equations 14 and 23

 
• F₀

  reference force defined in Fig. 11 (not the characteristic force)

 
• f_d(i)

  force acting on body segment i due to drag

 
• f_p(i)

  force acting on body segment i due to fluid pressure

 
• f

  function of (x,y,t) used in x coordinate transformation, equation A1

 
• g

  function of (x,y,t) used in y coordinate transformation, equation A2

 
• I

  moment of inertia, equation 26

 
• i

  body segment (indicates position of body segment relative to the nose)

 
• K

  Jacobian, equation A13

 
• L

  an operator, equation A50

 
• l

  body length along body centre-line (dimensionless)

 
• M(t)

  moment of forces acting on the creature’s body about the centre of mass, equations 15 and 28

 
• N²

  a second-order operator, equation A36

 
• n

  number of body segments

 
• O

  local origin

 
• P

  pressure gradient ∂p/∂α in the test case

 
• 

  normalised P, equation A81

 
• p

  dimensionless pressure

 
• 

  normalised pressure in the test case, equation A83

 
• Q

  mean velocity in the test case, equation A79

 
• Re

  Reynolds number, equation 9

 
• r(i),ϕ(i)

  polar coordinates of body segment i relative to the centre of mass

 
• r,α

  polar coordinates in the test case

 
• r_O,α_O

  local origin in r,α

 
• r_b

  radius of outer boundary in the test case

 
• 

  normalised radial coordinate in the test case, equation A78

 
• S₁, S₂, …, S₇

  symbols representing combinations of derivative terms, equations A61–A66

 
• s

  distance along body centre-line from nose (dimensionless)

 
• T₁, T₂

  second derivative terms, equations A34 and A35

 
• t

  dimensionless time

 
• t₀

  locomotor cycle duration (characteristic reference time)

 
• 

  time in transformed coordinates ⁠, equation A3

 
• U, V

  velocity components in the x̄ and directions, respectively, equations A51 and A52

 
• 

  velocity components perpendicular to and perpendicular to x̄respectively, equations A46 and A47

 
• u, v

  dimensionless velocity components in the x and y directions, respectively

 
• V_r, V_α

  velocity components in the r and α directions, respectively, used in the test case

 
• V_b

  tangential velocity of the outer boundary in the test case

 
• 

  normalised tangential velocity in the test case, equation A77

 
• v_f(t), v_l(t)

  velocity components of the centre of mass in the forward and lateral directions, respectively, equations 19 and 20

 
• 

  velocity component corresponding to ⁠, equation A88

 
• v_r(t)

  rotational velocity component, equation 27

 
• w_s

  body width at distance s from the nose, equation 3

 
• w₀

  body width at the nose

 
• w₁

  body width at the tail

 
• 

  resolved components of the rate of change in position of the local origin, O, equations A31 and A32

 
• x,y

  Cartesian coordinates (dimensionless and fixed in the enclosure)

 
• x_m(t),y_m(t)

  position of the creature’s centre of mass in x,y

 
• x_s,y_s

  Cartesian coordinates defining body shape independently of the enclosure (see equations 1 and 2)

 
• 

  transformed coordinates, equations A1 and A2

 
• 

  components of the rate of change in position of the body segment i in the x and y directions, respectively

 
• 

  components of the rate of change in position of the local origin, O, in the x and y directions, respectively, equations A25 and A26

 
• ŷ_sy_s

  scaled by a factor β

 
• α

  see r,α

 
• α(i)

  angular position of body segment i relative to the x axis (positive anti-clockwise), equation 12

 
• 

  rate of change of α(i) with respect to time

 
• β

  amplitude scaling factor, equation 30

 
• γ_f, γ_l, γ_r

  weighting parameters for velocity in the forward, lateral and rotational directions, respectively

 
• Δs

  body segment length (along body centre-line)

 
• Δt

  time step size

 
• θ_m(t)

  angular position of the x_s axis relative to the x axis (positive anti-clockwise)

 
• 

  inclination of the x̄ grid line relative to the x axis (positive anti-clockwise)

 
• μ

  dynamic viscosity of water

 
• ρ

  density of water

 
• ϕ (i)

  see r(i),ϕ(i)

 
• 

  rate of change of ϕ(i) with respect to time

 
• 

  inclination of the grid line relative to the y axis (positive clockwise)

 
• ω

  weighting factor for force, e.g. for F_f(t) in equation A87

 
• ∇²

  Laplacian operator, equations 7, A33 and A41

 
• 

  Laplacian operator in transformed coordinates, equations A45, A60 and A70

i, j integer pair depicting a cell in the computational grid (see Fig. 4)

o a variable at a previous time step (see equations 16–18, 25 and 28)

Superscript

(k) iteration level (see equations A86–A89)

The authors wish to thank Mrs Ella Town for her assistance in preparing the manuscript. This work was supported by the Biotechnology and Biological Sciences Research Council (UK).