2.1. The vortex sheet model

The $n$ bound vortex sheets representing the flat plates have positions

(2.1) \begin{equation} \boldsymbol{x}_b^j(\alpha,t)=(x_b^j,y_b^j)(\alpha,t),\quad \alpha \in [0,\unicode{x03C0}],\quad j=1,\ldots,n, \end{equation}

and sheet strengths $\gamma ^j(\alpha,t)$ . The sheets are initially on the horizontal $x$ -axis, with position

(2.2) \begin{align} x_b^j(\alpha,0)&=s(\alpha) - (j-1)(d_0+1),\quad s(\alpha)=\cos (\alpha)/2, \end{align}

(2.3) \begin{align} y_b^j(\alpha,0)&=0\,. \end{align}

That is, each plate has unit length, with $\alpha =0$ at the leading edge (right), $\alpha =\unicode{x03C0}$ at the trailing edge (left), and $d_0$ the tail-to-tip distance between plates. The horizontal evolution of $x^j_b$ is determined by the displacement $X^j(t)$ of plate $j$ from its initial position, while the oscillatory vertical component $y^j_b(t)$ is prescribed:

(2.4a) \begin{align} x^j_b(\alpha,t)&=X^j(t) + x_b^j(\alpha,0), \quad X^j(0)=0 \end{align}

(2.4b) \begin{align} y^j_b(\alpha,t)&=A\sin (\unicode{x03C0} t). \end{align}

The plate position $X^j(t)$ evolves according to the induced thrust and drag forces acting on the plate (see § 2.2). The plate vortex sheet strength is determined numerically by enforcing the no-penetration boundary condition on the plate walls.

The free vortex sheets representing the wakes are parametrised by the circulation parameter $\Gamma$ , where $\Gamma _T^j$ is the total circulation around the $j$ th sheet:

(2.5) \begin{equation} \boldsymbol{x}_f^j(\Gamma,t),\quad \Gamma \in [0,\Gamma _T^j],\quad j=1,\ldots,n. \end{equation}

Together, the free and bound vortex sheets induce the fluid velocity at any point $\boldsymbol{x}_o$ :

(2.6) \begin{align} \notag\boldsymbol{u}(\boldsymbol{x}_o,t) &= \frac {1}{2\unicode{x03C0} }\sum _{j=1}^n \left [ \int _0^{\unicode{x03C0} } \frac { (\boldsymbol{x}_o-\boldsymbol{x}_b^j(\alpha,t))^{\perp } } {|\boldsymbol{x}_o-\boldsymbol{x}_b^j(\alpha,t)\vphantom{x_f}|^2}\,\gamma ^j(\alpha,t)\, s^{\prime}(\alpha)\,\mathrm{d}\alpha \right. \\ &\quad\left. +\int _0^{\Gamma _T^j} \frac { (\boldsymbol{x}_o-\boldsymbol{x}_f^j(\Gamma,t))^{\perp } } {|\boldsymbol{x}_o-\boldsymbol{x}_f^j(\Gamma,t)|^2+\delta ^2}\,\mathrm{d}\Gamma \right], \end{align}

where $(x,y)^{\perp }=(-y,x)$ . Here, the parameter $\delta$ regularises the flow for points $\boldsymbol{x}_o$ near or on the free sheet, using the vortex blob method (Krasny Reference Krasny1986). If $\boldsymbol{x}_o$ is a point on the plate, then the integrals in (2.6) are evaluated in the principal value sense, and yield the average of the velocity above and below the plate.

The free sheet evolves with the fluid velocity (2.6) as

(2.7) \begin{equation} \frac {\mathrm{d}\boldsymbol{x}_f^j}{\mathrm{d} t}(\Gamma,t) =\boldsymbol{u}\left (\boldsymbol{x}_f^j(\Gamma,t),t\right). \end{equation}

Its circulation $\Gamma$ is determined using the circulation shedding algorithm of Nitsche & Krasny (Reference Nitsche and Krasny1994). Vorticity is shed at time $t$ by releasing a particle from the plate’s trailing edge, with vertical velocity equal to the plate’s vertical velocity, and horizontal velocity equal to the average horizontal fluid velocity $u(\boldsymbol{x}^j(\unicode{x03C0},t),t)=\overline {U}_e^j$ at the trailing edge, as obtained from (2.6). The particle’s circulation is such that the Kutta condition is satisfied, with

(2.8a) \begin{align} \frac {\mathrm{d}\Gamma _T^j}{\mathrm{d} t} &= -\frac {1}{2}[(u^j_+)^2-(u^j_-)^2], \end{align}

(2.8b) \begin{align} \overline {U}_e^j&=\frac {1}{2}(u^j_++u^j_-), \end{align}

(2.8c) \begin{align} \gamma ^j(\unicode{x03C0},t)&=u_-^j - u_+^j, \end{align}

Here, $u^j_+$ and $u^j_-$ are the horizontal velocity components above and below the $j$ th plate at the trailing edge, obtained from $\gamma ^j(\unicode{x03C0},t)$ and $\overline {U}_e^j$ .

2.2. Thrust and drag

The plates experience a thrust force due to the suction at their leading edges, which arises because the fluid velocity is singular there. Specifically, the thrust is obtained from the leading edge bound vortex sheet strength. The sheet strength is unbounded at the leading edge, with $\gamma ^j(\alpha,t)\sim \Delta s^{-1/2}$ as $\Delta s\rightarrow 0$ , where $\Delta s=1/2-s(\alpha)$ is the distance from the leading edge. However, the desingularised vortex sheet strength $\gamma ^j(\alpha)\,s^{\prime}(\alpha)$ remains bounded. The horizontal thrust force acting on the plate is given by (Fang et al. Reference Fang, Ho, Ristroph and Shelley2017; Eldredge Reference Eldredge2019)

(2.9) \begin{equation} F^j_x=\frac {\unicode{x03C0} }{4}\widetilde {\gamma }_j^2, \quad \text{where } \widetilde {\gamma }_j=\lim \limits _{s\to 1/2}\gamma ^j(s)\,\sqrt {\frac {1}{4}-s^2}, \end{equation}

and, with our choice of $s(\alpha)=\cos \alpha /2$ ,

(2.10) \begin{equation} \widetilde {\gamma }_j=-\lim \limits _{\alpha \to 0}\gamma ^j(\alpha)\,s^{\prime}(\alpha). \end{equation}

The relation (2.9) holds for steady or unsteady plate motion. For steady motions, the result reduces to $\widetilde {\gamma }=-2V$ , where $V$ is the vertical velocity of the leading edge, which is used by Alben (Reference Alben2010). For completeness, both the unsteady and steady results are derived here, in Appendices A and B, respectively. Appendix B also shows a comparison between the flow field around a plate in a steady vertical flow, and the unsteady flow field around an impulsively started plate moving vertically upwards with constant velocity, the latter being computed using the method described in § 2.3. We detail how the flow field (figure 14) and circulation (figure 15) in the unsteady case approach that of the steady case in the long-time limit.

As discussed at the beginning of § 2, the drag force acting on the plate is modelled to be proportional to $(U^j(t))^{3/2}$ , where $U^j(t)$ is the plate velocity. The horizontal motion of the $n$ plates is thus governed by the system of ordinary differential equations

(2.11a) \begin{align} \frac {\mathrm{d} X^j}{\mathrm{d} t}&= U^j, \quad X^j(0)=0 \end{align}

(2.11b) \begin{align} M\frac {\mathrm{d} U^j}{\mathrm{d} t}&=\frac {\unicode{x03C0} }{4}\widetilde {\gamma }_j^2-C_{{d}} (U^j)^{3/2},\quad U^j(0)=U_0. \end{align}

The flow evolution is determined by the system of ordinary differential equations (ODEs) (2.7), (2.8) and (2.11), in addition to a Fredholm integral equation for $\widetilde {\gamma }$ . Its numerical solution is described next.

2.3. Numerical method

The plates (2.1) are each discretised by $N_b+1$ point vortices corresponding to $\alpha _k=k\unicode{x03C0} /N_b$ , $k=0,\ldots,N_b$ . The free vortex sheets are discretised by $N_f+1$ vortex blobs with total circulation $\Gamma _T^j$ , created by shedding one blob at each time step; thus $N_f$ increases with time. Time is discretised by $t_m=m\,\Delta t(t)$ , where typically $\Delta t$ is smaller at early times near the beginning of the motion, and increases linearly from $\Delta t_1$ to a final time step $\Delta t_2$ at a time $t_1$ , after which it remains constant. The small initial time steps are needed to resolve the large initial velocities about the trailing edge. The system of ODEs (2.7), (2.8) and (2.11) is solved using the fourth-order Runge–Kutta method. At each stage, the following steps are taken (Sheng et al. Reference Sheng, Ysasi, Kolomenskiy, Kanso, Nitsche, Schneider, Childress, Hosoi, Schultz and Wang2012).

1. (i) Set the vertical plate velocity $V(t)=A\unicode{x03C0} \cos (\unicode{x03C0} t)$ .

2. (ii) For each plate $j=1,\ldots,n$ , find the sheet strength $\gamma ^j_k=\gamma ^j(\alpha _k,t)$ such that the vertical fluid velocity in (2.6) evaluated at the midpoints $\alpha _k^m=(\alpha _{k-1}+\alpha _k)/2$ equals the prescribed vertical plate velocity,

   (2.12) \begin{equation} \boldsymbol{u}\left (\boldsymbol{x}^j_b(\alpha ^m_k,t),t\right)\boldsymbol{\cdot} \boldsymbol{n} = V(t),\quad k=1,\ldots,N_b, \end{equation}
   where $\boldsymbol{n}=(0,1)$ is normal to the plate. Here, the singular integrals in the velocity (2.6) are computed with what is in essence the alternating point vortex method of Shelley (Reference Shelley1992). In addition, enforce zero total circulation around the plate and its wake. The resulting $(N_b+1)\times (N_b+1)$ linear system for $\gamma ^j_k$ , $k=0, \ldots, N_b$ , is inverted using a precomputed LU decomposition.

3. (iii) Set $\gamma ^j(\unicode{x03C0},t)= \gamma ^j_{N_b}$ and $\overline {U}_e^j= u (\boldsymbol{x}_b^j (\alpha _{N_b}^m,t),t)$ , and obtain $u_j^+$ and $u_j^-$ using (2.8b ) and (2.8c ). Then shed a particle $\boldsymbol{x}^j_{{new}}$ from each trailing edge with velocity

   (2.13) \begin{equation} \frac {\mathrm{d} \boldsymbol{x}^j_{{new}}}{\mathrm{d} t}= (\overline {U}_e^j,V) \end{equation}
   and circulation such that (2.8a ) is satisfied. As explained by Nitsche & Krasny (Reference Nitsche and Krasny1994, § 2.4(iii)), only separated flow is considered when determining $u_j^+$ and $u_j^-$ .

4. (iv) Move all previously shed particles with velocity

   (2.14) \begin{equation} \frac {\mathrm{d} \boldsymbol{x}^j_{f,k}}{\mathrm{d} t}= \boldsymbol{u}\left (\boldsymbol{x}^j_{f,k},t\right),\quad k=0,\ldots,N^j_f, \end{equation}
   and add $\boldsymbol{x}^j_{f,N^j_f+1}=\boldsymbol{x}^j_{{new}}$ , then increase $N^j_f$ by 1. Here, the integral over the bound vortex sheets in (2.6) is near-singular when the target point is near a plate. The near-singular integrals are computed using the corrected trapezoid rule of Nitsche (Reference Nitsche2021).

5. (v) Determine the desingularised leading edge sheet strength $\widetilde {\gamma }_j$ using (2.10), and move the plates by (2.11).

Typically, one new particle is shed per time step. In some cases, a particle is shed every other time step or more for the sake of computational efficiency.

3.1. Single plate, $n=1$

Here, we consider one plate moving with a prescribed oscillatory heaving motion of amplitude $A$ , given an initial horizontal velocity $U_0$ . The plate moves forward by the thrust force (2.9), as vorticity is shed from its trailing edge. Figure 2 shows the resulting horizontal plate velocity for a range of amplitudes $A$ and initial velocities $U_0$ . The plate horizontal velocity oscillates with the same frequency as the heaving motion. For each case, the figure shows both the actual oscillatory velocity $U(t)$ , and the cycle-averaged velocity obtained by averaging $U(t)$ over a moving window of size equal to the oscillation period. The figure shows that in all cases, the cycle-averaged plate velocity approaches a steady-state velocity $U_{\infty }$ that depends solely on the amplitude $A$ and is independent of the initial velocity $U_0$ . The steady-state velocity $U_{\infty }$ increases as $A$ increases, with values $U_{\infty }=0.665$ , 1.85, 3.58 and 5.83 for $A=0.1$ , 0.2, 0.3 and 0.4, respectively. These values of $U_{\infty }$ are approximately consistent with those obtained in experiments on heaving wings in a water tank (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), which motivates our choice of the drag coefficient $C_{{d}}$ .

Figure 3(a) shows the shed vortex sheet behind one plate oscillating with amplitude $A=0.2$ , at time $t=25$ . The given initial velocity is that of the corresponding steady state, $U_0=1.85$ . With each upward (downward) motion of the plate, positive (negative) vorticity is shed from the trailing edge, which leads to a sequence of counter-rotating vortex pairs shed during each oscillation period. By taking the curl of (2.6), we obtain the regularised vorticity corresponding to the shed vortex sheet, which is shown in figure 4(a). Soon after being shed, the vorticity settles into a sequence of almost equally spaced vortices of alternating sign, with the positive (negative) vortex positioned slightly below (above) the centreline $y=0$ .

Figure 3. Plate and vortex sheet positions at $t=25$ , for the flapping amplitude $A=0.2$ and $n=1$ , 2, 3 and 4 plates (from top to bottom). All plates are initially equispaced with a distance near $d_0=4.2$ .

3.2. Several plates, $n=2,3,4$

Figure 3(b) shows the two vortex sheets shed behind two plates, both heaving with amplitude $A=0.2$ . The colour scheme is the same as that in figure 1, with the vortex sheet shed by the leader (follower) plate indicated in blue (red). There is a large amount of stretching far downstream of the plates as the two sheets interact. We plot the sheets using a linear interpolant of the shed point vortices if neighbouring points are sufficiently close to each other. Otherwise, after a sufficiently large amount of stretching, we simply plot the individual point vortex positions. The initial distance between the plates is $d_0=4.2$ , which is, as will be seen later, near an equilibrium position. The two sheets are seen to roll up into four vortices per period, a phenomenon that is also evident in the top panel of supplementary movie 1. This is seen more clearly in figure 4(b), which shows a sequence of two positive vortices followed by two negative vortices per period. This pattern repeats, although in an irregular fashion. Groups of four vortices are discernible, but their relative positions do not settle into a clear steady pattern.

Figure 3(c) shows the three vortex sheets shed behind three plates heaving with amplitude $A=0.2$ . The initial separation between all plates is as before, $d_0=4.2$ . A pattern seems to appear farther downstream of the plates. Most notably, near $x=20$ , the vortices have spread further from the midline $y=0$ . Figure 4(c) reveals more clearly a repeating pattern of groups of six vortices, i.e. three positive followed by three negative. Sufficiently far downstream, for $x\lt 25$ , the vorticity is clearly concentrated away from the midline, as it has spread significantly farther than in the $n=2$ case shown in figure 4(b). A simulation with $n=3$ plates is shown in supplementary movie 2.

Figures 3(d) and 4(d) show the vortex sheets and corresponding vorticity shed behind $n=4$ moving plates. It is difficult to detect a pattern in the downstream vortex sheet positions. The vorticity is expected to consist of groups of eight vortices, i.e. four positive followed by four negative, but the pattern is less clear and more random. Furthermore, the spread away from the midline is not as big as for $n=3$ . A simulation with $n=4$ plates is shown in supplementary movie 3.

Figure 4. Same as figure 3, but the regularised vorticity is plotted instead of the vortex sheet. The plates are indicated in black. Absolute vorticity values larger than 4 are represented by the darkest blue and red colours.

The effect of the trailing plates on the leader is seen most clearly in figure 3. In all four cases, the given initial velocity is $U_0=1.85$ , and the initial distances between the plates are $d_j(0)= 4.2$ . The shape of the sheet between the first and second plates is basically unchanged in all cases from the single-plate case. That is, the second plate has little influence on the wake in front of it. Similarly, the wake behind the second and third plates is unchanged from the $n=2$ case, so the third plate also has little influence on the wake in front of it. It is evident that as the number of plates behind the front plate increases, the front plate has travelled slightly further by $t=25$ . However, after an initial transient, we observe that $U_{\infty }$ remains essentially unchanged, with only a small increase of less than 3 % as $n$ increases. We attribute this increase in the leader’s velocity, albeit small, to the fact that the leader’s bound vorticity is modified by the downstream plates’ bound and shed vorticity, as is evident from step (ii) of the numerical method presented in § 2.3. The leader’s bound vorticity in turn affects its thrust, as given by (2.9).

3.3. Steady-state positions for $A=0.2$ , $n=2$

The distances $d_j(t)$ between the plates evolve in time. We are interested in the long-term behaviour of $d_j$ and the associated possible steady-state configurations. This subsection presents results for a pair of wings ( $n=2$ ) with flapping amplitude $A=0.2$ , with a range of initial separation distances $d_0$ . Since the steady-state cycle-averaged plate velocity $U_{\infty }$ is found to always be independent of $U_0$ , from here on, the initial horizontal velocity of each plate is taken to be the steady-state cycle-averaged velocity of a single plate flapping with the same amplitude. We observe that the two plates evolve towards a steady configuration in which the distance between them, after undergoing several oscillations, approaches a constant, but that this equilibrium distance $d_1^{\infty }$ depends on the initial distance $d_0$ . Figure 5 shows snapshots of simulations performed for three different values of $d_0$ , at a time when the separation distance $d_1$ between the plates has reached the equilibrium values $d_1^{\infty,1}$ , $d_1^{\infty,2}$ and $d_1^{\infty,3}$ . Supplementary movie 1 also shows examples of simulations performed for different values of $d_0$ , which illustrates how changing $d_0$ causes the pair of plates to attain different equilibria.

Figure 5. Snapshot at $t=25$ of three of the schooling modes obtained for a pair of plates ( $n=2$ ) with flapping amplitude $A = 0.2$ . The plates are initially located near the first, second and third equilibria, which, from top to bottom, correspond to the distances $d^{\infty,1}_1=4.18$ , $d^{\infty,2}_1=8.01$ and $d^{\infty,3}_1=11.82$ , respectively.

Note that in figure 5(a), the distance $d_1^{\infty,1}$ between the plates is approximately the distance travelled by the plate in one period of the heaving oscillation, as evident by the single upward and downward wave of the shed sheet. In figure 5(b), the equilibrium distance is approximately twice that of figure 5(a), and in figure 5(c) it is approximately three times. The distance travelled by a plate with velocity $U_{\infty }$ in one oscillation period is the wavelength $\lambda =U_{\infty }/f$ , where $f=1/2$ . Using the value $U_{\infty }=1.85$ corresponding to $A=0.2$ , we find that the equilibrium positions satisfy

(3.1) \begin{align} d^{\infty,1}_{1}\approx 4.18\approx 1.13\lambda,\quad d^{\infty,2}_{1}\approx 8.01 \approx 2.16\lambda \quad \text{and}\quad d^{\infty,3}_{1}\approx 11.82 \approx 3.19\lambda. \end{align}

Clearly, the plates have settled at equilibrium configurations near integer values of the non-dimensional schooling number (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015)

(3.2) \begin{equation} S_k=\frac {d^{\infty,k}_{1}}{\lambda }. \end{equation}

Note that while the first equilibrium has schooling number $S_{1}=1.13$ , i.e. slightly larger than 1, the following schooling numbers $S_2=2.16$ and $S_3=3.19$ increase approximately by integer values, $S_k\approx S_{k-1}+1$ .

3.4. Equilibrium schooling states and their loss of stability

Figure 6 shows the dependence of the distances $d_j(t)$ between plates on the initial separation distance $d_0$ , for in-line formations of $n=2$ , 3 and 4 wings (columns from left to right) and the three flapping amplitudes $A = 0.4$ , 0.3 and 0.1 (top to bottom rows). These flapping amplitudes complement the value $A=0.2$ considered in § 3.3. We proceed by discussing each column in turn.

The first column shows results for $n=2$ plates. It shows the distance $d_1(t)$ between the plates as a function of time, for 35 or 36 different initial separation distances $d_0$ . During an initial time interval of length approximately $d_0/U_{\infty }$ , the vorticity of the leader has not yet reached the follower, and both plates travel with the same velocity, each without the influence of the other. During that time interval, the distance between them stays constant. The larger the initial distance $d_0$ , the longer this initial time interval. Once the vortex wake of the leader reaches the follower, that vorticity changes the evolution of the follower, and the distance between the plates changes. In all cases, the distance $d_1(t)$ initially oscillates, then the amplitude of these oscillations decreases in time as $d_1$ approaches a constant steady-state value. Notice that the steady-state distances $d_1^{\infty,k}$ decrease in proportion to $U_{\infty }$ as $A$ decreases from 0.4 (top row) to 0.1 (bottom row). That is, as the heaving amplitude is decreased progressively, the plates are closer together in their steady configurations. Note also that the oscillations in $d_1$ are larger, relative to the value of $d_1$ , for smaller values of $A$ .

Figure 6. Time evolution of the distances $d_1(t)$ (left-hand column), $d_2(t)$ (middle column) and $d_3(t)$ (right-hand column) for in-line formations of $n=2$ , 3 and 4 plates, respectively. The plots in the top, middle and bottom rows correspond to the heaving amplitudes $A = 0.4$ , 0.3 and 0.1, respectively. The different curves in each plot are obtained by varying the initial distances $d_j(0)$ , and are colour-coded according to the equilibrium schooling mode reached for two plates ( $n=2$ , left-hand column). The schooling numbers $S_k=d^{\infty,k}_1/\lambda$ corresponding to each equilibrium distance $d^{\infty,k}_1$ are written in black in the first column, with the differences between them written in red. Here, $\lambda = 2U_{\infty }$ is computed using the steady-state $U_{\infty }$ for the pair of plates.

As already observed in figure 5, there are several steady states, depending on the value of $d_0$ . The steady states are measured by their schooling number (3.2), using the value for $U_{\infty }$ corresponding to the given value of $A$ . The distinct values of $d_0$ yield six distinct corresponding steady states, $S_1,\ldots,S_6$ . The schooling numbers associated with each steady state are listed in black in the first column of figure 6, and the differences between them are in red. For all three values of $A$ , $S_1$ is slightly larger than 1, with larger schooling numbers increasing by integer values as observed in § 3.3.

The curves are colour-coded according to the steady state reached. The curve with smallest $d_0$ is black and reaches zero in finite time, which corresponds to collision of the two plates. While generally the leader travels without being influenced much by the follower, the follower may accelerate and collide with the leader if the initial separation distance is sufficiently small.

The second column in figure 6 illustrates the case of $n=3$ plates, and shows the distance $d_2(t)$ between the second and third plates. As for $n=2$ , the addition of an $n$ th plate behind a group of $n-1$ plates leaves the motion of the first $n-1$ plates mostly unchanged. That is, for $n=3$ , the dynamics of $d_1(t)$ is roughly the same as that shown in the first column of figure 6, except in those cases of collision between the second and third plates, after which the computation stops. The curves are coloured by the equilibrium reached in the absence of the third plate, i.e. using the same colour for a given $d_0$ as in the first column of figure 6.

In the case of $n=3$ plates, the second and third plates move with the same velocity until the vortex wake from the first plate reaches the last plate, since, before then, they both are affected only by the (equal) wake from the plate directly ahead. The time interval of constant $d_2$ is therefore approximately $2d_0/U_{\infty }$ , twice as long as the time interval of constant $d_1$ . However, after this interval, the motion is more irregular than in the case of two plates, as is evident by comparing the first and second columns of figure 6. For $A=0.4$ (top row), one trajectory close to $S_2$ (red) jumps to $S_1$ . For $A=0.3$ (middle row), there are more trajectories that move to a different, typically lower, equilibrium schooling number. There are also more collisions between the plates. The situation is noticeably worse for $A=0.1$ (bottom row). Only a few curves reach their corresponding equilibrium, with the rest either exhibiting collisions or being on the path to collision. We thus conclude that the equilibria $S_k$ appear to lose stability when $n$ is increased from 2 to 3, and when $A$ is decreased from 0.4 to 0.1. We note that some of the curves showing $d_2$ , e.g. the black curves at the bottom of the panels in the second column, terminate abruptly despite staying strictly positive. This behaviour indicates that the first two plates have collided, while the third has not.

The case of $n=4$ plates, with $d_3$ plotted in the third column of figure 6, confirms this pattern. There are more irregular paths for $A=0.4$ (top row) than there were with $n=3$ , with one, in green, possibly approaching collision. The results for $A=0.3$ (middle row) show many curves that end at a finite time, corresponding to collisions in $d_2$ that are shown in the middle column. In addition, there are more collisions and many more curves that leave their original equilibrium. The basins of attraction of these equilibrium schooling modes appear to have shrunk. The results for $A=0.1$ are the most striking, as most of the curves either exhibit collisions or are on the path to collision, and only a few reach the final time $t=60$ (30 flapping periods).

Taken together, the simulation results shown in figure 6 demonstrate that both as the number $n$ of plates in an in-line formation increases, and as the heaving amplitude $A$ decreases, the basins of attraction of the equilibria shrink, with most of the stability being lost for $A=0.1$ and $n=4$ . The loss of stability can be attributed to the increasing relative size of the oscillations in $d_k$ observed as $n$ increases and $A$ decreases.

3.5. Stabilisation of schooling modes

To prevent the loss of stability and collisions observed in figure 6, we propose a simple stabilisation algorithm: if the distance decreases (increases) between a given plate and the one directly ahead of it, then slow it down (speed it up) by decreasing (increasing) its flapping amplitude. We implement this mechanism by allowing the heaving amplitude $A$ in (2.4b ) to depend on each plate, and evolving the amplitude $A^j$ of plate $j$ according to the rule

(3.3) \begin{equation} A^j(t)=A_0[1+\beta\, d_j^{\prime}(t)], \end{equation}

where $A_0$ is the initial amplitude, and $d_j^{\prime}=U^{j-1}-U^j$ is the rate of change of the distance between the plates. The updated values of $A^j$ are entered in step (i) of the numerical method (§ 2.3), where $V$ is replaced by $V_j(t)=A^j(t)\,\unicode{x03C0} \cos (\unicode{x03C0} t)$ , and thereby affect all remaining steps (ii)–(v) at each time step.

Figure 7. Same as figure 6, but with the stabilisation rule (3.3) and $n=4$ throughout. The stabilisation factors are $\beta =0.15$ , 0.4 and 5.0 for the flapping amplitudes $A=0.4$ , 0.3 and 0.1, respectively. Here, the curves are colour-coded simply by the equilibrium state that they reach.

The stabilisation parameter $\beta$ determines how fast the amplitude changes occur as the distance between plates changes. Small values of $\beta$ may not be sufficient to prevent collisions in finite time. Large values of $\beta$ lead to fast adjustment, but can also lead to stiffness of the numerical scheme. For all runs shown below, we used $0.15\leqslant \beta \leqslant 5$ .

Figure 7 shows the distances $d_1$ , $d_2$ and $d_3$ between each pair of plates for a formation of $n=4$ wings, computed with the stabilisation mechanism. The figure shows that the stabilisation either eliminates or significantly reduces the oscillations in $d_j$ , which quickly approach a steady state. Apart from the smallest value of $d_0$ (black curves), all of the collisions observed in figure 6 are avoided, with steady-state solutions approached even in the most unstable case, $n=4$ and $A=0.1$ . We note that $d_1$ is not shown with $n=2$ , nor $d_2$ with $n=3$ , since in the absence of collisions, the results are the same as with $n=4$ . For the sake of clarity, the curves in figure 7 are colour-coded by the equilibrium state that they reach, as opposed to the equilibrium state reached by the corresponding wing pair ( $n=2$ ) in figure 6. Supplementary movie 4 contains examples of simulations with $n=3$ plates, and shows how the stabilisation algorithm prevents a collision between the two trailing plates.

The stabilisation has a significant impact on the vortex wake behind the plates, as an equilibrium state is reached more quickly. In figure 8, we compare snapshots of the vortex sheet behind two plates without (figure 8 a) and with (figure 8 b) stabilisation. Figure 8(a) is the same as shown earlier in figure 3(b), but at a different scale. The formation of a wake of quadrupoles is seen early on (near the plates) in figure 8(a), but no clear pattern emerges far from the plates. On the other hand, the results obtained with stabilisation (figure 8 b) clearly show the emergence of a stable pattern of quadrupoles behind the plates. A positive (negative) vortex is created behind each plate by the plate’s upward (downward) motion. With stabilisation, these four vortices per oscillation period are positioned in a regular stable configuration, which grows slightly in time and space as the plates move forward. Figure 9 shows the vorticity in the wake, corresponding to the regularised vortex sheet shown in figure 8. Figure 9(b) clearly shows the stable pattern in the wake, growing in time and space. The contrast with figure 9(a), which shows much less organised wake vorticity without any growth, is striking.

Figure 8. Wake at $t=25$ behind a pair of plates ( $n=2$ ) with heaving amplitude $A_0=0.2$ and initial separation distance $d_0 = 4.2$ . The plots are (a) without stabilisation, and (b) stabilised according to (3.3) with $\beta =2$ .

Figure 9. Vorticity at $t=50$ , with the same parameters as in figure 8.

Figure 10. Vorticity at $t=50$ behind three plates ( $n=3$ ) with heaving amplitude $A_0=0.2$ and initial separation distance $d_0 = 4.2$ . The plots are (a) without stabilisation, and (b) stabilised according to (3.3) with $\beta =2$ .

Figure 10 shows the vorticity in the wake of $n=3$ plates, with and without stabilisation. Here, each heaving cycle of the three plates generates six vortices per period. Without stabilisation (figure 10 a), the vortices exhibit a coherent pattern only at short times after separation, but they become less organised further downstream. By contrast, with stabilisation (figure 10 b), a regular configuration of six vortices per period clearly emerges. The width of the wake grows in time and space, as was the case for a pair of plates ( $n=2$ , figure 9 b).

We were unable to compute such stable patterns of eight vortices per period in the wake of $n=4$ plates due to the computational cost of doing so. While the stabilisation mechanism allows stable schooling modes to be achieved relatively quickly, the vortex configurations take longer to reach equilibrium. Long time simulations using fast tree-algorithms (Greengard & Rokhlin Reference Greengard and Rokhlin1987; Wang, Krasny & Tlupova Reference Wang, Krasny and Tlupova2020) can be performed in the future to study this phenomenon further.

3.6. Reduced model based on thin-aerofoil theory

We now propose and analyse a reduced model to provide rationale for two of the phenomena observed in our simulations (§ 3.4), namely, that collectives of heaving plates in an in-line formation are increasingly stabilised as the flapping amplitude is increased progressively, and such collectives become less stable as the number of plates increases. To that end, we deploy the seminal thin-aerofoil theory of Wu (Reference Wu1961), who considered the motion of a single flat plate executing a small amplitude heaving motion ( $A \ll 1$ ) while translating horizontally at velocity $U$ in an inviscid and incompressible two-dimensional fluid. To simplify the problem, the plate is assumed to shed a flat vortex sheet from its trailing edge, and the sheet strength is determined by the Kutta condition. By solving the linearised Euler equations using conformal mapping, the (dimensionless) thrust $T_0$ on the plate is found to be

(3.4) \begin{align} T_0 = \frac {\unicode{x03C0} }{2}(\unicode{x03C0} A)^2\,|C(\sigma)|^2,\end{align}

where

(3.5) \begin{align} C(\sigma)=\frac {\mathrm{K}_1(\mathrm{i}\sigma)}{\mathrm{K}_0(\mathrm{i}\sigma)+\mathrm{K}_1(\mathrm{i}\sigma)} \end{align}

is the so-called Theodorsen function, $\mathrm{K}_{0,1}$ are the modified Bessel functions of the second kind of orders zero and one, and $\sigma = \unicode{x03C0} /(2U)$ is the reduced frequency. An analogous theory has been developed to describe how a follower plate located downstream of the leader and executing the same heaving kinematics would be affected by the leader’s wake (Wu & Chwang Reference Wu and Chwang1975; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016). In this theory, the no-penetration boundary condition is not satisfied simultaneously on both wings; instead, an approximation for the hydrodynamic thrust on the follower is found by assuming that the leader’s wake is unaffected by the follower’s presence. This assumption is consistent with our numerical simulations before the leader’s wake has reached the follower, as discussed in § 3.2. This thrust is found to be $T = T_0 [1+F(d)]$ , where $d$ is the distance between the leader’s trailing edge and the follower’s leading edge (as in figure 1), and

(3.6) \begin{align} F(d) &= G(\sigma)^2+2G(\sigma)\cos \left [2\unicode{x03C0} \left (\frac {d}{2U}-g(\sigma)+\frac {\sigma }{\unicode{x03C0} }\right)\right]. \end{align}

Here, the real-valued functions $G(\sigma)$ and $g(\sigma)$ are defined through the relation

(3.7) \begin{align} G(\sigma)\,\mathrm{e}^{2\unicode{x03C0} \mathrm{i}\, g(\sigma)}&= \frac {{K}_1(-\mathrm{i}\sigma)\,{K}_0(\mathrm{i}\sigma)+{K}_0(-\mathrm{i}\sigma)\,{K}_1(\mathrm{i}\sigma)}{{K}_1(\mathrm{i}\sigma)\,\left [{K}_0(\mathrm{i}\sigma)+{K}_1(\mathrm{i}\sigma)\right]}. \end{align}

We thus propose the following model for a collective of $n$ plates that execute small-amplitude heaving kinematics in an in-line formation:

(3.8) \begin{align} M\ddot {x}_1+C_{{d}}\dot {x}_1^{3/2}&=T_0,\nonumber \\ M\ddot {x}_m+C_{{d}}\dot {x}_m^{3/2}&=T_0\left [1+F(x_{m-1}-x_{m}-1)\right], \quad 2\leqslant m\leqslant n, \end{align}

where $x_m$ is the horizontal position of the centre of the $m$ th plate. We have assumed that each plate is influenced only by its nearest upstream neighbour; that is, the leader wing ( $x_1$ ) is unaffected by all of the followers, and moves at a steady speed $U$ determined by the balance between drag and thrust forces on it:

(3.9) \begin{align} C_{{d}}U^{3/2}&=\frac {\unicode{x03C0} }{2}\left (\unicode{x03C0} A\right)^2\,|C(\sigma)|^2. \end{align}

Every other wing is influenced by its upstream neighbour through the function $F$ , which depends on the steady speed $U$ . That is, our theory is valid in the quasi-steady regime in which the wings’ velocities are all close to $U$ . The nearest-neighbour assumption may be justified using our simulations that employ stabilisation (figure 7); specifically, $d_1^{\infty,k}\approx d_2^{\infty,k}\approx d_3^{\infty,k}$ , which suggests that the equilibrium distances are mostly insensitive to the number of upstream wings. The assumption is also supported by recent experiments on up to five flapping hydrofoils in an in-line formation (Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024). Table 1 shows a comparison between the steady speeds $U$ obtained in our vortex sheet simulations and those predicted from (3.9), with surprisingly good agreement given the simplified nature of our reduced model.

Table 1. The values of the velocity $U$ of a single wing, as obtained in our vortex sheet simulations (§ 3.1) and predicted by thin-aerofoil theory from (3.9), are shown for different values of the flapping amplitude $A$ .

We proceed by seeking steadily translating solutions of (3.8), in which the wings move at a constant speed $U$ with a constant separation distance $d$ between them, $x_{m}(t)-x_{m+1}(t)=d+1$ . By substituting the form $x_m(t)=Ut-(m-1)(d+1)$ into (3.8), we find that $d$ satisfies the algebraic equation

(3.10) \begin{align} F(d)=0. \end{align}

This equation has infinitely many positive solutions $d^{\infty,k}$ for $k\in \mathbb{Z}_{\gt 0}$ . For the parameter values adopted in this paper, $M = 1$ , $C_{{d}}=0.1$ and $A = 0.1{-}0.4$ , we find that the corresponding schooling numbers $S_k=d^{\infty,k}/(2U)$ are equal to $S_1+k-1$ , where $S_1$ is between 0.95 and 1.1, in good agreement with values obtained in our vortex sheet simulations (figure 6). We note that a similar result was obtained for two plates by Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016) using a similar thin-aerofoil theory, and by Baddoo et al. (Reference Baddoo, Moore, Oza and Crowdy2023) using a more sophisticated thin-aerofoil theory based on conformal maps of multiply connected domains. Reduced models of flapping wings that approximate the hydrodynamic interactions using time-delay differential equations also recover near-integer values of the equilibrium schooling numbers (Heydari, Hang & Kanso Reference Heydari, Hang and Kanso2024; Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024). We also note that (3.10) is satisfied by a second family of solutions with schooling numbers $\tilde {S}_k=\tilde {S}_1+k-1$ , where $\tilde {S}_1$ is between 0.6 and 0.7 for the parameter values adopted in this paper. We neglect these solutions in the forthcoming discussion because they are linearly unstable to perturbations and thus not observed in our simulations or in laboratory experiments (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016).

We now examine the linear stability of the state in which the distance between each pair of wings corresponds to the schooling number $S_1$ , the analysis of the other schooling numbers $S_k$ being identical due to the periodicity of the function $F$ . To that end, we substitute the form $x_m(t)=Ut-(m-1) (d^{\infty,1}+1)+\epsilon\, \tilde {x}_m(t)$ into (3.8), and retain terms at leading order in $\epsilon$ . The linearised equations are thus (dropping the tildes)

(3.11) \begin{align} M\ddot {x}_1+\frac {3C_{{d}}}{2}U^{1/2}\dot {x}_1&=0,\nonumber \\ M\ddot {x}_m+\frac {3C_{{d}}}{2}U^{1/2}\dot {x}_m&=T_0\xi \left (x_{m-1}-x_{m}\right),\quad 2\leqslant m\leqslant n, \end{align}

where $\xi = F^{\prime } (d^{\infty,1})$ . We assume that an impulsive perturbation is applied to the leader:

(3.12) \begin{align} x_1(0)=1,\quad \dot {x}_1(0)=s_1\equiv -\frac {3C_{{d}}U^{1/2}}{2M}\quad \mathrm{and} \quad x_m(0)=\dot {x}_m(0)=0\ \text{for}\ m\geqslant 2. \end{align}

It is straightforward to obtain the solutions

(3.13) \begin{align} x_1(t)&={\rm e}^{s_1t}\quad\mathrm{and}\quad x_2(t)={\rm e}^{s_1t} + \mathrm{Re}\big[B_0^{(2)}{\rm e}^{s_2t}\big],\quad \mathrm{where}\quad B_0^{(2)}=-1+\frac{{\rm i} s_1}{2\omega},\nonumber \\\omega &= \frac{|s_1|}{2}\sqrt{\frac{16MU^{1/2}\xi}{9C_{\mathrm{d}}}-1}\quad\mathrm{and}\quad s_2=\frac{s_1}{2}+{\rm i}\omega \end{align}

is a root of the quadratic polynomial $p(s)\equiv Ms^2+(3C_{{d}}/2)U^{1/2}s+T_0\xi$ . We find that $\omega \in \mathbb{R}$ for the parameter values adopted in this paper. Each wing after the second oscillates in resonance with its upstream neighbour, so the general solution to (3.11) has the form

(3.14) \begin{align} x_m(t)=\mathrm{e}^{s_1t}+\mathrm{Re}\left [\sum _{k=0}^{m-2}B_k^{(m)}t^k\,\mathrm{e}^{s_2t}\right]. \end{align}

The coefficients $B_0^{(m)}$ may be obtained using the initial conditions:

(3.15) \begin{align} B_0^{(m)}=-1+\frac {\mathrm{i}}{\omega }\left (\frac {s_1}{2}+\mathrm{Re}\left [B_1^{(m)}\right]\right),\quad m\geqslant 3. \end{align}

To obtain the remaining coefficients, we write (3.11) in the compact form $\mathcal{L}x_m=T_0\xi x_{m-1}$ for $m\geqslant 2$ , where $\mathcal{L}=p(\mathrm{d}/\mathrm{d} t)$ . Using the facts that

(3.16) \begin{align} \mathcal{L}\,\mathrm{e}^{s_1t}&=T_0\xi\, \mathrm{e}^{s_1t},\quad \mathcal{L}\,\mathrm{e}^{s_2t}=0,\quad \mathcal{L}t\,\mathrm{e}^{s_2t}=2\mathrm{i}\omega M\,\mathrm{e}^{s_2t}\nonumber \\ \mathrm{and}\quad \mathcal{L}t^k\,\mathrm{e}^{s_2t}&=\left [2\mathrm{i}\omega Mkt^{k-1}+Mk(k-1)t^{k-2}\right]\mathrm{e}^{s_2t}\quad \mathrm{for}\ k\geqslant 2, \end{align}

it can be shown that the coefficients $B_k^{(m)}$ satisfy the recursion relation

(3.17a) \begin{align} &B_{m-2}^{(m)}=-\frac {\mathrm{i} T_0\xi B_{m-3}^{(m-1)}}{2\omega M(m-2)}\quad \mathrm{for}\ m\geqslant 4, \end{align}

(3.17b) \begin{align} &2\mathrm{i}\omega MB_{k+1}^{(m)}(k+1)+MB_{k+2}^{(m)} (k+1)(k+2)=T_0\xi B_k^{(m-1)}\quad \mathrm{for}\ k=0,\ldots,m-4. \end{align}

Taken together, (3.15) and (3.17) determine all of the coefficients in the solution (3.14).

We proceed by making two observations about this solution. First, the equilibrium schooling states are linearly stable because the perturbations decay in time, $x_m(t)\rightarrow 0$ as $t\rightarrow \infty$ . Moreover, the decay rate $s_1/2$ increases with $U$ , as is evident from (3.12), and thus with the flapping amplitude $A$ , providing a rationale for the observation from our vortex sheet simulations (§ 3.4) that a larger flapping amplitude stabilises the schooling state. Second, by plotting the solution (3.14) for the parameter values adopted in this paper, we observe that the perturbations are amplified downstream, in that at intermediate times they exhibit oscillations that grow in time and have progressively larger amplitude as $m$ increases. To see why this is the case, we observe that the exponential decay factors $\mathrm{e}^{s_2t}$ in (3.14) are multiplied by increasingly large powers of $t$ as $m$ increases. More precisely, note that (3.17a ) can be solved explicitly for $B_{m-2}^{(m)}$ , from which the long-time asymptotic behaviour of $x_m(t)$ can be deduced:

(3.18) \begin{align} x_m(t)\sim \text{Re}[X_m(t)]\ \text{as}\ t\rightarrow \infty,\quad \text{where}\ X_m(t)=\frac {B_0^{(2)}}{(m-2)!}\left (\frac {-\mathrm{i} T_0\xi }{2\omega M}\right)^{m-2}t^{m-2}\,\mathrm{e}^{s_2t}. \end{align}

We find that for the parameter values adopted in this paper, $|X_m(t)|$ provides an adequate approximation for the envelope of the oscillating solution $x_m(t)$ in (3.14). Moreover, it is easy to see that $|X_m(t)|\propto t^{m-2}\,\mathrm{e}^{s_1t/2}$ attains its maximum value

(3.19) \begin{align} X_m^*\equiv |X_m(t_m)|=\frac {\left |B_0^{(2)}\right |}{(m-2)!}\left (\frac {T_0\xi (m-2)}{\omega M|s_1|\mathrm{e}}\right)^{m-2}\quad \text{at the time}\ t=t_m\equiv \frac {2(m-2)}{|s_1|}. \end{align}

By using Stirling’s formula, $m!\sim \sqrt {2\unicode{x03C0} m}\,(m/\mathrm{e})^m$ as $m\rightarrow \infty$ , we conclude that the peak of the envelope $X_m^*$ increases with $m$ for sufficiently large $m$ provided that $T_0\xi /(\omega M\,|s_1|) \gt 1$ , as is the case for the parameter values adopted in this paper. Moreover, the time $t_m$ at which the peak is attained increases with $m$ , which further illustrates that the oscillations propagate downstream. The linear theory thus provides a qualitative rationale for the observation from our vortex sheet simulations (§ 3.4) that the schooling states become less stable as the number of wings increases.

Figure 11. Results of numerical simulations of the (nonlinear) reduced model (3.8). The initial velocity of the leader is $1.3U$ , and the remaining plates are initialised in an equilibrium schooling mode with velocity $U$ (given by (3.9)) and inter-plate distance $d^{\infty,1}$ (given by (3.10)). The middle panel shows the inter-plate distances for $n=3$ wings flapping with amplitude $A=0.3$ ; after a transient, the system returns to the equilibrium schooling state. The right-hand panel shows the effect of adding a plate ( $A = 0.3$ , $n=4$ ), and the left-hand panel shows the effect of decreasing the flapping amplitude ( $A=0.2$ , $n=3$ ). It is evident that both lead to a collision between the trailing plates.

While the foregoing linear theory predicts that small-amplitude perturbations to the equilibrium schooling state eventually decay in time, we find that transient linear growth combined with nonlinear effects can cause finite-amplitude perturbations to dislodge the trailing plates from their equilibrium positions. Specifically, we performed numerical simulations of the nonlinear reduced model (3.8), initialising the plates with the positions $x_m(0)=-(m-1) (d^{\infty,1}+1)$ and velocities $\dot {x}_m(0)=U$ for $m \gt 1$ and $\dot {x}_1(0)=1.3U$ ; that is, we assume that the trailing wings initially move with the steady speed of a single wing, as determined by (3.9), and subject the leader to a finite-amplitude velocity perturbation. The results are shown in figure 11. The middle panel shows the inter-plate distances $d_i$ for $n = 3$ wings and flapping amplitude $A = 0.3$ : after a transient, the trailing wings eventually settle back into their equilibrium positions. By contrast, increasing the number of wings to $n = 4$ (right-hand panel) or decreasing the flapping amplitude to $A = 0.2$ (left-hand panel) causes the trailing wings to become dislodged from their equilibrium positions and collide. A qualitatively similar behaviour is observed for the other flapping amplitudes considered in this paper. We thus conclude that despite its simplifications, the reduced model presented in this section captures the two salient features of our vortex sheet simulations, namely, that schooling modes become progressively less stable as the number of wings is increased or the flapping amplitude is decreased.

We note that a number of physical effects are neglected in our reduced model (3.8) for the sake of analytical simplicity. Specifically, the assumptions are that the wake shed by each plate is flat and unaffected by modulations in the plate’s horizontal speed, both of which are at odds with our vortex sheet simulations. The model also assumes nearest-upstream-neighbour interactions; while our simulation results seem to confirm this assumption (§ 3.2), strictly speaking each plate feels the vortex sheets shed by all of the other plates.