## Abstract

We provide initial evidence that a structure formed from an articulated series of linked elements, where each element has a given stiffness, damping and driving term with respect to its neighbours, may ‘swim’ through a fluid under certain conditions. We derive a Lagrangian for this system and, in particular, we note that we allow the leading edge to move along the *x*-axis. We assume that no lateral displacement of the leading edge of the structure is possible, although head ‘yaw’ is allowed. The fluid is simulated using a computational fluid dynamics technique, and we are able to determine and solve Euler–Lagrange equations for the structure. These two calculations are solved simultaneously by using a weakly coupled solver. We illustrate our method by showing that we are able to induce both forward and backward swimming. A discussion of the relevance of these simulations to a slowly swimming body, such as a mechanical device or a fish, is given.

## 1. Introduction

In this article, we wish to show by using a numerical approach that a simple articulated structure may propel itself through a fluid, and we provide preliminary investigations to illustrate our method. We shall assume that the fluid may be viewed as two-dimensional (2D) continuum (such as in a soap film) and that our structure has no depth, and so may be assumed to be one-dimensional (1D). We restrict our simulations to those cases in which the motion of the mechanical structure lies in the 2D plane. Another explicit assumption of the calculations is that the head of the structure always lies on the *x*-axis and points in the direction of motion. We wish to present the method in this article and illustrate this approach by looking at cases that have relatively low Reynolds number (*Re*) of approximately 250. We note that this choice for *Re* aids numerical stability, and we intend to consider the more challenging situation of higher *Re* in future calculations.

Experiments concerning the properties of soap films have recently been carried out (Martin & Wu 1995; Chomaz & Cathalau 1996; Rutgers *et al*. 1996; Zang *et al*. 2000), and these systems have been proposed as the experimental versions of theoretical 2D liquids. Indeed, we note that particularly interesting fluidic properties were observed when a silk filament was introduced into the flowing soap film (Zang *et al*. 2000). Such flapping filaments or flags have been pointed out to have a clear relevance to swimming fishes (Huber 2000).

We note that problems that are conceptually similar to the problem of our ‘swimming structure’ have previously been treated with considerable success by using a so-called ‘weakly coupled’ solver. We shall adopt this approach here, and examples of such problems are a single filament (Farnell *et al*. 2004*a*) and coupled filaments in a flowing soap film (Farnell *et al*. 2004*b*), and the motion of the leaflets in artificial heart valves (Horsten 1990; Fenlon & David 2001*a*,*b*; Farnell *et al*. 2004*c*). Our approach has therefore been extensively tested for these systems; thus, we shall not extensively test the methodology again here, and the interested reader is referred to these references for more information. Furthermore, the validity and accuracy of our approach has also been proven in these extensive earlier calculations. We now wish to extend the range of applicability of our approach and hope only to show that, in theory, a simple articulated structure is able to propel itself through a fluid.

We note that two main methods have been used in the past to investigate theoretically the motion of elongated animals in fluids. Taylor's theory (Taylor 1952) was based essentially on a resistive model, whereas that of Lighthill (1960) used a ‘reactive’ analysis, where the investigation concentrated on the reactive forces between a small volume of fluid and the animal's surface in contact with the fluid. However, a review by Lighthill (1969) showed that the resistive method seemed to provide the best results for elongated vertebrae. The swimming action of fishes has also been studied both experimentally (Gray 1933, 1955; D'Août & Aerts 1997, 1999; Drucker & Lauder 1999; Sfaliotakis *et al*. 1999) and theoretically (Gray & Hancock 1955; Williams *et al*. 1995; Liu *et al*. 1996; Sigvardt & Williams 1996; Carling *et al*. 1998; Liu & Kawachi 1999; Long *et al*. 2002; Huber 2000).

We begin our treatment of this system by firstly considering our mathematical and numerical methodology in detail. We then illustrate our methodology by presenting results for our initial studies of this system, and show that we may induce either backward or forward motion by simply changing the form of the driving terms acting along the length of the mechanical structure. We then discuss the relevance of our studies to swimming action of fishes and/or mechanical devices (Sfaliotakis *et al*. 1999). We conclude by summarizing our results and considering other possible future applications.

## 2. Method

### 2.1 The Navier–Stokes equations

We assume that the fluid is incompressible, with a viscosity that is Newtonian in form, and that the density is constant. In addition, we use similar assumptions to those used in experiments for measuring the viscosity of soap films (Martin & Wu 1995). The flow is treated as being isothermal and laminar, and an argument for this assumption is given below. The relevant time-dependent equations are thus given by the Navier–Stokes equations,(2.1a)

(2.1b)

Here, is the 2D vector field corresponding to the *x*- and *y*-directions as shown in figure 1, is the fluid pressure, *ρ* is the density of the fluid and *ν* is the kinematic viscosity of the fluid. We non-dimensionalize the fluid equations such that(2.2)where *L* is the actual (dimensional) length of the mechanical structure and *U* is a (dimensional) characteristic speed (e.g. the swimming speed of the mechanical structure). Equation (2.1*a*,*b*) now becomes(2.3a)(2.3b)where the *Re* is given by . We assume the fluid domain is as given in figure 1; namely, a channel with an inlet and outlet. Non-slip boundary conditions are assumed on the walls of the channel, and ‘natural’ (or stress-free) boundary conditions are assumed at the inlet and outlet. We note that the width *h* of the channel is also set to one (non-dimensional) unit. The width of the mechanical structure is assumed to be 0.04 non-dimensional units and we set the *Re* to be *Re*=250, as this aids numerical stability (see §4). The length *H* of the channel is 15 (non-dimensional) units, far enough away from the mechanical structure so that the outlet boundary has negligible effect. The mechanical structure is placed initially approximately at the mid-point of the channel. Details (see also Fenlon *et al*. 2001;Farnell *et al*. 2004*a*,*b*) of performing a computational fluid dynamics (CFD) calculation are presented in §5. A fuller investigation of the effects of varying the values for these parameters on our calculations will be considered in a future article.

## 3. The Euler–Lagrange equations

Similar to earlier calculations (Williams *et al*. 1995, Sigvardt & Williams 1996; Carling *et al*. 1998), we assume that the mechanical structure is composed of *N* equal length, homogeneous, rigid elements. The (dimensional) length of each element is and the (dimensional) mass of each element is . However, in our case, the elements are fixed to one another at their hinge (or fulcrum) points, and the mechanical structure is thus assumed to be approximated by a form of *N*-tuple pendulum, in which each hinge has a positive spring stiffness coefficient and a positive damping coefficient associated with it. We assume that no lateral displacement of the head in the *y*-direction may occur where this assumption is supported by experiment (Gray 1933, 1955). However, head yaw is allowed, as the angle that the first element makes with the *x*-axis does not have to be zero. We note that the leading edge of the mechanical structure is allowed to move in the *x*-direction only. Each element *s* makes an angle *θ*_{s} with the horizontal *x*-axis and the position of the leading edge of the mechanical structure on the *x*-axis is given by . In analogy with previous studies (Farnell *et al*. 2004*a–c*), the Langrangian *Ψ* may be shown to be given by(3.1)where the terms due to the stiffness at each hinge are valid only for small values of (and that ). Note that the subscript *s* for the element angles runs from 1 to *N*, although we may define—notationally only—an angle *θ*_{0}, which aids our definition of the potential energy in equation (3.1). (However, we do assume that the minimum potential energy position is given by a profile in which the angle for the leading edge element is zero.)

The Lagrangian given by equation (3.1) may now be used to obtain the Euler–Lagrange equations by writing(3.2)where and represent non-conservative external torques and forces, respectively. The function represents the externally applied torques exerted on a given element *s*, and these torques are non-conservative. We subdivide the terms within this function into two distinct pieces, given by(3.3)

We note that contains the driving-term to our equations and corresponds to the ability of the mechanical structure to cause torques along its length (see §4 for more details). The damping terms are thus encoded in the terms , and favour a difference between angular velocities of successive elements (both above and below a particular element *s*) that is small. An inflexible system is thus approximated by having large values for both and . The term in equation (3.3) refers to external pressure forces (as distinct from the internal stiffness and damping terms or, indeed, the driving terms) acting on the mechanical structure elements creating a torque of strength , which in this case are caused by the fluid pressures forces generated by the liquid. For the present case, we can show from an order of magnitude argument that the pressure forces are much larger than those exerted by viscosity and hence the viscous shear forces acting parallel to the *N*-tuple pendulum are considered negligible.

We also note that the non-conservative forces for the leading edge generalized variable are given by(3.4)

We solve these equations using adaptive step-sized Runge–Kutta method (see Farnell *et al*. 2004*a–c* for more details). Equation (2.3*a*,*b*) is solved using a commercial code (FIDAP). Although our solutions are mesh independent, we do not fully investigate the variation of our solution with respect to the width of mechanical structure and *Re*. We wish to present only the initial calculations and thereby investigate whether our new approach might be used to give reasonable results. We now non-dimensionalize the Euler–Lagrange equations by introducing (as before) characteristic scales for mass, length, velocity and time, such that(3.5)

The quantities , and may also be non-dimensionalized in a similar manner to give(3.6)

The non-dimensionalized Euler–Lagrange equations for the element angles are given by(3.7)and the non-dimensionalized Euler–Lagrange equation for the leading-edge term is given by(3.8)

The pressure force acting on an element of the mechanical structure in the dimensional system in our Euler–Lagrange equations (namely ) must be related to the pressure difference (between those pressures acting on the bottom, , and the top, , of a mechanical structure element) such that(3.9)

The non-dimensional form of this equation is given by(3.10)where *ρ*_{f} is the 2D density of the fluid. Furthermore, we note that we may express the mass of the mechanical structure in the experimental (dimensional) system in terms of a *linear density*, *σ*, such that . Again, note that we assume that our idealized mechanical structure is 1D.

## 4. Values for the parameters

We choose values for the parameters such that they are consistent with earlier simulations for different, although related, systems (Fenlon & David 2001*a*,*b*; Farnell *et al*. 2004*a–c*). In particular, we choose a value for *F* of unity. As we wish to prove the principle of whether we are able to simulate any of the observed properties of the system, we feel that this is a reasonable assumption. In a similar manner, we set the *Re* to be *Re*=250 as this aids numerical stability. Note that inertial forces are dominant and viscous forces are usually neglected for 10^{3}<*Re*<10^{5} (Sfaliotakis *et al*. 1999).

The fluid is initially at rest and non-slip boundary conditions are imposed on the walls of the channel, while the inlet and outlet have stress‐free boundary conditions. We choose a fluid mesh such that there are at least four mesh nodes within the width of the mechanical structure. All of the parameters above are chosen to be in reasonable agreement with those experimental values given elsewhere (D'Août & Aerts 1997, 1999; Drucker & Lauder 1999; Sfaliotakis *et al*. 1999).

The values for the stiffness, damping and driving terms are chosen to induce either forward or backward swimming. For the case of the backward swimming mode, we keep the stiffness and damping terms constant, where(4.1)and we set the driving terms to have a travelling wave form, given by(4.2)

In this case, *A*=0.03 and *ω*=1.0. The stiffness and damping terms for the forward swimming mode are set to be constant but of greater magnitude (to reduce head yaw), where(4.3)

The driving terms for the forward swimming mode are given by(4.4)with *A*=0.2 and *ω*=1.0. We see from equation (4.4) that the driving terms of the mechanical structure are much stronger at the tail than at the head for the forward swimming mode. Our calculations were tested for *λ*=0.25 and *λ*=1.5, and qualitatively similar results were seen for both cases for backward and forward swimming. Results quoted below (especially in figures 2–6) are for the case *λ*=0.25.

## 5. Computational solution

FIDAP (a finite element code used to solve the Navier–Stokes equations) has a moving boundary condition application that can be used to apply a superimposed moving boundary internally in the fluid domain. It continues to solve over the full fluid domain, including that superimposed by the moving boundary, but a user-defined subroutine is used to impose boundary conditions to the elements and nodes contained within the imposed moving boundary: in this case, the mechanical structure. The situation is displayed graphically in figure 2.

The structure, given an initial position and defined by a fixed end of the first rod, has a finite thickness, *δϵ*, is considered about the rods, and it is on elements within this thickness that boundary conditions are imposed. For instance, in figure 2, the downward movement of the structure results in velocities normal to the movement (arrows) being applied to the nodes that lie within the specified thickness surrounding the leaflet rods.

At each time-step of the fluids solver, for every node that falls within a specified moving region of the modelled fluid domain, FIDAP solves the Euler–Lagrange equations via the user subroutine. The pressure (calculated by FIDAP) at the elements that are situated within a certain distance above and below each rod is passed to the user subroutine and used for calculation of the movement. Using the structure motion information, FIDAP then solves the Navier–Stokes equations at the next time-step, with the added moving boundary constraints placed on the nodes that fall within the moving domain. This means that the interaction of the fluid and leaflet lags the fluid solution by one time-step; however if time-steps are made small enough in FIDAP, then this approximation has been shown to be accurate to within the error of the fluids solver.

The algorithm essentially uses sub-cycling for each time-step for the fluids solver. The Euler–Lagrange equations are solved via an LU decomposition and a fourth-order Runge–Kutta algorithm, with each Runge–Kutta time-step at least an order of magnitude smaller than the fluid time-step. We use a weakly coupled solver for this system. Note that the total (summed) error for the Runge–Kutta technique for the position *x*_{0} and angles (and their velocities) was kept at 10^{−8} by using adaptive step-sizing.

The total error for the velocity and pressure values was set to 10^{−3}, which was previously found to be adequate for related and similar simulations for the related study of flapping flags and heart valve leaflets (Farnell *et al*. 2004*a–c*). The fluid mesh was formed of regular square-shaped elements in the region in which the structure could swim, and we note that at least four nodes were contained in the width of the mechanical structure (namely, 0.04 non-dimensional units). This approach has previously been shown to be enough to adequately simulate similar structures (flapping filaments and leaflets in heart valves) in articles by Farnell *et al*. (2004*a–c*) and Fenlon & David (2001*a*,*b*).

## 6. Results

We may induce forward swimming by setting the amplitude of the driving terms in equation (4.4) to be smaller at the head of the mechanical structure than at the tail (see §4; for more details). Figure 3*a* illustrates the velocity vector field for a typical profile for the forward swimming case (see also figure 5 for the forward swimming case). We note that vortices are clearly created at the trailing edge of the mechanical structure as the latter end of the structure (which we shall refer to as the ‘tail’) moves from side to side. Backward swimming in figure 3*b* may be induced by keeping the magnitude of the stiffness and damping terms constant for each element angle, given by in equation (4.1), and also by keeping the amplitude of the driving term *A* in equation (4.2) constant (see §4 for more details). We note that there are net pressure differentials acting along the *x*-axis that cause the mechanical structure to swim either backwards or forwards (see below). These pressure differentials are created, in circular fashion, by the swimming action of the mechanical structure over the course of the simulations.

Results for a typical profile seen in this simulation are given in figure 3*b*, in which velocity vectors for the moving fluid are also presented, and it can be seen that vortices are built up along the length of the mechanical structure. Clearly, however, these vortices will again eventually dissipate owing to viscous effects in the fluid. An important point to note is that we are fully able to simulate the fluid properties, such as the fluid pressures and the velocity vector field in the whole of the channel. An advantage of our approach is that we are able to simulate this system (e.g. the pressure scalar field of the fluid) from an *a priori* point of view.

Instantaneous streamlines for typical profiles are plotted in figure 4. We see that complex vortex structures are built up by the mechanical structure for both the backward and forward swimming modes. We note that the streamlines indicate a different type of vortex structure for forward swimming than for backward swimming. Figure 5 shows figures for the vortex structure at the trailing edge of the structure for the forward-swimming mode in more detail.

An alternative method of inducing forward swimming is to make the stiffness and damping terms much higher at the head of the mechanical structure than at the tail and to keep the magnitude of the element‐driving terms constant, as in equation (4.2). Indeed, this was also seen to lead to a forward swimming mode. Again, we note that we do not wish to perform an exhaustive survey of the (rather large) parameter space at this early stage. Rather, we wish to identify cases in which regular swimming motion (both backwards and forwards) is observed, and thus illustrate our method. We believe that our approach is able to capture at least some of the essential physics occurring in a true swimming system (see below). Further investigations might consider the efficiency of the various swimming modes, although this is not considered here.

The evolution of the profiles of the swimming mechanical structure with respect to time is investigated in figure 6. Figure 6*a* illustrates a number of typical profiles for the mechanical structure for the forward swimming mode. The forward motion of the leading edge of the mechanical structure is now clear. We note that after an initial phase in which the amplitude is large, the amplitude of oscillation remains roughly constant throughout the 25 non-dimensional time units of the simulation. The amount of head yaw is thus initially large, but then settles into a lower value once the system has reached the stable swimming mode.

Figure 6*b* illustrates a number of typical profiles for the backward swimming mode of the mechanical structure. The backward motion of the leading edge of the mechanical structure is clear. We also note that the amplitude of oscillation remains approximately constant throughout the 25 non-dimensional time units of the simulation, and that this amplitude for backward swimming is larger than the case for the forward swimming mode throughout the simulation.

Figure 7*a* plots the positions of the centre of mass of the mechanical structure in the forward and backward swimming modes. We see that an initial phase of little motion is observed, and that the speed of the swimming mechanical structure, both backward and forward, increases with time. The two cases are still accelerating by the end of the simulation, although we would expect the mechanical structure to reach a terminal velocity at some point at which the fluid dissipation of energy would balance that fed into the system from the mechanical structure's driving terms.

Figure 7*b* presents the *y*-values (in non-dimensionalized units) for the trailing edge of the mechanical structure with respect to time. Interestingly, the forward swimming mode appears to have more tail flutter than the backward swimming mode, although it is unclear whether this is a generic characteristic or is caused by the particular choice of parameters for the forward swimming mode. Furthermore, forward swimming seems to take a little longer to reach its stable oscillating mode than backward swimming, and the backward swimming mode has a larger amplitude tail motion than the forward swimming mode in their respective stable oscillating modes.

Results for the element angles in figure 8 clearly show that the amount of head yaw, indicated by the magnitudes of the first element angles, is reduced for forward swimming compared with backward swimming. This is in agreement with previous experiments (e.g. D'Août & Aerts 1997, 1999; Sfaliotakis *et al*. 1999). (Note that the leading edge of the structure is constrained to lie on the line *y*=0.)

Figure 9 presents a plot of the Fourier transforms of the time evolution of the values of the element angles of the trailing edge of the mechanical structure with respect to (non-dimensional) frequency. This plot shows a strong peak at one unit of frequency (non-dimensional) for both forward and backward swimming. This is to be expected as we have set *ω*=1 in equations (4.2) and (4.4). Interestingly, both forward and backward swimming also demonstrate lower frequency modes. This mode occurs at a slightly larger frequency value (and is relatively larger in amplitude compared with the peak at one unit of non-dimensional frequency for the forward swimming case) than for the backward swimming case. In both cases, the lower frequency peak is more pronounced for larger values of *s* (i.e. towards the tail rather than at the head of the mechanical structure). It is unclear why a lower frequency mode should occur although we note that this might be a ‘resonance mode’ of the *N*-tuple pendulum itself owing to the particular values for the stiffness and damping terms. This would also explain why the peaks occur at different frequency values because the stiffness and damping coefficients have been increased for the forward swimming mode compared with backward swimming.

The backward swimming mode also shows side-bands near to 1 cycle unit^{−1} time. One could speculate that this behaviour is linked to the direction in which the mechanical structure will swim, and also the lower efficiency of swimming observed in experiment for backward swimming, although this issue is beyond the scope of this preliminary article.

## 7. Comparison to experiments

We shall now try to connect our results to those of previous experiments. Clearly, our value for the *Re* of 250 is much lower than that of a real swimming (e.g. anguilliform or eel-like) adult fish, which may have values of 5000 *Re* or higher. We note, however, that our value for the *Re* of 250 might correspond to larvae, very slow-swimming or very young fishes. Our results suggest that a slow-swimming articulated mechanical device (Sfaliotakis *et al*. 1999) might be constructed by using motors and springs at the hinges for the driving terms and stiffness terms, respectively, and that this device could swim both backwards and forwards.

The behaviour of the fluid for adult anguilliform swimming fishes was seen here to be somewhat different from that observed in experiment (Müller *et al*. 2001; Tytell & Lauder 2004). In particular, the vortices were seen to be long-lived in these experiments for forward swimming, whereas our vortices were seen to dissipate quickly. However, vortices that were seen to be created at the tail end of fishes in both of these experiments were also seen in our preliminary calculations, although they were not as long-lived in our calculations. Furthermore, we argue that a similar streamline structure as seen in figures 3*a*, 4*a* and 5 is also seen in a previous experiment (Müller *et al*. 2001). It would be interesting to see if even better correspondence with past experiments might occur for further simulations at higher values of *Re*. Note that our results do not appear to show the large stable re-circulation zones behind the tail of the fishes seen in numerical calculations by Carling *et al*. (1998), although the validity of such large stable re-circulation zones remains to be conclusively proven (Tytell & Lauder 2004).

We may also consider the typical shapes that eels make when swimming and compare these with experimental results (Gillis 1996, 1997, 1998*a*,*b*). Our results tend to show higher amplitude motion towards the front of our structure or ‘fishes’ in the forward swimming mode than those seen in experiment (Gillis 1996, 1997, 1998*a*,*b*). This may be owing to the particular values for the driving, damping and stiffness terms, and the fact that the *Re* chosen here was not a perfect fit for an anguilliform fish. For example, we choose *λ*=0.25 in equations (4.2) and (4.4), and this may well be too small for anguilliform fishes in particular. However, figures 6 and 9 also indicate that the head yaw for forward swimming is reduced when compared with backward swimming, as is seen in previous experiments (e.g. D'Août & Aerts 1997, 1999; Sfaliotakis *et al*. 1999). We also note that calculations in which the stiffness terms in the structure were much stronger at the head of the structure than at the end, and in which the driving terms were of constant amplitude (as for the backward swimming case), also led to forward swimming with even greater relative reduced amplitude at the front of the structure.

It must be said that the initial numerical simulations provided here are such that the mechanical structure is still accelerating to its equilibrium value. The results of Taylor's earlier works (Taylor 1952) for forward swimming using a resistive analysis show that(7.1)where *V* is the velocity at infinity relative to the material of the sheet or mechanical structure. The sheet thus moves with velocity −*V* relative to the fluid at infinity when waves of lateral displacement travel with velocity +*U* down the sheet. Importantly, the analysis indicated (using a swimming diagram) that the ratio of the equilibrium velocity of the mechanical structure to the velocity of the travelling wave could not exceed unity.

Our results indicate that the mechanical structure after 20 non-dimensional seconds has reached a velocity ratio of magnitude approximately 0.1 in both the forward and backward swimming cases. Numerical results by Carling *et al*. (1998), where the mechanical structure's motion was prescribed (but its axial and lateral velocity was unknown), showed that the velocity ratio was 0.7. Clearly, the velocity ratio of the simulations presented here will saturate at some value above 0.1, and so our results are at least consistent with these earlier results. Future calculations will enable a more detailed comparison of our results with these previous simulations.

It is difficult to compare the results for our backward swimming results to experiment as this case has previously been considered to a lesser extent. However, our simulations for backward swimming were carried out with the same assumptions as for the forward swimming case. The only parameters changed were internal parameters for the driving, stiffness and damping terms in the Euler–Lagrange equations for the structure, and the parameters for the fluid were not changed. We see radically different vortex structures in this case. Indeed, we observe vortices built up along the length of the mechanical structure for backward swimming. It would be interesting to determine whether similar structures are seen when the structure is at constant velocity (and not accelerating) and/or in experiment.

## 8. Conclusions

We have outlined a theoretical/numerical procedure for simulating the swimming action of a 1D mechanical structure. This structure is an *N*-tuple pendulum that has appropriate driving terms as well as stiffness and damping terms at each hinge. We use a combination of CFD to solve the fluid and Euler–Lagrange equations for the mechanical structure. By using a weakly coupled solver, we are able to solve both the mechanical structure and fluid on an ‘equal footing’. In particular, we are able to determine the fluid pressures due to the motion of the mechanical structure (itself driven by the driving terms), which cause the whole of the mechanical structure to move through the fluid. We have identified stable forward and backward swimming modes.

We note that the simulations presented here are merely preliminary and more detailed calculations are necessary. For example, the effects of changing the simulation parameters, such as mesh fineness and mechanical structure width, and of changing the fundamental parameters governing the fluid properties and mechanical properties of the mechanical structure need to be considered. Undoubtedly, there are many different values for the parameters that lead to both forward and backward swimming, and mapping out the properties of this coupled system would be interesting.

We have discussed the potential relevance of our procedure in order to study the swimming action of a mechanical device or swimming fish from the point of view of such a numerical model. The above calculations suggest that a slow-swimming articulated mechanical device might be constructed by using motors and springs at the hinges for the driving terms and stiffness terms, respectively.

Note that our value for the *Re* of 250 is much lower than that of a real swimming (e.g. anguilliform or eel-like) fish, which may have values of 5000 *Re* or higher. For *Re* in this range, we would expect a much different picture for the fluid behaviour to emerge; for example, we might reasonably expect the vortices thus created to exist for longer times than those seen here.

It would therefore be interesting to determine whether our approach might be extended to such higher Reynolds numbers so as to make more exact mappings to actual values in a real-life scenario. We might thus be able to study which of the two modes (backward or forward swimming) would be the most efficient. We note that experimental investigations suggest that forward swimming is the more efficient and that backward swimming is used infrequently (e.g. in reaction to a predator).

## Acknowledgments

We acknowledge and thank Dr G. Blyth and Mr E. Allwood of University of Leeds for their technical advice regarding UNIX and their support in running FIDAP simulations.

## Footnotes

- Received September 24, 2004.
- Accepted January 17, 2005.

- © 2005 The Royal Society