3.1. Schooling state

When a lateral forced heaving motion is applied to the leading edge of the flexible plate, the plate begins to move in the horizontal direction. The cruising speed of the $i$ th plate is defined as its time-averaged forward speed at equilibrium, i.e. $U_{c,i}= {1}/{T}\int _{t'}^{t'+T} u_i(t) {\rm d}t$ , where $t'$ is the time after reaching equilibrium, and $u_i(t)$ denotes the instantaneous speed of plate $i$ . The leader’s cruising speed $U_{c,1}$ is influenced by its self-propulsion kinematics, the channel width $H$ and the hydrodynamic interaction with the following plate (specifically, whether they are in a compact configuration). The dynamical states of the multi-plate systems are characterised by the followers’ relative motion to the leader, which can be quantified by the schooling number (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), $\tilde {G}_{i,i+1}=G_{i,i+1}/\lambda$ , where $\lambda =U_{c,1}/f$ represents the wavelength of the leader’s swimming trajectory.

Our simulations show three distinct schooling states: hybrid mode, regular mode and unstable mode, achieved by varying $N$ and $H$ . Figure 1( $b$ ) shows snapshots of the configurations for the stable hybrid and regular modes. In the stable ordered configuration at equilibrium, the followers possess the same cruising speed as the leading plate, with the whole group moving at a speed of $U_c=U_{c,i}$ , $i=1,2, \ldots , N$ . For the hybrid mode, as the group size $N$ increases, the plates following the leader may spontaneously split into several subgroups or individual plates. This mode typically occurs when the distance $H$ is small. For the regular mode with a sparse configuration, the gap spacings between any two neighbouring plates are nearly identical, approximating one wavelength, and no subgroups emerge. When the channel width is sufficiently large, the wake vortices become increasingly chaotic as the group size $N$ increases (notably when $N \gt 8$ ). This chaotic wake destabilises the alignment of trailing plates, causing the separation distance between them to progressively decrease until collisions occur (Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024; Nitsche, Oza & Siegel Reference Nitsche, Oza and Siegel2025). We classify this state as an unstable mode. In contrast, for sufficiently small channel widths and large group sizes (e.g. $H\leqslant 2.0$ and $N=72$ ), a divergent unstable mode emerges, resulting in irreversible gap widening between specific neighbouring plates (Heydari, Hang & Kanso Reference Heydari, Hang and Kanso2024).

Figure 2.

( $a$ , $b$ ) Normalised longitudinal distance $\tilde {D}_{1i}=D_{1i}/\lambda$ $(i=2,\ldots ,N)$ between the following plates and the leader and ( $c$ , $d$ ) the schooling number $\tilde {G}_{i-1,i}=G_{i-1,i}/\lambda$ as a function of time. ( $a$ , $c$ ) Regular mode ( $N=36, H=3.5$ ). ( $b$ , $d$ ) Hybrid mode ( $N=36, H=2.0$ ). Colour gradient indicates plate index $i$ : blue ( $i=2$ ) to green ( $i=N$ ).

Two typical cases are analysed to quantify the evolution of ordered configurations: a regular-mode case ( $N = 36$ , $H = 3.5$ ) and a hybrid-mode case ( $N = 36$ , $H = 2.0$ ). For the two cases considered, figure 2( $a$ , $b$ ) shows the time-dependent locations of the following plates in the group, where the normalised gap is $\tilde {D}_{1i}=D_{1i}/\lambda$ ( $i=2,\ldots ,N$ ). In the regular mode at $H=3.5$ , the plates maintain an ordered formation while the group expands uniformly at a constant speed. The settling time for $\tilde {D}_{1i}$ to reach equilibrium scales linearly with its index $i$ . The formation fully stabilises at $t/T\gt 125$ , when all distances $\tilde {D}_{1i}$ reach steady state. In the hybrid mode ( $H=2.5$ ), the dynamics of $D_{1i}(t)$ differs significantly from that in the regular mode, exhibiting irregular oscillations. The system takes $t/T\gt 300$ to reach steady state.

Figures 2( $c$ ) and 2( $d$ ) illustrate the schooling number $\tilde {G}_{i-1,i}=G_{i-1,i}/\lambda$ plotted against time for the regular mode ( $N=36$ , $H=3.5$ ) and the hybrid mode ( $N=36$ , $H=2.0$ ), respectively. In the regular mode at equilibrium state, the gap distance between the first two plates is the largest among all adjacent plates ( $G_{12}\approx 1.13\lambda$ ), slightly exceeding one wavelength. As the plate number $i$ increases, $G_{i-1,i}$ decreases and converges to a unique curve (see the inset in figure 2 $c$ ), with minor fluctuations around its equilibrium value. In contrast, the hybrid mode exhibits a wide range of equilibrium gaps, including vanishing distances ( $G_{i-1,i}\approx 0$ ) that imply subgroup formation.

We systematically study the emergence of different schooling modes as functions of group size ( $N$ ) and channel width ( $H$ ). For specific $(H,\ N)$ pairs, bistability between hybrid and regular modes may emerge depending on initial conditions such as the initial gap spacing $G_0$ and flapping phase difference (Peng et al. Reference Peng, Huang and Lu2018c ; Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024). In the present configuration, bistability is precluded by enforcing uniform initial spacing ( $G_0 = 3$ ) and synchronised flapping kinematics, ensuring a uniquely determined system state. In the $N$ – $H$ plane, as shown in figure 3, the regular mode, hybrid mode and unstable state are represented by circles, diamonds and triangles, respectively. It is observed that the channel confinement can enhance group cohesion. When the group size is less than $9$ , the schooling configuration remains stable, exhibiting a hybrid mode for $H\lt 2.5$ (except for $N=3$ ) and a regular mode for $H\geqslant 2.5$ . Notably, the hybrid mode can be recovered for the case of $H = 1.5$ , $N = 3$ when a smaller $G_0$ (e.g. $G_0=1$ ) is employed. For $N\geqslant 9$ , an unstable region is observed above the regular-mode region, and as $N$ increases further, the boundary between the stable and unstable regions gradually shifts towards smaller values of $H$ , approaching $H=4$ . This indicates that as the number of individuals in the schooling group increases, a stronger channel confinement is required to maintain the group’s cohesion. When the channel confinement is further intensified ( $H\leqslant 2.0$ ), the stable ordered configurations exhibit a hybrid mode for $3\lt N\leqslant 36$ , where subgroups form within the schooling group with inter-individual gaps within subgroups approaching zero. However, for $N=72$ , the collective state transitions back to an unstable mode characterised by excessive spacing between certain neighbouring individuals, where the gap exceeds $5$ wavelengths.

3.2. Schooling performance

To quantitatively evaluate how group size ( $N$ ) and channel width ( $H$ ) affect stable formations, we analysed the normalised equilibrium gap spacing $\bar {G}_{i,i+1}/\lambda$ across the flexible plates. Here, $\bar {G}_{i,i+1}$ is defined as $ {1}/{T}\int _{t'}^{t'+T}G_{i,i+1} {\rm d}t$ . Figure 4( $a$ ) displays the normalised equilibrium gap distance $\bar {G}_{i,i+1}/\lambda$ versus gap index $i$ in the regular mode ( $H=3.5$ ) for group sizes $N\in \{4,9,18,27,36,72\}$ . All profiles with $N\in \{4,9,18,27,36\}$ collapse onto the $N=72$ profile, indicating that upstream formations are hydrodynamically insulated from group size variations. Near the leader, the spacing exhibits $i$ -dependent variations, while the downstream spacing asymptotically approaches $\bar {G}_{i,i+1}/\lambda = 0.89$ for $i\geqslant 18$ . Figure 4( $b$ ) shows $\bar {G}_{i,i+1}/\lambda$ versus $i$ for the hybrid mode ( $H=2.0$ ) with varying group sizes $N\in \{4,9,18,27,36\}$ . The results demonstrate that the spacing between the $i$ th and its following plate, $\bar {G}_{i,i+1}$ , is influenced by the group size $N$ , suggesting upstream propagation of disturbances induced by followers. Furthermore, figure 4( $c$ ) illustrates the normalised equilibrium gap distance $\bar {G}_{i,i+1}/\lambda$ as a function of the gap index $i$ for channel widths $H$ ranging from $1.5$ to $4.0$ , with the group size fixed at $N=36$ . For smaller $H$ (e.g. $H=1.5$ ), the spacing $\bar {G}_{i,i+1}/\lambda$ exhibits strong dependence on the gap index $i$ , characteristic of the hybrid mode. As $H$ increases beyond $2.5$ , the $\bar {G}_{i,i+1}/\lambda$ vs. $i$ curves gradually flatten into a nearly straight line, marking the transition to a regular mode with adjacent plates uniformly spaced at roughly one wavelength.

Figure 3.

Parameter space for assessing the stability of ordered configuration in multi-plate schooling, with respect to varying the group size ( $N$ ) and the channel width ( $H$ ). Symbols indicate: circles – regular mode; diamonds – hybrid mode; triangles – unstable mode.

Figure 4.

Normalised equilibrium gap distance $\bar {G}_{i,i+1}/\lambda$ as a function of gap index $i$ ( $i=2,3,..,N-1$ ) for three configurations: (a) regular mode at $H=3.5$ ; (b) hybrid mode at $H=2.0$ ; (c) system with fixed group size, $N=36$ .

Figure 5 shows the performance metrics(cruising speed $U_{c,i}$ , input power $P_i$ and swimming efficiency $\eta _i$ ) and their normalised forms ( $U_{c,i}/U_{c,1}$ , $P_i/P_1$ and $\eta _i/\eta _1$ ) for the regular mode ( $N=36$ , $H=2.0$ ) and hybrid mode ( $N=36$ and $H=3.5$ ), respectively. Here, $U_{c,1}$ , $P_1$ and $\eta _1$ represent the cruising speed, input power and swimming efficiency of an isolated plate in a channel of identical width. Notably, with the initial spacing set at $G_0=3$ , neither the regular nor hybrid mode induces subgroup formation between the leading and second plates. As a result, the leading plate’s swimming dynamics remains unaffected by trailing plates, demonstrating performance identical to that of an isolated plate under equivalent channel conditions. The input power of the $i$ th plate is calculated by $P_i=- {1}/{T}\int _{t'}^{t'+T}\int _0^1 {\boldsymbol{F}}_{s,i}(s,t)\boldsymbol{\cdot }\partial _t {\boldsymbol{X}}_i(s,t) {\rm d}s {\rm d}t$ , where ${\boldsymbol{F}}_{s,i}$ and ${\boldsymbol{X}}_i$ are the Lagrangian force and position vector on the $i$ th plate. The swimming efficiency $\eta _i$ is defined as the ratio of its kinetic energy to input work, i.e. $\eta _i=0.5MU_{c,i}^2/W_i$ , where $W_i=P_i T$ . It is observed that, at equilibrium, the cruising speed of each plate is uniform (figure 5 $a{,}b$ ) and matches the speed of an isolated plate. Figure 5( $c$ , $e$ ) shows that, in the regular mode, for small $i$ , the input power $P_i$ and swimming efficiency $\eta _i$ vary irregularly with $i$ because each individual experiences a distinct vortex environment. For larger $i$ , however, the energy consumption decreases monotonically with $i$ , while the swimming efficiency follows a nearly monotonic increasing trend, eventually converging to a constant value. Notably, the efficiency of downstream plates can reach $1.9$ times that of the leading plate. The hybrid mode, featuring subgroups and spatially aperiodic vortex flow, exhibits differing input powers ( $P_i$ ) and swimming efficiencies ( $\eta _i$ ) among individuals within the system (figure 5 $d{,}f$ ), despite their uniform swimming speeds. This variation in efficiency arises from their distinct fluid-mediated interactions.

Figure 5.

Schooling performance of the groups. ( $a{,}b$ ) Cruising speed $U_{c,i}$ and its normalised form $U_{c,i}/U_{c,1}$ ; ( $c{,}d$ ) input power $P_i$ and its normalised form $P_i/P_1$ ; ( $e,f$ ) efficiency $\eta _{i}$ and its normalised form $\eta _i/\eta _1$ . Results are shown for the ( $a{,}c{,}e$ ) regular mode ( $N=36$ , $H=3.5$ ) and ( $b{,}d{,}f$ ) hybrid mode ( $N=36$ , $H=2.0$ ). Here, $U_{c,1}$ , $P_1$ and $\eta _1$ correspond to the cruising speed, input power and efficiency of an isolated plate swimming in a channel of the same width.

Figure 6.

( $a$ ) Swimming speed $U_c$ , ( $b$ ) group efficiency $\eta$ and ( $c$ ) the relative efficiency difference $\eta _r = (\eta -\eta _1)/\eta _1$ between the group and an isolated plate in channels of identical width $H$ , plotted versus channel width for various group sizes $N=1,3,6,9,18,36$ and $72$ .

For a stable configuration, the swimming efficiency of the whole group is defined as $\eta = ({1}/{2})M\sum _{i=1}^N U_{c,i}^2/\sum _{i=1}^N W_i$ . The relative difference between the group’s swimming efficiency ( $\eta$ ) and that of an isolated plate in a channel of identical width is given by $\eta _r=(\eta -\eta _1)/\eta _1$ . Figure 6 shows the swimming speed $U_c$ , the group efficiency $\eta$ and the relative efficiency difference $\eta _r$ as functions of the channel width $H$ for varying group sizes $N$ . Figure 6( $a$ ) reveals that the cruising speeds of the groups match that of an isolated plate in channels of identical width. This occurs because the wake vortices generated by the leading plate increase inter-plate separation during motion, leading to a sparse configuration between the leading and trailing plates. In this regime, the leading plate’s motion becomes effectively independent of the trailing plates, mimicking the behaviour of an isolated plate in a channel of the same width. For stable configurations, this results in $U_{c}=U_{c,1}$ . Furthermore, the cruising speed decreases monotonically with increasing channel width $H$ (Jeong et al. Reference Jeong, Lee and Park2021), asymptotically approaching the free-swimming speed as the channel confinement diminishes.

Figure 6( $b$ ) demonstrates that the swimming efficiency $\eta$ of an isolated individual decreases monotonically with increasing channel width $H$ (Jeong et al. Reference Jeong, Lee and Park2021). The schooling groups exhibit generally higher swimming efficiency compared with isolated plates, except in certain cases with small group sizes $N$ (figure 6 $b$ , $c$ ). Notably, large schooling groups achieve significant energetic advantages when $H \in [3, 3.5]$ , which correlate with the formation of regularised vortex flow patterns (figure 9). These results suggest that well-organised vortex structures in the flow field not only facilitate schooling cohesion but also substantially enhance swimming performance.

3.3. Flow field

Figure 7.

The instantaneous vorticity contour for the case with $N=36$ and $H=3.5$ , with plates depicted by red curves. The flow field of the $36$ -plate system is segmented into four parts, represented in panel ( $a$ ) by four rows of small graphs. The x and y axes represent the position vector components. Panels ( $b$ ), ( $c$ ) and ( $d$ ) are zoomed-in views of regions V1, V2 and V3 in panel ( $a$ ), respectively, where blue arrows indicate velocity vectors.

The regular configurations, characterised by spontaneously selected near-integer values of the normalised gap spacing, $\tilde {G}_{i,i+1}$ , arise as a consequence of wake–plate–wall interactions. Figure 7 displays the vorticity contour for the regular mode that emerges when $N=36$ and $H=3.5$ . The wake vortex from the first plate exhibits a 2P mode, shedding two vortices of the same sign within one half-cycle period (Zhu et al. Reference Zhu, He and Zhang2014), and interacts with the second plate, influencing its vortex formation and merging with its wake vortex. The first row of figure 7( $a$ ) illustrates that the vortex environments surrounding the first few following plates are distinct, explaining the varying gap distances observed between some adjacent plates in figure 2( $c$ ). As the plate index $i$ increases, the vortex environment for each plate gradually converges to a spatially periodic state. Figures 7( $b$ ), 7( $c$ ) and 7( $d$ ) depict the vortex structures behind the 3rd, 12th and 30th plates, respectively. These figures illustrate the evolution of the wake vortex from a mixed state to a 2S mode. As $i$ increases, the plate’s trajectory in the vortex street eventually transitions to a vortex-slalom mode. (see also figure 9).

Figure 8.

Quantitative measurements of the regular-mode case ( $N=36$ , $H=3.5$ ). ( $a$ ) Schematic of the vortical structures generated by the plates. ( $b$ ) Temporal evolution of the horizontal force for selected plates: plates 1, 9, 18, 27, 36. ( $c$ ) Normalised longitudinal distance $\tilde {D}_{i,i+1}=D_{i,i+1}/\lambda$ as a function of plate index $i$ ( $i=1$ , $2$ , …, $N-1$ ). ( $d$ ) Vorticities at marked points $(x_{L,i},0)$ and $(x_{T,i},0)$ as a function of plate index $i$ . ( $e$ ) Integrals of vorticity ( $\omega _{\textit{Int}}$ ), absolute vorticity ( $|\omega |_{\textit{Int}}$ ) and enstrophy $E_{\textit{Int}}$ in domain $\varOmega _i$ as functions of plate index $i$ . ( $f$ ) Horizontal components of impulse force $I_i$ and vortex force $L_i$ in domain $\varOmega _i$ as functions of plate index $i$ . Panels ( $c$ )–( $f$ ) present results at the non-dimensional time $t/T = 197$ .

To quantitatively examine the swimming performance of the plates and the spatial-periodicity characteristics of the surrounding flow field, we analyse the system with increasing plate index $i$ in the large- $N$ regular mode ( $N=36$ , $H=3.5$ ). Figure 8 presents measurements of several key parameters. To clarify the calculation methodology, figure 8( $a$ ) presents a schematic diagram of flexible plates moving within vortex streets. Here, $x_{L,i}$ and $x_{T,i}$ denote the horizontal coordinates of the leading and trailing edges of the $i$ th plate, respectively; $D_{i,i+1}=x_{L,i+1}-x_{L,i}$ represents the longitudinal distance between the leading edges of the $i$ th and $(i+1)$ th plates; and $\varOmega _i$ indicates the integral domain bounded by $[x_{L,i},x_{L,i+1}]\times [-H/2, H/2]$ . Figure 8( $b$ ) displays the temporal evolution of horizontal hydrodynamic forces for representative plates (plates 1, 9, 18, 27 and 36). While plates 1, 9 and 18 exhibit significant force variations, Plates 18–36 demonstrate nearly identical force profiles, explaining the convergence of normalised gap distance $\tilde {G}_{i,i+1}$ to a single curve in figure 2( $c$ ). Figure 8( $c$ – $f$ ) captures a snapshot of the flow field at $t/T = 197$ , corresponding to an instant after the schooling system attains a quasi-steady state. The normalised longitudinal distance $\tilde {D}_{i,i+1}$ in figure 8( $c$ ) shows convergence at $i=18$ , consistent with the force results in figure 8( $b$ ). Vorticity measurements at marked points $(x_{L,i},0)$ and $(x_{T,i},0)$ in figure 8( $d$ ) reveal convergent behaviour with increasing $i$ , with faster convergence at $(x_{L,i},0)$ .

A detailed analysis of the vorticity function integrals is conducted within the vortex-pair domain $\varOmega _i$ , including: circulation $\omega _{\textit{Int},i}=\int _{\varOmega _i}\omega {\rm d}V$ , absolute vorticity $|\omega |_{\textit{Int},i}=\int _{\varOmega _i}|\omega | {\rm d}V$ , enstrophy $E_{\textit{Int},i}=\int _{\varOmega _i}\omega ^2/2 {\rm d}V$ , horizontal components of impulse force $I_i=\int _{\varOmega _i}{\boldsymbol{x}}\times {\boldsymbol \omega }\boldsymbol{\cdot }{\boldsymbol{e}}_x {\rm d}V$ and vortex force $L_i=\int _{\varOmega _i}{\boldsymbol{u}}\times {\boldsymbol \omega }\boldsymbol{\cdot }{\boldsymbol{e}}_x {\rm d}V$ , where ${\boldsymbol{e}}_x$ is the unit vector in the horizontal direction. The observed convergence of circulation ( $\omega _{\textit{Int},i}$ ), absolute vorticity ( $|\omega |_{\textit{Int},i}$ ) and enstrophy ( $E_{\textit{Int},i}$ ) with increasing $i$ indicates that a dynamic equilibrium is achieved between the vortex generation/evolution induced by the flapping motion of flexible plates and the viscous dissipation of the vortex field under wall effects. According to impulse theory (Kang et al. Reference Kang, Liu, Su and Wu2017), the convergence of impulse $I_i$ and vortex force $L_i$ directly explains the force profile convergence in figure 8( $b$ ). The results in figure 8 demonstrate the intrinsic coupling between flow field structure convergence and force characteristic convergence with increasing $i$ .

Figure 9.

The instantaneous vorticity contours for cases with varying channel width (e.g. $H=1.5$ , $2.0$ , $2.5$ , $3.0$ , $3.5$ , $4.0$ ) and fixed group size ( $N=36$ ). The plates are depicted by green curves. The flow field region extends $55$ plate length forward from the vicinity of the last plate, corresponding to $x^*\in [0,55]$ . For the regular mode at $H\geqslant 2.5$ , the black curves superimposed on each vorticity contour plot represent the trajectory of the leftmost flexible plate within the vortex street.

Figure 9 compares the vorticity fields of a schooling system ( $N=36$ ) under varying channel widths ( $H = 1.5$ , $2.0$ , $2.5$ , $3.0$ , $3.5$ and $4.0$ ). Based on the results in figure 7 – where the spatial periodicity of vortices and the convergence of plate kinematics primarily occur for larger plate indices ( $i$ ) – we display the flow field region extending $55$ body lengths forward from the vicinity of the last plate. For $H \leqslant 2.0$ (hybrid mode), the vorticity field appears relatively disordered, explaining the observed kinematic variability among flexible plates. For $H\geqslant 2.5$ , the vortex flow field exhibits a distinct spatial periodicity similar to that shown in figure 7. Moreover, as $H$ increases, the vortex structure transitions progressively from a $2S$ mode to a $2P$ mode as the channel width increases. For the regular mode at $H\geqslant 2.5$ , the black curve on each vorticity contour depicts the trajectory of the leftmost flexible plate in the vortex street. The trailing plates downstream follow this trajectory, demonstrating a vortex-slalom mode. The vortex-slalom mode here describes a flexible plate slaloming between vortex cores. Unlike the energetically unfavourable vortex-slaloming case reported by Zhu et al. (Reference Zhu, He and Zhang2014) (where tandem self-propelled plates avoided this mode), our observed path deviates by approximately half a wavelength in the vortex street. As demonstrated in § 3.4, this configuration actually reduces energy consumption.

Figure 10.

( $a$ ) Schematic of a bounded unsymmetric double vortex row in a channel of width $H$ , where $a$ denotes the inter-row spacing and $L_{vor}$ the intra-row vortex spacing. ( $b$ ) Vortex spacing $L_{vor}$ versus channel width $H$ . Red solid circles correspond to cases with $N=36$ and $H=2.5$ , $3.0$ , $3.5$ , $4.0$ . The green line represents the function $L_{vor}=1.419H$ .

The stability of a confined vortex street was theoretically investigated by Rosenhead (Reference Rosenhead1929), as illustrated in figure 10. The configuration consists of two infinite parallel vortex rows with anti-symmetric circulations, where each vortex in one row is aligned with the midpoint between two consecutive vortices in the other row. Here, $a$ is the distance between the rows and $L_{vor}$ is the spacing between consecutive vortices, with the channel width being $H$ . His analysis demonstrates that when $L_{vor}\geqslant 1.419 H$ , the system is stable for all values of $a$ . As is shown in figure 9, the vortex structures surrounding followers in large-N groups under regular mode exhibit spatial periodicity, closely resembling the bounded unsymmetrical double-row configuration described by Rosenhead (Reference Rosenhead1929) (although not forming a perfectly infinite vortex street). We measured the vortex spacing $L_{vor}$ in four regular-mode cases ( $N=36$ and $H=2.5$ , $3.0$ , $3.5$ , $4.0$ ) shown in figure 10( $b$ ). It is seen that the vortex spacing $L_{vor}$ gradually decreases with increasing $H$ in the four measured cases. The green line in figure 10( $b$ ) represents the function $L_{vor} = 1.419H$ , with the left side of the line indicating the presence of stable vortex streets and the right side indicating the inability to form them. Our measurements confirm that for $H\leqslant 4.0$ all cases fall within the stable regime. Extrapolation suggests that cases with $H\geqslant 4.5$ would enter the unstable regime, demonstrating that excessively large $H$ values prevent periodic vortex-street formation and consequently inhibit large-scale collective cohesion. This explains why in figure 3, for groups with larger- $N$ , the upper boundary of the regular mode and unstable mode occur at $H\in (4.0, 4.5)$ .

Figure 11.

( $a$ ) Schematic of the wake structures generated by the plates in the $N=36$ system. Vortices shed directly from the plates are denoted by $\omega _I$ , while the induced vortices near the upper and lower walls are labelled $\omega _{\textit{II}*}$ and $\omega _{\textit{II}**}$ , respectively. Dashed boxes in the wake of the last plate mark regions $V1$ , $V2$ , $V3$ and $V4$ . ( $b$ ) Instantaneous vorticity contours in the wake of the last plate after the flow has reached a quasi-steady state, for varying channel widths ( $H=1.5$ , $2.0$ , $2.5$ , $3.0$ , $3.5$ , $4.0$ ) at a fixed group size ( $N=36$ ). ( $c$ ) Circulation as a function of channel width. The curve labelled $\varGamma _{I,Vi}$ is computed via $\int _{Vi}|\omega _I|{\rm d}V$ , and $\varGamma _{\textit{II},\textit{Vi}}$ represent $\int _{Vi}|\omega _{I*}+\omega _{I**}|{\rm d}V$ .

In a stable collective system with group size $N$ , the vortex structure in the trailing plate’s wake directly governs the dynamical behaviour of the $(N+1)$ th plate, as manifested by its capacity to maintain stable spacing and form ordered configurations. As shown in figure 11, we analyse the wake vortex structures of stable $N=36$ configurations across varying channel widths. As demonstrated by Jeong et al. (Reference Jeong, Lee and Park2021), the vortex streets shed by the plates maintain near-perfect alignment along the centreline. For each vortex ( $\omega _I$ ) in the street, counter-rotating induced vortices ( $\omega _{\textit{II}*}$ , $\omega _{\textit{II}**}$ ) emerge near the upper and lower walls, as illustrated in figure 11( $a$ ). To quantify circulation variations with channel width and downstream position, we examined four equidistant regions immediately downstream of the trailing plate (figure 11 $b$ ). Each region, spanning half the horizontal spacing between vortex pairs ( $L_{vor}/2$ ), contains one primary vortex ( $\omega _I$ ) and two counter-rotating induced vortices ( $\omega _{\textit{II}*}$ , $\omega _{\textit{II}**}$ ). The circulation strength of the primary vortex is calculated as $\varGamma _{I,Vi}=\int _{Vi}|\omega _I| {\rm d}V$ , while that of the induced vortices is determined by $\varGamma _{\textit{II},\textit{Vi}}=\int _{Vi}|\omega _{\textit{II}*}+\omega _{\textit{II}**}|{\rm d}V$ . It is seen from figure 11 that vortex strength decays downstream for both primary vortex and induced vortices. Except at $H=1.5$ (where trailing plate sub-grouping occurs), $\varGamma _I$ increases with channel width. For the circulation of induced vortices, two distinct regimes are observed: (i) near-equivalent strength to the primary vortices when $H \leqslant 2.5$ , and (ii) monotonic decay with increasing channel width for $H \geqslant 2.5$ .

3.4. Schooling stability

To investigate how channel width and group size affect collective stability, we introduce a $(N+1)$ th plate with identical flapping kinematics (2.1) into the wake of a stable $N$ -plate system and examine its tracking stability. By varying the spacing between the $(N+1)$ th plate and the $N$ th plate, distinct vortex-street interactions are induced. Through investigation of the net horizontal forces ( $F_x$ ) acting on the $(N+1)$ th plate under different wake–plate interactions, we analyse its stability within the vortex street generated by the $N$ -plate system. In our numerical simulations, it is challenging to maintain precise spacing control between the additional $(N+1)$ th plate and the $N$ th plate under self-propelled free-motion conditions. Since the tracking stability of the $(N+1)$ th plate is predominantly governed by its hydrodynamic interactions with the $N$ th plate’s vortex street, we employ the quantified instantaneous wake vorticity field (figure 11) as a prescribed background flow. By varying the distance between this additional plate and the $N$ th plate, we can modify its trajectory within the vortex wake. Furthermore, we estimated the relative velocities and hydrodynamic forces during these trajectory changes with respect to the vortex street.

Figure 12.

( $a$ ) Schematic of possible trajectories for a trailing plate in the $N$ th plate wake. Slalom mode (navigating between vortices, lateral-flow aligned) and vortex-locking mode (penetrating vortices, lateral-flow opposed). Vortex polarity: red/blue arrows (positive/negative). ${\boldsymbol{V}}_{\!B}=(u_B, v_B)$ is trailing plate mid-point velocity; ${\boldsymbol{V}}_{\!\varGamma} =(u_\varGamma , v_\varGamma )$ vortex-induced velocity; $D_m$ spacing between the $N$ th plate and the trailing plate. ( $b$ ) Velocity vectors $(0, v_R)$ and $(0,u_\varGamma )$ corresponding to slalom-mode and vortex-locking-mode trajectories. The background vortex field is extracted from the $N=36$ , $H=3.5$ case. Black curves represent the plate trajectories within the vortex street. The upper subpanel displays the $(0, v_R)$ vector field along trajectories, where $v_R=v_B-v_\varGamma$ (black arrows: slalom mode; green arrows: vortex-locking mode), while the lower subpanel shows the $(0, u_\varGamma )$ distribution. ( $c$ ) Estimated mean thrust $\bar {F}_T$ and drag $\bar {F}_D$ on the hypothetical $(n+1)$ th plate versus its spacing from the $N$ th plate. $c_t$ and $c_d$ are corresponding thrust and drag coefficients, respectively. The normalised spacing $\hat {D}_m$ is defined as $D_m/L_{vor}$ . ( $d$ ) Time-dependent forces $F_x$ on the $N$ th plate (simulated) and $(N+1)$ th (estimated) plate at $\hat {D}_m=1$ ( $N=36,H=3.5$ ).

Figure 12( $a$ ) illustrates two possible trajectories of the $(N+1)$ th plate in the vortex street: the vortex-slalom mode and the vortex-locking mode. In the slalom mode, the $(N+1)$ th plate navigates between vortices, replicating the trajectory of the $N$ th plate. This synchronisation maintains an inter-plate spacing $D_m$ that equals integer multiples of the vortex spacing $L_{vor}$ . In contrast, the vortex-locking mode features vortex-core penetration. Here, we define ${\boldsymbol{V}}_{\!B}=(u_B, v_B)$ as the velocity of the $(N+1)$ th plate’s centre point, and ${\boldsymbol{V}}_{\!\varGamma} =(u_\varGamma , v_\varGamma )$ as the vortex-induced velocity at the same position. Here, ${\boldsymbol{V}}_{R} \equiv {\boldsymbol{V}}_{\!B} - {\boldsymbol{V}}_{\!\varGamma} = (u_R, v_R)$ represents the relative velocity between the $(N+1)$ th plate and the local surrounding fluid. Notably, the plate’s lateral motion aligns with the flow direction in slalom mode but opposes it in vortex-locking mode.

Figure 12( $b$ ) systematically compares the kinematic signatures of vortex-slalom and vortex-locking modes for the hypothetical $(N+1)$ th plate in a $N=36$ , $H=3.5$ configuration. The slalom mode manifests flow–plate synchronisation, evidenced by attenuated lateral relative velocity ( $|v_R|^{\textit{max}}/|v_B|^{\textit{max}} \approx 0.23$ ) and persistent streamwise induction ( $u_\varGamma /u_B \approx 0.5$ ). This mode may promote energy recovery through constructive vortex interaction. The vortex-locking mode generates strong lateral flow opposition ( $|v_R|^{\textit{max}}/|v_B|^{\textit{max}} \approx 1.85$ ) with disruptive $u_\varGamma$ fluctuations ( $|u_\varGamma |/\bar {u}_B\lt 0.17$ , where $\bar {u}_B$ is the mean value of $u_B$ over one period), which may cause larger energy consumption compared with an isolated fish swimming.

For hydrodynamic force quantification on the $(N+1)$ th plate, we established a reduced-order model of the wake–plate interaction dynamics. In this model, the $(N+1)$ th plate is kinematically reduced to a point mass, with its prescribed motion governed by the position and velocity of the plate’s centre point. We assumed that the $(N+1)$ th plate experiences a drag ( $F_D$ ) proportional to the $1.5$ power of its horizontal velocity relative to the ambient fluid at the plate’s centre point location, i.e. $F_D=c_d u_{R}^{1.5}$ (Oza et al. Reference Oza, Ristroph and Shelley2019). Similarly, the thrust ( $F_T$ ) is assumed to be proportional to the square of its vertical velocity component relative to the ambient fluid at that location, i.e. $F_T=c_t v_{R}^2$ (Newbolt et al. Reference Newbolt, Zhang and Ristroph2019). Here, $c_t$ and $c_d$ are constant thrust and drag coefficients, respectively. The average thrust and drag over a flapping period are $\bar {F}_T= {1}/{T}\int _{t'}^{t'+T}F_T {\rm d}t$ , $\bar {F}_D= {1}/{T}\int _{t'}^{t'+T}F_D {\rm d}t$ , respectively. Figure 12( $c$ ) illustrates the variation of $\bar {F}_T/c_t$ and $\bar {F}_D/c_d$ with respect to the normalised spacing $\hat {D}_m$ , where $\bar {F}_D/c_t= {1}/{T}\int _{t'}^{t'+T} u_R^{1.5} {\rm d}t$ , $\bar {F}_T/c_t= {1}/{T}\int _{t'}^{t'+T} v_R^{2} {\rm d}t$ , $\hat {D}_m=D_m/L_{vor}$ . As $\hat {D}_m$ increases, both $\bar {F}_T/c_t$ and $\bar {F}_D/c_d$ exhibit periodic variations with a unit period, showing synchronised trends. When $\hat {D}_m$ takes integer values (i.e. coinciding with the vortex-slaloming trajectory of the $N$ th plate), the thrust and drag approach their minima. We define these integer- $\hat {D}_m$ configurations as the equilibrium path (where thrust and drag balance each other), which serves as the basis for determining the thrust and drag coefficients.

To validate the reliability of our simplified force prediction model, we estimated the temporal variation of the horizontal force $F_x = -F_T + F_D$ acting on the $(N+1)$ th plate along the equilibrium path at $\hat {D}_m=1$ , shown in figure 12( $c$ ). We compared this result with the computed force time history of the $N$ th plate from numerical simulations. The results show close agreement between the two curves in their temporal variation patterns, confirming the reliability of our reduced-order model.

Figure 13.

Schooling stability comparison for the regular modes: ( $a$ , $b$ ) hydrodynamic force $\bar {F}_x$ and ( $c$ ) potential $\varPhi (\hat {D}_m)=-\int \bar {F}_x d(\hat {D}_m)$ as functions of $\hat {D}_m$ . For each case, the stable trajectory (blue curve) that satisfies $\bar {F}_x=0$ is superimposed on the inset vorticity plots of panels ( $a$ ) and ( $b$ ); $\Delta \varPhi _L$ and $\Delta \varPhi _R$ in panel ( $c$ ) represent the potential well depths on the left and right sides of the equilibrium point at $\hat {D}_m=1$ , respectively.

For regular modes, we calculated the horizontal hydrodynamic force ( $\bar {F}_x=-\bar {F}_T+\bar {F}_D$ ) of the $(N+1)$ th plate as a function of the normalised spacing ( $\hat {D}_m$ ) in several representative wake vortex fields (see the insert vorticity plots in figure 13 $a$ , $b$ ). Referencing the case with $N=36, H=3.5$ , we determined the thrust and drag coefficients based on the assumption that $\bar {F}_x=-\bar {F}_T+\bar {F}_D=0$ when $\hat {D}_m=1$ . Figure 13( $a,b$ ) shows $\bar {F}_x$ on the plate varying with $\hat {D}_m$ in these representative wake flow fields. At $\hat {D}_m=1$ (corresponding to the trajectory represented by the blue curve in the inset vorticity plots), all cases satisfy $\bar {F}_x\approx 0$ . Near the equilibrium position, the hydrodynamic force versus distance will satisfy a spring-like relationship, i.e. $\bar {F}_x\approx -k(\hat {D}_m-1)$ .

In addition, the emergence of equilibrium positions may be illustrated in terms of a hydrodynamic potential, which is defined as the integral of force with respect to distance, i.e. $\varPhi (\hat {D}_m)=-\int \bar {F}_x d (\hat {D}_m)$ (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016; Peng et al. Reference Peng, Huang and Lu2018c ). Figure 13( $c$ ) shows the hydrodynamic potential $\varPhi$ as a function of $\hat {D}_m$ . It can be observed that there are stable wells in the potential energy landscape. The bottoms of these wells correspond to the stable equilibrium states, with $\hat {D}_m$ . The potential well depth $\Delta \varPhi _L$ (between the $\hat {D}_m=1$ minimum and left peak) characterises the energy barrier against forward displacement beyond equilibrium due to wake-induced repulsion. As shown in figure 13( $c$ ), these values show negligible variations across all cases. Correspondingly, the downstream well depth $\Delta \varPhi _R$ (between the well bottom and adjacent right peak) quantifies the restoring potential energy when hydrodynamic attraction recovers a lagging plate’s position. It can be seen that for a fixed group size $N$ , an optimal channel width of $H=3.5$ is most effective for schooling stability, as excessively larger or smaller widths of channel hinder the formation of stable regular modes. Moreover, when $H$ is fixed at 3.5, increasing the cluster size enhances stability for the regular modes – a finding that highlights a fundamental contrast with unconfined flow (Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024; Nitsche et al. Reference Nitsche, Oza and Siegel2025). Whereas increasing the group size in unbounded flow amplifies downstream perturbations and triggers flonon instabilities through complex vortex interactions (Newbolt et al. Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024; Nitsche et al. Reference Nitsche, Oza and Siegel2025), confinement within an optimal channel width ( $H \approx 3.5$ ) promotes spatially periodic vortex streets. In this structured wake, the plates spontaneously adopt a vortex-slalom trajectory with inherent self-restoring stability, thereby explaining how larger groups can enhance collective cohesion.

Figure 14.

Propulsive performance comparison for the regular modes: ( $a$ ) estimated power function $\hat {P}(\hat {D}_m)= {1}/{T}\int _{t'}^{t'+T}|v_R|^3 {\rm d}t$ for the $(N+1)$ th plate versus $\hat {D}_m$ at varying channel widths ( $H=2.5$ , $3.0$ , $3.5$ , $4.0$ ) with fixed group size ( $N=36$ ). ( $b$ ) Estimated power function $\hat {P}(\hat {D}_m=1)$ versus input power $P_N$ for the $N$ th plate. ( $c$ ) Estimated efficiency function $\hat {\eta }= ({1}/{2})MU_c^2/(\hat {P}T)$ at $\hat {D}_m=1$ versus efficiency $\eta _N$ for the $N$ th plate.

To determine whether the vortex-slalom path formed spontaneously in the stable regular mode is energetically favourable, we estimated the energy cost versus $\hat {D}_m$ . Here, we assume the energy consumption is proportional to the cube of the relative flapping velocity (Newbolt et al. Reference Newbolt, Zhang and Ristroph2019; Liu, Liu & Huang Reference Liu, Liu and Huang2022), i.e. $P(\hat {D}_m)\sim \hat {P}(\hat {D}_m)= {1}/{T}\int _{t'}^{t'+T} |v_R|^3 {\rm d}t$ . Figure 14( $a$ ) compares the energy consumption function $\hat {P}$ of the $(N+1)$ th plate in the wake of an $N$ -plate system ( $N=36$ ) versus its spacing $\hat {D}_m$ from the $N$ th plate, for different channel widths ( $H=2.5$ , $3.0$ , $3.5$ , $4.0$ ). The results show that $\hat {P}$ varies periodically with $\hat {D}_m$ , reaching near-minimum values at equilibrium positions ( $\hat {D}_m=1,2$ ), confirming the energetic advantage of spontaneously formed vortex-slalom paths. Figure 14( $b$ ) presents the energy consumption function $\hat {P}(\hat {D}_m=1)$ at equilibrium versus channel width $H$ . For the regular mode ( $H\in [2.5, 4]$ ), $\hat {P}$ first decreases then increases, matching the energy consumption trend ( $P_N$ ) observed in numerical simulations for the $N$ th plate. Furthermore, we calculated the efficiency function $\hat {\eta }= ({1}/{2})MU_c^2/(\hat {P}T)$ at $\hat {D}_m=1$ versus $H$ , finding its variation consistent with the efficiency trend ( $\eta _N$ ) of the $N$ th plate in simulations. This explains why the high-efficiency phenomenon for large- $N$ regular-mode groups occurs at $H\in [3,3.5]$ in figure 6.

Figure 15.

( $a$ ) Temporal evolution of the power consumption $P$ for the last plate ( $i=36$ ) in the 36-plate system for the regular-mode state. Labels $t_1$ – $t_4$ denote the peaks of the $P$ – $t$ curve for $H = 2.5,\;3.0,\;3.5,\;4.0$ , respectively, and $t_5$ – $t_8$ denote the corresponding troughs. ( $b$ ) Instantaneous vorticity contours for the eight instants marked in panel ( $a$ ). The plate-averaged velocity $\tilde {{\boldsymbol{u}}} = \int _0^1 {\boldsymbol{u}}(s,t)\text{d}s$ and averaged force $\tilde {{\boldsymbol{F}}} = -\int _0^1 {\boldsymbol{F}}_s(s,t) \text{d}s$ are shown as vectors, where ${\boldsymbol{F}}_s$ is the Lagrangian fluid force on the plate. Also, $\theta$ is the angle between $\tilde {{\boldsymbol{u}}}$ and $\tilde {{\boldsymbol{F}}}$ . ( $c$ ) Statistics of $|\tilde {{\boldsymbol{u}}}|$ , $|\tilde {{\boldsymbol{F}}}|$ , $\theta$ and $\tilde {{\boldsymbol{u}}}\boldsymbol{\cdot }\tilde {{\boldsymbol{F}}}$ for the eight cases $t_1$ – $t_8$ .

Figure 15( $a$ ) shows that for $H=3.5$ both the peak and the trough of the instantaneous power of the 36th foil in the array are the lowest among all tested channel widths, confirming the minimal mean power consumption reported in figure 14( $b$ ). Figure 15( $b$ ) shows the vorticity fields together with the foil-averaged velocity $\tilde {{\boldsymbol{u}}}$ and force $\tilde {{\boldsymbol{F}}}$ vectors at the instants of peak and trough power. The instantaneous power can be estimated by $P\sim \tilde {{\boldsymbol{u}}}\boldsymbol{\cdot }\tilde {{\boldsymbol{F}}}$ , which is positive when the angle $\theta$ is acute. Figure 15( $c$ ) quantifies $|\tilde {{\boldsymbol{u}}}|$ , $|\tilde {{\boldsymbol{F}}}|$ , $\theta$ and $\tilde {{\boldsymbol{u}}}\boldsymbol{\cdot }\tilde {{\boldsymbol{F}}}$ for all cases; the product $\tilde {{\boldsymbol{u}}}\boldsymbol{\cdot }\tilde {{\boldsymbol{F}}}$ is minimised at both the peak and the trough for $H\approx 3.5$ . The propulsive speed typically increases with the narrowing of the channel (reducing $H$ ), while the power shows a distinct minimum close to $H\approx 3.5$ (figure 14 $b$ ), resulting in maximum efficiency at this confinement width (figure 14 $c$ ). Physically, the minimum in power consumption and maximum in the efficiency at $H\approx 3.5$ stem from a favourable kinematic–dynamic synergy between the foil motion and the induced flow, which minimises the effective lateral speed of the foil relative to the fluid and thereby reduces the power required for flapping.

Figure 16.

Hybrid mode comparison: ( $a$ , $b$ ) hydrodynamic force $\bar {F}_x$ and ( $c$ ) potential $\varPhi (\hat {D}_m)=-\int \bar {F}_x d(\hat {D}_m)$ as a function of $\hat {D}_m$ . The insets in panels ( $a$ ) and ( $b$ ) correspond to the wake vorticity contours for $N=12$ , $H=1.5$ and $N=12$ , $H=2.0$ , respectively.

The foregoing analysis suggests that while the $N$ th plate in a hybrid-mode $N$ -plate system generates a quasi-regular vortex street (figure 16 $a$ , $b$ insets), the $(N+1)$ th plate’s equilibrium position in this wake field may exhibit irregular spacing due to emergent subgroup dynamics (figure 9). Crucially, the hydrodynamic force and potential on the $(N+1)$ th plate remain measurable. As shown in figure 16( $a$ ), with $H$ fixed at $1.5$ , a stable equilibrium exists near $\hat {D}_m \approx 0.1$ . This stable position implies that the $(N+1)$ th plate tends to form a subgroup with the $N$ th plate. At larger values of $\hat {D}_m$ , the $(N+1)$ th plate experiences a positive hydrodynamic force for all tested cases of $N$ (i.e. $3$ , $12$ and $36$ ). This indicates a persistent drag effect that drives the plate rearward until vortex dissipation becomes sufficiently weak. This explains why the unstable region below the regular-mode zone in figure 3 corresponds to continuously increasing gap spacing. The enhanced dissipation at $H=1.5$ is evident in figure 11( $c$ ), where the circulation of the main vortex in $V1$ exceeds twice that in $V4$ . This implies that while the $(N+1)$ th plate cannot achieve equilibrium at $\hat {D}_m=1$ , it may stabilise at a position several equilibrium distances downstream. Figure 16( $c$ ) further confirms the absence of stable equilibrium positions for the $(N+1)$ th plate near $\hat {D}_m$ . For $H=2.0$ (figure 16 $b$ ), the $(N+1)$ th plate exhibits alternating positive/negative horizontal forces, with the $\bar {F}_x=0$ equilibrium position not necessarily coinciding with $\hat {D}_m=1$ . For instance, stable equilibrium positions exist near $\hat {D}_m \approx 0.1$ , where the plate tends to form a subgroup with the preceding plate. Notably, the force amplitudes in hybrid mode (figure 16 $a$ , $b$ ) are significantly smaller than in conventional mode ( $N=36$ , $H=3.5$ ), resulting in shallower potential wells (figure 16 $c$ ). This reduced potential well depth may promote subgroup formation and unstable mode.