Experimental approach

Our experimental design is guided by the goals of achieving self-propulsion for prescribed flapping motions and realizing an effectively infinite array of bodies to explore their interactions. As shown in Fig. 1a and detailed in the Methods section, the device consists of one or more wings or hydrofoils that are heaved up and down, and consequently swim in orbits around a cylindrical water tank33,34,35. A video showing the device in operation is available in Supplementary Movie 1. Free swimming is achieved by mounting the wings to a vertical axle via a low-friction rotary bearing. Thus, the forward locomotion is not prescribed but is an outcome of the interaction with the fluid. This self-propulsion condition as well as the inertial or high Reynolds number flows are important properties shared with bird flight and fish swimming. Our system, however, involves fixed separations between bodies, whereas inter-individual spacing is free or dynamically determined in animal groups, and this difference will likely lead to differences in the emergent locomotion.

Figure 1: Flapping wings swimming in rotational orbits mimic an infinite array of locomotors. (a) A motor heaves a wing or wing pair up and down at prescribed frequency f and peak-to-peak amplitude A, resulting in swimming of rotational frequency F around a cylindrical water tank. (b) The rotational geometry allows for interactions with the flows generated in previous orbits. Full size image

Previous studies have shown that the dynamics and flows observed in rotational systems compare well with those in translational geometries33,34,35,36. To ensure the phenomena we observe are not a peculiarity of geometry, we accompany our three-dimensional (3D) rotational experiments with computational simulations in two-dimensional (2D) translational domains, as discussed in detail below. Most importantly, the rotational geometry induces each wing to swim within the wakes of others, as indicated by the trajectories of Fig. 1b, and this arrangement thus mimics the conditions in an infinite array. The device also allows us to explore swimming behaviour for varying strength of hydrodynamic interactions between swimmers. Increasing the number of wings connected to the assembly decreases the spacing between members in an array and amplifies interactions. Further, increasing flapping frequency or amplitude induces faster swimming and thus stronger flows to be encountered in subsequent orbits. The swimming dynamics are characterized by measuring the rotational speed around the tank, and use of a clear-walled tank allows for visualization of the flow fields.

Dynamics of interacting wings

Previous studies have used rotational systems to study the dynamics and hydrodynamics relevant to an isolated wing33,34,35. This is achieved by relatively slow swimming of widely spaced hydrofoils, in which case interactions are weak and the swimming speed follows simple dependencies on flapping frequency and amplitude. Here we aim to induce strong interactions by considering fast swimming of a wing pair, and indeed we observe markedly different dynamics. To systematically characterize the locomotion, we vary both the peak-to-peak amplitude A and flapping frequency f for the pair and measure the resulting rotational frequency F around the tank. For each A, we incrementally increase f, and this upward sweep is followed by a downward sweep to low values. The data of Fig. 2a show that faster flapping leads to faster swimming, although F(f) does not always change continuously. For A=10 cm, for example (red curve), F increases smoothly until a critical frequency near f=1 Hz, at which F abruptly increases, in this case doubling its speed. As f is increased further, F again increases continuously. For decreasing f, F remains high before abruptly dropping near f=0.2 Hz. Thus, the measured dynamics display a hysteresis loop, and different flapping amplitudes lead to other loops. This hysteretic behaviour is associated with bistability of states: for identical flapping kinematics, the array can take on one of two ‘gears’ corresponding to slow and fast swimming modes.

Figure 2: Dynamics of interacting wings. (a) Swimming speed, as measured by the rotational frequency F, versus flapping frequency f. For each peak-to-peak amplitude A, an upward sweep of f is followed by a downward sweep (as indicated by arrows), and the data form a hysteresis loop. (b) Schooling number S, which represents the number of wavelengths separating successive wings. Each hysteresis loop is bounded by in-phase (integer value of S) and out-of-phase states (half-integer S). (c) A polar histogram of S mod 1 shows peaks corresponding to preferred spatial phase relationships between successive wings. Full size image

To organize these observations, we first note that the quantity f/F represents the number of flapping strokes in one orbit of the tank, and thus the number of strokes separating the two wings is given by S=f/2F. Equivalently, this schooling number S represents the separation distance between successive swimmers measured in wavelengths of motion, and this quantity serves as a general way to characterize arrays of synchronized swimmers. Specifically, S encodes the spatial phase shift between successive swimmers: whole integer values of S denote spatially in-phase states in which neighbouring wings trace out the same path through space, and half-integer values indicate out-of-phase states. In Fig. 2b, we show how S depends on f for all the data of Fig. 2a. At each amplitude, low f leads to slow swimming or large S, which indicates weak interactions. Increasing f leads to decreasing S followed by saturation near a whole integer value, thus residing nearly in-phase for a broad range of frequencies. Further increasing f induces an abrupt downward jump in S. Conversely, when f is decreased, the system approaches half-integer S (out-of-phase) before an abrupt upward jump.

Viewing all amplitude data sets together reveals preferred and avoided ranges of S for strongly interacting wing arrays. We quantify this observation in Fig. 2c, in which the histogram of S mod 1 reveals the distribution of spatial phase shifts. Values between 0.5 and 1 are almost entirely absent, and among the prevalent phases, distinct peaks appear at values of 0, 0.25 and 0.5. Thus, these swimmers preferentially occupy spatial phase relationships corresponding to purely in-phase, leading by a quarter cycle, and purely out-of-phase. The appearance of these particular phase relationships suggests that schooling in this system involves bodies interacting coherently through their wave-like flow fields.

Flow visualization

To shed light on the hydrodynamics of these modes, we seed the fluid with particles and use high-speed video and image analysis to extract the flow fields (Supplementary Movies 2 and 3). First, we consider a single wing swimming at relatively low amplitude and frequency, in which case interactions are weak. Indeed, the flows produced are largely dissipated between successive passes, and thus the wing swims into quiescent fluid, as shown by the velocity vector field of Fig. 3a and the schematic of Fig. 3b. The wing leaves behind a stream-wise array of counter-rotating vortices that are staggered vertically. This so-called reverse von Karman wake has been observed in biological and physical experiments, and is a signature of thrust production during flapping33,34,35,36,37,38. The vortices are separated by strong up and down flows, which we highlight in Fig. 3a using a colour map of the vertical component of velocity. Each upstroke generates an upward flow (red) and each downstroke a downward flow (blue), and thus the wake can be viewed as a wave in space that reflects the trajectory of the wing.

Figure 3: Flow visualization. (a,b) Flow field around a non-interacting wing of chord length 3 cm extracted using particle image velocimetry and rendered in a schematic. The colour map indicates the vertical component of the velocity vector field, with red indicating upward and blue downward flows. The upstroke produces an upward flow and the downstroke a downward flow. (c,d) Slow mode of interacting wings: the downstroke of a wing (red path) occurs within the downward flow of its predecessor (blue path). (e,f) Fast mode: the downstroke occurs within the upward flow of its predecessor. Full size image

To induce stronger interactions, we attach multiple wings to the assembly but maintain the same flapping motion. Under these conditions, the system exhibits a hysteresis loop associated with the slow and fast swimming modes. Figure 3c and d show the flow field for the slow mode, in which successive wings trace out similar paths through space. Captured here during the downstroke, a wing is seen to sink within the downward (blue) flow generated by its predecessors. Subsequently, the upstroke occurs within an upward (red) flow. This visualization reveals that not only are successive wings spatially in-phase with one another but also any given wing is also in-phase with the existing flow structure into which it swims. Thus, this constructive wing–wake interaction mode involves the repeated reinforcement of existing flows. Conversely, the fast mode involves antagonistic motions between each wing and the flow it encounters. As shown in Fig. 3e,f, the downstroke of the wing moves against the upward flow of its predecessor, and the subsequent upstroke occurs within a downward flow. This destructive mode thus involves the inversion of flow fields, in which the passing of a wing replaces an existing upward flow with a downward one, for example.

Simulations

To gain additional insights into the intrinsic dynamics of swimmer arrays, we conduct computational fluid dynamics simulations that solve for the flow field and locomotion of a flapping body. Similar to the experiments, these simulations involve free swimming of a wing, where the emergent speed reflects hydrodynamic interactions (see Methods section as well as Supplementary Movies 4 and 5). The swimming dynamics is determined from the computed fluid forces, and an infinite array of swimmers is replicated by having a single wing repeatedly traverse a domain of length L with periodic boundary conditions39. Unlike the experiments, the flow field is 2D and the swimming motion is translational rather than rotational. Thus, observations common to both experiments and simulations are expected to be generic and not due to system dimensionality or geometry. We systematically explore the dynamics through a procedure similar to that of the experiments: Flapping frequency is incrementally increased, with each step initialized with the final data of the previous step and then allowed to reach a terminal swimming speed. This upward sweep of frequency is then followed by a downward sweep.

These simulations reveal a number of features in common with the experiments, including the coexistence of different locomotion states for identical flapping motions. Figure 4a shows the complete characterization of the dynamics using S=f/F, where F=U/L is the traversal frequency. Here, the simulation results depend only on flapping Reynolds number Re f , which plays a role analogous to f in experiments (see Methods). Dynamical signatures of strong interactions again include bistability of states, a hysteresis loop, and avoided values of S, suggesting that these are general features of locomotor arrays. Similar to the experiments, these simulations display a downward jump at whole integer S, but unlike the experiments the upward jump seems to occur at a quarter rather than half integer, leading to a wider band of avoided phases. Nonetheless, the flow fields at these transition points reveal constructive and destructive characters, respectively. For the slow mode near Re f ≈60 and S≈2, the vorticity plot of Fig. 4c shows that the wing slaloms between vortices, contributing vorticity of the same sign to each. For the fast mode near Re f ≈20 and S≈1.25, Fig. 4d shows that the wing intercepts and destroys each oncoming vortex core. Interestingly, previous 2D experiments of flapping foils in unsteady flows have also shown vortex slaloming and interception modes40 (see Discussion section).

Figure 4: Simulations of interacting wings. (a–d) An infinite array of synchronized or temporally in-phase wings is simulated by a single airfoil driven to flap up and down, and allowed to swim freely left to right across a periodic domain. (a) Schooling number S for increasing (blue) and decreasing (red) flapping Reynolds number, Re f . A non-interacting wing (dashed curve) swims at a speed intermediate between the two schooling modes. (b) Input power normalized by that of an isolated wing. (c) Computed vorticity field for the slow mode (blue circle in a): The wing slaloms between vortices. (d) Fast mode (red circle in a): the wing intercepts each vortex core. (e,f) Schooling dynamics and power consumption for an array in which nearest neighbours flap temporally out-of-phase with one another. Full size image

Importantly, these simulations allow us to compare the swimming speed and power input for interacting bodies versus an isolated, non-interacting wing. The latter case is simulated in a longer domain to ensure that terminal speed is reached in the absence of self-interactions (see Methods section). To compare speed, we plot as a dashed curve in Fig. 4a the computed S for an isolated wing, which is found to lie between the curves for the slow and fast modes. Thus, with respect to speed, hydrodynamic interactions can be beneficial or detrimental. To compare energetic cost, we use the computed vertical force and prescribed vertical velocity to form the average input power at steady state, and Fig. 4b shows the power investment in the presence of interactions normalized relative to that of a non-interacting wing at the same Re f . Surprisingly, both the slow and fast modes save power for nearly all Re f , and the presence of interactions can lead to savings of >50%. The slow mode always saves more than the fast, and the fast mode displays a moderate energetic increase relative to an isolated wing for a narrow range of Re f near the upward transition point.

These simulations also allow us to explore how differences in temporal phases between swimmers in an array affect their dynamics and power consumed. We consider a linear array in which each wing flaps out-of-phase in time with its nearest neighbours (see Methods section for details). To compare these results to the previous case of a temporally in-phase array, we use a natural generalization of the schooling number based on the total phase shift due to both the temporal phase ϕ T and spatial phase contributions: S=ϕ T /2π+f/F. Temporally in-phase arrays have ϕ T =0 and out-of-phase arrays have ϕ T =π, and for all cases this definition preserves the property that integer values of S represent spatially in-phase trajectories while half-integer S indicate spatially out-of-phase trajectories. In Fig. 4e, we characterize the locomotion of out-of-phase arrays using this modified schooling number, and we find that all the salient features of in-phase arrays are also present in this arrangement. Interestingly, important differences between in-phase and out-of-phase arrays arise in the power consumed, as seen by comparison of Fig. 4b,f. In particular, the fast mode for out-of-phase arrays is associated with higher power consumption, and indeed it is typically higher than that of an isolated swimmer. The origin for this difference is unclear but may be related to the more erratic and intense flows observed for the case of temporally out-of-phase arrays.

Mathematical model

Our experiments and simulations motivate a minimal model that describes the collective dynamics of a linear array of swimmers. As shown in Fig. 5a, an infinite array of bodies flapping in synchrony and spaced by a distance L is represented by a single body that repeatedly traverses a domain of length L that is specified by periodic boundary conditions. In our conception, the body’s horizontal speed is perturbed because it encounters the wake produced in its previous pass through the domain. The perturbation strength depends on the traversal time t−t′, which is the time elapsed since the body was last at the same location: X(t)=X(t′)+L. Models of this type take the form of a delay differential equation for the swimming speed . Here, U 0 is the speed in the absence of interactions—that is, the speed of a single, isolated swimmer—and ΔU represents the perturbation due to wing–wake interactions. The effect of memory is explicitly incorporated through the time delay t′, which is not a constant but rather depends on the dynamical history.

Figure 5: Mathematical model. (a) An infinite linear array of synchronized swimmers is represented by a single particle undergoing repeated passes across a domain specified by periodic boundary conditions. (b) Schooling number for a model with parameter values s=1.5, P=1.5, ɛ=1, τ=1 (see text for details). Full size image

Here we consider a specific model of this type given by the equation:

The first term describes how the speed of an isolated swimmer increases with flapping frequency f, where s and p are parameters. This power law dependence of speed on frequency is consistent with our measurements for a single wing. The second term represents the perturbation to the speed, where ɛ is the wing–wake interaction strength. Importantly, the perturbation depends on the difference 2πf(t−t′) in the current phase in the flapping cycle and the phase when last at the same location. One might expect that the forcing is a periodic function of this phase difference, and the cosine form, in particular, is found to yield model solutions that closely correspond to the experimental data (see below). Finally, the dissipation of flows, and thus weakening of interactions for longer traversal times, is captured by the exponential term with a decay timescale of τ.

We then seek steady swimming solutions corresponding to , where F=1/(t−t′) is the frequency with which the body crosses the domain. Putting these relationships in the above dynamical equation and taking L=1, we obtain a nonlinear algebraic equation relating F and f:

To illustrate the structure of the solutions, which are solved numerically, we display in Fig. 5b the schooling number S=f/F for a model with order-one parameter values, as given in the caption of Fig. 5. The solution curve S(f) displays a fold that consists of upper and lower stable branches (solid curves) connected by an unstable branch (dotted curve). The non-interacting swimmer (dashed curve) serves as a point of comparison.

At low f, the wing progresses slowly, S is large, and the solution resembles the non-interacting case. The fold appears at higher f and corresponds to the coexistence of a slow mode (S≈1) and a fast mode (S≈0.5), and at yet higher f the slow mode disappears. The speed of an isolated wing lies intermediate between these modes, as found in simulations. The model thus reproduces remarkably well the key observations from the experiments and simulations, suggesting that relatively simple interaction laws that include memory underlie the complex hydrodynamics of swimmer arrays.