4.1 Emergent dynamics and configurations

In our numerical settings, the in-phase flapping scenario ( $\unicode[STIX]{x1D719}=0$ ) was first investigated. The bending stiffness $K=1$ is chosen at first because it lies in the optimal range of stiffness for the largest propulsive speed for the isolated plate. The side-by-side initial configurations with $D_{0}=0$ and $H\in (0.2,6)$ are examined. Four typical states of collective pattern are identified. They are formed spontaneously by the two-plate system, i.e. the staggered-following (SF) mode, alternate-leading (AL) mode, non-periodic (NP) mode and non-interfering (NI) mode.

Figure 4.

Instantaneous vorticity contours at $t=(0/4)T$ , $T/4$ , $T/2$ and $(3/4)T$ from top to bottom. Panels (a–d) and (e–h) are the SF state ( $H=0.3$ , $K=1$ and $\unicode[STIX]{x1D719}=0$ ) and the AL state ( $H=0.9$ , $K=1$ and $\unicode[STIX]{x1D719}=0$ ), respectively. Blue and red colours denote negative (clockwise) and positive (anticlockwise) vorticity, respectively.

Figure 5.

The longitudinal spacing $D$ (a,c,e) and the relative propulsive speed $U_{r}$ of the two plates (b,d,f) as a function of time. Panels (a,b), (c,d), and (e,f) correspond to the cases of the SF, AL and NP states, respectively. The key parameters $H$ and $K$ are labelled in (a,c,e).

Figure 6.

Time-averaged longitudinal spacing $\overline{D}$ (a) and the peak-to-peak value $D_{pp}$ (b) as functions of $H$ for the in-phase flapping. The effect of $D_{0}$ is also given.

Figure 4(a–d) shows the vorticity contours at four typical phases in one flapping cycle for the case with $H=0.3$ and $K=1$ , corresponding to the SF mode. In the SF state, the following plate moves forward following the leading plate with a staggered configuration. The two plates behave like one single and longer plate and the vortex street is behind the follower. The time histories of the longitudinal spacing $D(t)$ and relative propulsive speed $U_{r}(t)=u_{1}(t)-u_{2}(t)$ for the SF mode are shown in figure 5(a,b), where $u_{1}(t)$ and $u_{2}(t)$ are the instantaneous horizontal velocity of plate 1 and 2, respectively. It is seen that $D(t)$ and $U_{r}(t)$ evolve to periodic states with $\overline{D}\approx 1.3$ and $\overline{U}_{r}=0$ after a few flapping periods, where $\overline{D}$ and $\overline{U}_{r}$ are average values over one flapping period.

Figure 4(e,f,g,h) shows the vorticity contours in the AL mode. The reversed Kármán vortex street is clearly shown in the wake. The shedding vortex pattern is different from that in the SF mode. In the AL mode ( $H=0.9$ ), $D(t)$ changes periodically with $\overline{D}=0$ (see figure 5 c,d), indicating that during the propulsion, the two plates leapfrog or move forward in an alternate-leading way. When the plates become harder, the non-periodic (NP) mode appears. Figure 5(e,f) shows that in the NP mode, $D(t)$ and $U_{r}(t)$ vary non-periodically with time. It is noted that the lateral flow interactions between the two plates are negligible when the lateral distance is large enough, e.g. $H=4.0$ , which is referred to as the non-interfering state.

The time-averaged longitudinal spacing $\overline{D}$ and peak-to-peak value $D_{pp}$ as functions of $H$ for $K=1$ are shown in figure 6. It is seen that the collective pattern of the two-plate system switches from the SF mode to the AL mode as $H$ at $H_{c}\approx 0.6$ . It is worth noting that the peak-to-peak value of $D$ , i.e. $D_{pp}$ , is significant in the SF and AL modes, implying dynamic adjustment in the horizontal direction between individuals due to the flow-mediated interactions. Additionally, $D_{pp}$ tends to be zero due to the lateral interference vanishing as $H\geqslant 4$ .

To further examine whether the motion state can be maintained when the system is perturbed, the initial configurations are adjusted by varying the initial longitudinal spacing $D_{0}$ . Figure 6 also shows the results for $D_{0}=0.5$ , 1.0 and 2.0. It is seen that the ultimate stable configurations are mainly determined by $H$ but independent of $D_{0}$ when the horizontal disturbance is small ( $D_{0}\leqslant 2.0$ ). Hence, the formed configurations due to the fluid mediation are stable.

Passive flexibility of flapping plates, which is closely associated with their propulsive performance and wake properties (Hua et al. Reference Hua, Zhu and Lu2013; Zhu et al. Reference Zhu, He and Zhang2014b ), may play a key role in collective locomotion through the flow interactions between individuals. Thus, we further examine the effect of the flexibility on the two-plate system. In our numerical setting, $K$ ranges from $0.1$ (very flexible) to $10$ (very hard). In the simulations, the side-by-side initial configurations are adopted, i.e. $D_{0}=0$ and $H\in (0.2,6)$ .

Figure 7.

Phase diagram for the four motion states ( $\unicode[STIX]{x1D719}=0$ ). Symbols ▵, ○, ▹ and ♢ represent the SF, AL, NP and NI states, respectively.

Our results indicate that the occurrence of the locomotion states of the plates depends mainly on the bending stiffness $K$ and the lateral spacing $H$ . The phase diagram for the four modes in the $H$ – $K$ plane is shown in figure 7. Each point in the figure represents a case we simulated. We do not intend to provide accurate borders between the modes, therefore it is only a schematic diagram. It is seen that the SF state occurs mainly in the region with small $H$ . Both the AL and NP states mainly appear at moderate $H$ but the NP state covers the region with larger $K$ (say $K>2.0$ ). When $H$ is large enough (say $H\geqslant 4.0$ ), only the NI state occurs because the flow-mediated interactions between the two plates are negligible.

4.2 Plate performances in the SF state and the AL state

According to our simulations, stable configurations can be spontaneously formed. The propulsive properties of the plates in the SF and AL modes are further investigated. To quantify the propulsive performance of the plates, three important indicators, namely, the mean propulsive speed, input work and propulsive efficiency (Hua et al. Reference Hua, Zhu and Lu2013; Zhu et al. Reference Zhu, He and Zhang2014a ) are defined.

The mean propulsive speed $U$ is defined as the time-averaged forward speed at the equilibrium state, i.e.

(4.1) $$\begin{eqnarray}\displaystyle U=\frac{1}{T}\int _{t^{\prime }}^{t^{\prime }+T}u(t)\,\text{d}t=-\frac{1}{T}\int _{t^{\prime }}^{t^{\prime }+T}\left(\left.\frac{\unicode[STIX]{x2202}X}{\unicode[STIX]{x2202}t}\right|_{s=0}\right)\text{d}t, & & \displaystyle\end{eqnarray}$$

where $T$ is the flapping period. The input work $W$ is required to maintain the self-propulsion of the flapping plates. It is computed as a time integral of the power $P$ performed by the surface of the body on the surrounding fluid over one flapping period, i.e.

(4.2) $$\begin{eqnarray}\displaystyle W=\int _{t^{\prime }}^{t^{\prime }+T}P(t)\,\text{d}t=\int _{t^{\prime }}^{t^{\prime }+T}\left[\int _{0}^{1}\boldsymbol{F}_{r}(s,t)\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{X}(s,t)}{\unicode[STIX]{x2202}t}\,\text{d}s\right]\text{d}t, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{F}_{r}$ represents the force on the surrounding fluid by the plates. To quantify the propulsive efficiency of a self-propulsive body, the ratio of the kinetic energy of the body and the input work has been employed. Therefore the propulsive efficiency is

(4.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}={\textstyle \frac{1}{2}}MU^{2}/W. & & \displaystyle\end{eqnarray}$$

The propulsive properties of the plates for the in-phase ( $\unicode[STIX]{x1D719}=0$ ) scenario are presented in figure 8. It is seen that there are three propulsion modes (the NI, AL and SF states, respectively) as $H$ decreases from 6.0 to 0.2 for $K<2.0$ . Figure 8(a) shows the mean propulsive speed $U$ versus $H$ at the typical bending stiffness $K=1$ , which is an optimal $K$ for an isolated flapping plate. The mean propulsive speed of an isolated plate with $K=1$ is also shown by the dot-dashed line for comparison. When $H>4.0$ (the NI state), the mean propulsive speed $U$ is very close to that of the isolated plate. In the AL mode, $U$ decreases with $H$ monotonically and reaches the lowest point at $H\approx 0.6$ . As $H$ decreases further, the plates switch to the SF state, and $U$ becomes significantly larger than that of the isolated plate.

Figure 8.

Propulsive speed $U$ (a), input work $W$ (b) and propulsive efficiency $\unicode[STIX]{x1D702}$ (c) as functions of $H$ for $K=1$ and $\unicode[STIX]{x1D719}=0$ . The green, grey and white regions represent the SF, AL and NI states, respectively.

Further, the input work and the propulsive efficiency as a function of $H$ are shown in figure 8(b,c). In the AL mode, the input work $W$ and propulsive efficiency $\unicode[STIX]{x1D702}$ of the plates are lower than those of the isolated one, while they jump to higher level in the region of the SF mode. It suggested that the flow-mediated interactions between the plates may be very different for the AL and SF states. Moreover, in the SF state, the plates have identical propulsive speed. Compared with the leading plate, the following one uses less input work (see figure 8 b) and therefore achieves higher propulsive efficiency (see figure 8 c). Hence, the SF configuration achieves better propulsive performance, and the follower can take advantage of the vortex shedding from the leading one in this compact configuration.

Figure 9.

Propulsive speed $U$ (a), average input work $(W_{1}+W_{2})/2$ (b) and average propulsive efficiency $(1/2)M(U_{1}^{2}+U_{2}^{2})/(W_{1}+W_{2})$ (c) as functions of $K$ for different $H$ ( $\unicode[STIX]{x1D719}=0$ ).

The effect of $K$ on the performance is shown in figure 9. The average values of input work and propulsive efficiency of the two plates are calculated to evaluate the performance of the two-plate system, i.e. $(W_{1}+W_{2})/2$ and $(1/2)M(U_{1}^{2}+U_{2}^{2})/(W_{1}+W_{2})$ , respectively. It is seen that, compared to the isolated plate ( $H=\infty$ ), the propulsive performance of the plates is significantly improved as $K$ increases when $K>0.6$ for $H=0.2$ and 0.5. However, for $H=0.7$ and 1.5, the propulsive speed and efficiency of the plates decrease with $K$ when $K>0.3$ . It is also noted that when $K>2.0$ , the mode switches to the non-periodic mode and it is not possible to calculate the average values of $U$ , $W$ and $\unicode[STIX]{x1D702}$ according to their definitions. Hence, there are no data for $K>2.0$ .

4.3 Plate performance in the SF mode and the tandem case

For the two in-phase flappings, in this section we will discuss whether the performance in the staggered following is better than the simple tandem case ( $H=0$ ). In the simulations of the simple tandem case, $H=0$ is considered. For the staggered cases, five cases with $H=0.1{-}0.5$ are simulated. In the cases of $H=0$ , 0.1 and 0.2, the initial longitudinal spacing $D_{0}=1.1$ is adopted to avoid collisions. As mentioned above, in all of the simulations, the ultimate longitudinal spacing $D$ is independent of $D_{0}$ (see figure 6). For the five cases with $H\in [0.1,0.5]$ , the plates adopt the staggered-following mode with the specified lateral distance $H$ . Figure 10(a,b) shows the propulsive speed $U$ and efficiency of the plate $\unicode[STIX]{x1D702}$ for the tandem case $(H=0)$ and staggered cases $(H\in [0.1,0.5])$ . It is seen from figure 10(a) that the propulsive speed in the case of $H=0.1$ , as well as that in the case of $H=0.2$ , is somewhat higher than that of the simple tandem case.

Figure 10.

Propulsive speed (a) and efficiency (b) of the plates for the simple tandem configuration $(H=0)$ and staggered configurations $(H=0.1{-}0.5)$ .

Figure 10(b) shows that at $H\approx 0.2$ the two plates achieve the optimal overall efficiency, which is better than that of the tandem case. It is also seen that in the tandem case ( $H=0$ ), the efficiency of the leading plate is higher than that of the following one. In this situation it seems that the leading plate can take advantage of the front–rear flow-mediated interaction. However, for the staggered cases with $H=0.2,0.3,0.4,0.5$ , the following plate can take more advantage than the leading one from the flow-mediated interaction. It is noted that there is an intermediate situation, i.e. at $H=0.1$ , the efficiencies of the two plates are almost identical.

For cases with multiple swimmers, e.g. four fish in Ashraf et al. (Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017), it seems that the swimmers adopt the SF mode rather than the tandem configuration. The possible reason is that in view of hydrodynamics, the overall efficiency is optimized at the compact staggered configuration, e.g. the case of $H=0.2$ . The follower energy saving due to the flow-mediated interaction in the SF mode also matches well with the actual observation about follower’s egoism (Landa Reference Landa1998). In summary, for the in-phase flapping, the compact staggered configuration with a small lateral spacing is preferred for the propulsive performance in comparison with the tandem case.

4.4 Forces and powers

In the following, the forces experienced by the plates and the produced power are analysed to reveal the propulsion mechanism. It is noticed that, although the plates experience the net drag or thrust in the $x$ -direction due to the fluid–structure interactions at a specific time, the time-averaged force over the cycle is zero. Previous studies (Thiria & Godoy-Diana Reference Thiria and Godoy-Diana2010; Ramananarivo, Godoy-Diana & Thiria Reference Ramananarivo, Godoy-Diana and Thiria2011) on the flapping flexible wing and plate have indicated that the flexible deformation plays an important role in the thrust and drag as a ‘shape factor’ that redistributes the contribution of the aero- or hydrodynamic forces in the direction of forward motion. Here we will explore quantitatively the flexible deformation influence on the distribution of thrust and drag and the role of the flow-mediated interactions in the collective motions.

The jump in the fluid force across the plate at a certain Lagrangian point, i.e. $\boldsymbol{F}_{s}$ can be decomposed into two parts: one is the normal force $\boldsymbol{F}^{n}$ in which the pressure component dominates, the other is the tangential force $\boldsymbol{F}^{\unicode[STIX]{x1D70F}}$ which comes from the viscous effects. These forces are defined as

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{s}=[-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D64F}]\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{F}^{n}+\boldsymbol{F}^{\unicode[STIX]{x1D70F}}, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}^{n}=(\boldsymbol{F}_{s}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}=(F_{x}^{n},F_{y}^{n}), & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}^{\unicode[STIX]{x1D70F}}=(\boldsymbol{F}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D749})\unicode[STIX]{x1D749}=(F_{x}^{\unicode[STIX]{x1D70F}},F_{y}^{\unicode[STIX]{x1D70F}}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the unit tensor, $\unicode[STIX]{x1D64F}$ the viscous stress tensor, $\unicode[STIX]{x1D749}$ the unit tangential vector toward the trailing edge, $\boldsymbol{n}$ the unit normal vector as shown in figure 11(a) in the schematic and $[~]$ denotes the jump in a quantity across the immersed boundary. The power $P$ can also be decomposed into the two parts: $P^{n}$ and $P^{\unicode[STIX]{x1D70F}}$ , which are contributed by $\boldsymbol{F}^{n}$ and $\boldsymbol{F}^{\unicode[STIX]{x1D70F}}$ , respectively, i.e.

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle P=\boldsymbol{F}_{s}\boldsymbol{\cdot }\boldsymbol{u}=P^{n}+P^{\unicode[STIX]{x1D70F}}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle P^{n}=\boldsymbol{F}^{n}\boldsymbol{\cdot }\boldsymbol{u}=P_{x}^{n}+P_{y}^{n}, & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle P^{\unicode[STIX]{x1D70F}}=\boldsymbol{F}^{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{u}=P_{x}^{\unicode[STIX]{x1D70F}}+P_{y}^{\unicode[STIX]{x1D70F}}, & \displaystyle\end{eqnarray}$$

where $P_{x}^{n}=F_{x}^{n}\cdot u$ , $P_{y}^{n}=F_{y}^{n}\cdot v$ , $P_{x}^{\unicode[STIX]{x1D70F}}=F_{x}^{\unicode[STIX]{x1D70F}}\cdot u$ and $P_{y}^{\unicode[STIX]{x1D70F}}=F_{y}^{\unicode[STIX]{x1D70F}}\cdot v$ . It is noted that to maintain the heaving motion at the leading edge of the plate, input power done by an external force is required to sustain the propulsion system. Suppose the input power is $P_{in}$ , we have $P_{in}+P_{x}^{n}+P_{y}^{n}+P_{x}^{\unicode[STIX]{x1D70F}}+P_{y}^{\unicode[STIX]{x1D70F}}=0$ .

Figure 11.

Schematic diagram for (a) force decomposition, (b) thrust generation, (c) phase difference between $\unicode[STIX]{x1D709}=\boldsymbol{F}_{s}\boldsymbol{\cdot }\boldsymbol{n}$ and $\tan \unicode[STIX]{x1D703}$ . Black curved lines represent the plate. $\unicode[STIX]{x1D749}$ and $\boldsymbol{n}$ denote the local tangential and normal vectors, respectively.

Figure 11(a) shows a schematic diagram of the force decomposition as mentioned above (see (4.4)–(4.6)). It is shown that the bending deformation of the flexible plate greatly influences the generation and distribution of $F_{x}^{n}$ due to different local geometry along the plate (Thiria & Godoy-Diana Reference Thiria and Godoy-Diana2010; Ramananarivo et al. Reference Ramananarivo, Godoy-Diana and Thiria2011). Actually, the normal force contribution to the thrust or drag depends on the phase relation between the normal force and the local slope of the flexible plate (Mysa & Venkatraman Reference Mysa and Venkatraman2016). A schematic diagram to illustrate this point is shown in figure 11(b,c). In figure 11(b), it is seen that when the normal force and local slope at a certain point of the plate are simultaneously positive ( $\unicode[STIX]{x1D709}=\boldsymbol{F}_{s}\boldsymbol{\cdot }\boldsymbol{n}>0$ and $\unicode[STIX]{x1D703}>0$ ) or negative ( $\unicode[STIX]{x1D709}<0$ and $\unicode[STIX]{x1D703}<0$ ), the normal force contributes to a thrust; otherwise, the normal force contributes to a drag. Moreover, figure 11(c) shows the curves of $\unicode[STIX]{x1D709}$ and $\tan (\unicode[STIX]{x1D703})$ for a segment on the plate in two periods. From the periodic curves, the phase difference between them ( $\unicode[STIX]{x1D6F9}$ ) can be estimated to characterize the phase relation between the normal force and the local slope. When $\unicode[STIX]{x1D6F9}=0$ , they vary synchronously and the normal force contributes to a thrust force throughout the cycle. When $\unicode[STIX]{x1D6F9}=\unicode[STIX]{x03C0}$ , the normal force contributes to a drag instead of the thrust throughout the cycle. When $\unicode[STIX]{x1D6F9}\in (0,\unicode[STIX]{x03C0})$ , as shown in figure 11(c), the normal force contributes to a drag in the grey region and a thrust in the white region throughout the cycle.

Figure 12.

Phase difference, normal force and local slope along the flexible plates for the in-phase $(\unicode[STIX]{x1D719}=0)$ scenario. (a) Phase difference between normal force and slope, (b) time-averaged absolute normal force, (c) time-averaged absolute slope and (d) time-averaged $F_{x}^{n}$ . The forces are normalized by $F_{ref}=(1/2)\unicode[STIX]{x1D70C}U_{ref}^{2}L$ .

For the in-phase $(\unicode[STIX]{x1D719}=0)$ scenario, the deformation and thrust in the isolated case, the SF mode and the AL mode are compared in the follows. The phase difference, normal force and local slope along the flexible plates are presented in figure 12(a–c), respectively. Similar to the situation in the isolated case, $F_{x}^{n}$ contributes to the thrust for the SF and AL states in figure 12(d).

For the SF mode, figure 12(a) shows that the phase difference is closer to zero than that of the isolated case. Compared to the isolated case, it is seen that $\overline{|\boldsymbol{F}^{n}|}$ along the plate is slightly enhanced (see figure 12 b) and the local bending deformation increases significantly (see figure 12 c). All these features result in a larger thrust contributed by the normal force in the SF state. Compared to the leading plate, although the following plate experiences smaller $\overline{|\boldsymbol{F}^{n}|}$ in the anterior part, a larger bending deformation in the posterior part leads to a larger $\overline{F_{x}^{n}}$ in the posterior part. The overall $\overline{F_{x}^{n}}$ of the following plate is larger than that of the leading one.

For the AL state ( $H=0.9,K=1.0$ ), it is seen from figure 12 that, compared to the isolated case, the phase difference of the plate in the AL state ( $\unicode[STIX]{x1D6F9}$ ) is far from zero and the bending deformation ( $\overline{|\text{tan}(\unicode[STIX]{x1D703})|}$ ) as well as the normal force ( $\overline{|\boldsymbol{F}^{n}|}$ ) are slightly smaller. They all lead to a lower contribution of the normal force to the thrust in the AL state.

Integrating the total work along the whole plate, we obtain the total work of the following, leading and the isolated plates, which are 4.35, 4.75 and 4.12, respectively. While the work done by $F_{x}^{n}$ (the useful work) is approximately 1.80, 1.66, 1.01, respectively. It is seen that more useful work is achieved for each plate in the SF mode compared to the isolated case. Hence, the propulsive efficiency of the plates in the SF state is higher than that of the isolated plate in figure 8(c).

4.5 Flow structure and the unsteady dynamics

To better understand the mechanisms underlying the propulsive performance, the vortical structures around the plates and pressure distributions on the surfaces are investigated in the following.

Figure 13.

Instantaneous pressure contours at (a) $t=(0/4)T$ , (b) $T/4$ , (c) $(2/4)T$ and (d) $(3/4)T$ for the SF mode ( $H=0.3$ , $K=1$ and $\unicode[STIX]{x1D719}=0$ ). Solid and dashed lines denote the positive and negative normalized pressure contours, respectively.

First, the instantaneous vorticity and pressure in the SF mode are analysed to understand how the flow-mediated interactions improve the propulsive performance. Figure 13 shows the pressure contours at typical instants of one flapping cycle. The two plates are propelled forward with a compact staggered configuration. From figure 13(a), it is seen that at the beginning of the down stroke, the thrust contributed by the normal force is negligible for the following plate because the pressure difference across the plate is small. The horizontal force in the SF mode is shown in figure 14. At $t=0$ , the forces $F_{x}$ of the two plates are almost identical.

Figure 14.

Horizontal force experienced by the plates (a,c) and power of the plates (b,d) as functions of time in one period. Panels (a,b) and (c,d) represent the SF mode and the AL mode, respectively. The integration of the force in one period is zero. The forces and powers are normalized by $F_{ref}=(1/2)\unicode[STIX]{x1D70C}U_{ref}^{2}L$ and $P_{ref}=\unicode[STIX]{x1D70C}U_{ref}^{3}L$ , respectively.

At $t=T/4$ , figure 13(b) shows that the pressure difference across the following plate is approximately equal to that of the leading one. Therefore their thrusts contributed by the normal force are almost identical. On the other hand, the positive vortex shed from the leading plate extends to the lower surface of the following plate. Due to the vortex inducing effect, there are negative pressure contours ahead the following plate (Vandenberghe, Zhang & Childress Reference Vandenberghe, Zhang and Childress2004), which also contribute to a thrust ( $F_{x}<0$ ). Hence, the overall thrust of the following plate is larger than that of the leading plate (see figure 14 a).

Figure 13(c) shows that the pressure difference across each plate and the deformation of each plate are small at $t=T/2$ . The pressure difference between the head and tail of each plate induces the drag force. From figure 14(a), we can see that both plates experience drag forces ( $F_{x}>0$ ) at this moment. At $t=(3/4)T$ , it is seen from figure 13(d) that there are more positive pressure contours ahead of the following plate than that ahead of the leading plate, which result in a larger drag force. Although the pressure difference across each plate produces thrust, the overall horizontal force for the leading (following) plates is the thrust (drag) force, as shown in figure 14(a). The power of each plate is shown in figure 14(b). It is seen that the total work done by the leading plate is slightly larger than that of the following one during one period. On the other hand, the plates move with an identical mean propulsive velocity. Hence, the following flapping plate achieves better efficiency and takes advantage of the wake of the leading one.

The instantaneous vorticity and pressure in the AL mode are analysed here. From figure 4, it is seen that in the AL mode, the leading-edge vortex is generated and convected downstream along the plates per half-cycle, which is similar to that of an isolated swimming plate. At $t=(0/4)T$ , the leading edges of the plates reach their maximum vertical positions, the upper and lower surfaces of the plates induce positive and negative vorticity in figure 4(e), respectively. There is a low pressure area near the vortices region in figure 15(a). Meanwhile, the leading edges of the plates are subject to a considerable drag due to the presence of the local high pressure region ahead of the plates. At this moment, both plates experience a net drag (see figure 14 c). As the plates start to stroke down ( $t=T/4$ ), i.e. figure 15(b), because the normal force due to the pressure difference acting on the lower plate is larger, the thrust of the lower plate is larger than that of the upper plate (see figure 14 c). Due to the larger thrust, the laggard lower plate begins to chase the upper leading one.

Figure 15.

Instantaneous pressure contours at (a) $t=(0/4)T$ , (b) $T/4$ , (c) $(2/4)T$ and (d) $(3/4)T$ for $H=0.9$ , $K=1$ and $\unicode[STIX]{x1D719}=0$ (the AL state). Solid and dashed lines denote the positive and negative normalized pressure contours, respectively.

When the plates stroke up, they experience the reverse procedure of the down stroke. In this way, the two plates would chase each other and achieve identical mean propulsive speed and efficiency over one flapping cycle. In the AL mode, compared to the isolated case, the normal force acting on each plate, and therefore the thrust, are reduced. Hence, the propulsive performance of each plate is poorer than that of the isolated case. The power of each plate is shown in figure 14(d). It is seen that the overall power is significantly smaller than that in the SF mode. The above mechanism is also simplified in § 4.6.2.

4.6 Antiphase $(\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0})$ scenario

The phase difference between the swimmers may play an important role in forming the stable flow-mediated configurations and achieving collective advantages. For example, Ashraf et al. (Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017) studied experimentally synchronization of tail-beating kinematics between neighbouring individuals within a group and found that fish swimming together at their free-ranging speed shows no consistent phase relationship between neighbours, which is consistent with the natural observations (Partridge & Pitcher Reference Partridge and Pitcher1979). On the contrary, fish in schools sustaining high swimming gaits shows strong synchronized dynamics, characterized by in-phase and out-of-phase (OP, or antiphase) modes (Ashraf et al. Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017) and the OP mode is favoured with respect to the in-phase state (Ashraf et al. Reference Ashraf, Godoy-Diana, Halloy, Collignon and Thiria2016). Moreover, the antiphase undulating motions are usually considered as one of the typical cases in schools and investigated in the literature (e.g. Dong & Lu Reference Dong and Lu2007; Boschitsch et al. Reference Boschitsch, Dewey and Smits2014; Dewey et al. Reference Dewey, Quinn, Boschitsch and Smits2014). In this section the effects of the phase difference on the emergent dynamics and configurations are investigated.

Figure 16.

Instantaneous vorticity (a–d) and pressure (e–h) contours at (a,e) $t=(0/4)T$ , (b,f) $T/4$ , (c,g) $(2/4)T$ and (d,h) $(3/4)T$ from the left column to right for $H=1.5$ , $K=3.0$ in the antiphase scenario. Blue and red colours denote negative (clockwise) and positive (anticlockwise) vorticity, respectively. Solid and dashed lines for pressure contours denote the positive and negative normalized pressure contours.

In the current setting, the key parameters in our simulations are $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0},D_{0}=0,H=1.2{-}2$ and $K=1{-}10$ . In the antiphase scenario, the minimum value of $H$ is slightly larger than $2A$ to avoid collision of two plates. Two motion modes are identified based on our numerical results, i.e. the moving abreast (MA) mode and the AL mode. The instantaneous vorticity and pressure contours at $t=(0/4)T$ , $T/4$ , $(2/4)T$ and $(3/4)T$ are shown in figure 16. It is seen that in the MA mode, $D(t)=0$ and the propulsive velocities of the two plates are identical, i.e. the two plates move side-by-side with an identical speed. Besides, the instantaneous snapshots of the flow field and the geometries of the two plates at any time are symmetric about $y=0$ . For the AL motion state, the two plates move forward in an alternate-leading way.

4.6.1 Phase diagram and propulsive properties

The phase diagram for the two motion modes are shown in figure 17. It is seen that in the $H$ – $K$ plane, the AL mode appears at the moderate bending stiffness $K$ with a small gap distance ( $H<1.5$ ). When the lateral distance $H$ is larger than 1.5, the AL mode disappears.

Figure 17.

Phase diagram for the MA (○) and AL (▵) modes in the $H$ – $K$ plane for $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ .

Figure 18.

Mean propulsive speed $U$ (a), input work $W$ (b), propulsive efficiency $\unicode[STIX]{x1D702}$ (c) as functions of $K$ for varying $H$ (antiphase flappings).

In the antiphase ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ) scenario, the mean propulsive speed $U$ versus the bending stiffness $K$ for $H\in (1.2,2.0)$ is shown in figure 18(a). For each $H$ , the speed $U$ increases rapidly with $K$ first and reaches a plateau region with a maximum $U$ , then it decreases rapidly at a critical bending stiffness and finally approaches a constant. The plateau region for the maximum $U$ becomes narrower as $H$ increases from 1.2 to 2.0. Note that for $H=\infty$ , i.e. the isolated plate case, the plateau region narrows into a peak. Hence, short-range flow-mediated interactions (smaller $H$ ) between the plates not only improve the propulsive speed, but also enlarge the range of optimal  $K$ . Moreover, figure 18(b,c) shows the input work and the propulsive efficiency as functions of $K$ for $H\in (1.2,2.0)$ . The input work increases with $K$ first and decreases rapidly to a constant at a critical $K$ . For the propulsive efficiency, it increases to a peak at $K\approx 0.4$ and then decreases to a constant at the critical $K$ . The decreasing section can be understood in the following way. In the plateau region of propulsion speed $U$ (approximately $K\in (1,4)$ ), the input work $W$ increases continuously to peak value. According to expression (4.3), it results in the decreasing propulsive efficiency. As a whole, figure 18 shows that the propulsive performance can be much improved when two closely spaced plates possess moderate flexibility.

4.6.2 Flow fields

Because the MA mode is able to achieve better propulsive performance, the vorticity and pressure contours of a typical MA case ( $H=1.5$ and $K=3.0$ ), that shown in figure 16, will be analysed first. In the following description, plate 1 (the upper plate) is taken as an example because the situation of plate 2 is symmetric to that of plate 1 about the line $y=0$ . From figure 16(a), it is seen that at the initial stage of down stroke of plate 1, the bending deformation increases to a maximum. The leading-edge positive vortex on the lower surface of plate 1 moves from the leading edge to the trailing edge, and interacts with the trailing-edge vortex. At $t=T/2$ in figure 16(b), all vortices are shed into the wake and the mutual induction among the vortices prevents them leaving vertically, which is helpful for the generation of thrust and improvement of propulsive efficiency (the detailed mechanism is shown in figure 19 and described in the next paragraph). Moreover, the high pressure region between the two plates grows when they stroke toward each other in figure 16(f), resulting in the enhancing normal force across the plates and the larger bending deformation. During the second half of the flapping cycle in figure 16(c,d,g,h), the vortices shed from the trailing edges move along the $x$ direction. The lower pressure region between two plates grows up gradually, which also enhances the normal force acting on the plate and the deformation of each plate, therefore the thrust increases. In the MA mode, the flow-mediated interactions contribute to a better propulsive performance as shown in figure 18.

Figure 19.

Instantaneous vorticity contours for $K=3.0$ (a) and $K=3.5$ (b) in the antiphase scenario. The other parameters are $H=1.5$ and $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ . Blue and red colours denote negative (clockwise) and positive (anticlockwise) vorticity, respectively. Red (blue) vectors represent the velocities induced by the neighbouring red (blue) vortices. Time-averaged jet in the wake for $K=3.0$ (c) and $K=3.5$ (d). The contours show the magnitude of adjusted horizontal velocity $\overline{u^{\prime }}=\overline{u}+U$ .

The vorticity contours at an instant for $K=3.0$ and $K=3.5$ are shown in figure 19. The induced velocity field around the vortices is also plotted in figure 19. It is seen that because the vortex ‘D1’ in figure 19(a) is weaker than that in figure 19(b) and the distance between vortices ‘U1’ and ‘D1’ in figure 19(a) is larger than that in figure 19(b), according to the Biot–Savart law, the induced lateral velocity acting on the ‘U1’ (the inclined vector) in the case of $K=3.0$ is smaller than that of ‘U1’ in the case of $K=3.5$ . Hence, the lateral movement of the vortices is slower in figure 19(a). In this way, it is seen clearly that the mutual induction among the vortices prevents them leaving vertically for the $K=3.0$ case. The comparison of the time-averaged jet in the wake for $K=3.0$ and $3.5$ is shown in figure 19(c,d). It is seen that the jet inclination angle for the $K=3.5$ case is approximately $15^{\circ }$ , which is larger than that in the $K=3.0$ case. This larger inclination angle is ‘a prognosticator of a further drop in thrust’ (Dong, Mittal & Najjar Reference Dong, Mittal and Najjar2006). Moreover, due to the small deformation of the plates in the case of $K=3.5$ , the thrust contributed by the normal force reduces significantly. Hence, the propulsive performance of $K=3.5$ becomes poorer.

Figure 20.

Schematic diagrams for the pressure contours in (a) the MA ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ) and (b) AL ( $\unicode[STIX]{x1D719}=0$ ) modes, respectively. Solid and dashed contours denote the higher pressure region (HPR) and lower pressure region (LPR), respectively.

In summary, the propulsive performance of the side-by-side configuration due to the in-phase and antiphase flappings can be understood in the following way. Figure 20(a) shows that, as the plates flap in antiphase, when they are flapping close to each other, the pressure between them is enhanced significantly because the fluid between them is squeezed. Hence, the pressure difference between the upper and lower surfaces of each plate is larger than that in the isolated case. The larger pressure difference results in a larger deformation and normal force, and therefore a larger thrust. Similarly, when they are flapping apart from each other, the low pressure region between them also enhances the normal force acting on each plate. Hence the propulsive performance of each plate is enhanced compared to the isolated case. For the in-phase flapping (see figure 20 b), similar to the isolated plate, higher pressure or lower pressure is generated below or above each plate as it flaps downward. Hence, in the region between them, the higher pressure (HPR) generated by the flapping of plate 1 will interact with the lower pressure (LPR) generated by the flapping of plate 2. The HPR and LPR would weaken each other. Hence, the pressure difference between the upper and lower surfaces of each plate is not so large as that in the isolated case. The smaller pressure difference results in a smaller normal force, and therefore a smaller thrust. Then the propulsive performance is weakened.

4.6.3 Rescaling propulsive speed and efficiency

As discussed in the above, at moderate bending stiffness and small $H$ (e.g. $H=1.2$ ), the thrust is mainly contributed by the normal force. Under this circumstance, the propulsive speed is enhanced and a better propulsive efficiency is achieved as shown in figure 18. It is also seen from figure 18(a,c) that for different $H$ , there exists a varying critical bending stiffness. When $K$ is beyond the critical value, the propulsive speed and efficiency of the plates decrease dramatically. For example, for $H=1.5$ the propulsive speed and efficiency decrease much as $K$ increases from 3.0 to 3.5. The possible reason is that the thrust contributed by the normal force is dominant in the case $H=1.5$ and $K=3.0$ . As $K$ increases to 3.5, the contribution of the normal force to the thrust reduces significantly due to the small deformation of the plates.

Figure 21.

Rescaled propulsive speed (a) and efficiency (b) of the flapping plates as functions of rescaled bending stiffness for varying $H$ . Legend is identical to that in figure 18.

Generally speaking, the bending deformation of the flapping plate is dominated by the fluid dynamic pressure on the surface, inertial force of the plate and elastic restoring force due to the flexibility. Because the mass ratio $M$ is fixed at 0.2 in our simulation, the inertial force of the plate is very minor due to the small $M$ . Hence, the fluid force and the elastic restoring force should be estimated to determine the deformation.

To examine the relative importance of the two forces, the ratio of the elastic force to the fluid force is often used. Usually the ratio is quantified by the effective stiffness $K=EI/\unicode[STIX]{x1D70C}U_{ref}^{2}L^{3}$ . However, here the fluid dynamic pressure, i.e. the normal force, is more prominent. Hence, the normal force instead of $\unicode[STIX]{x1D70C}U_{ref}^{2}L$ is proposed to normalize the elastic force $EI/L^{2}$ directly, i.e. $K^{\ast }=EI/\overline{|\boldsymbol{F}^{n}|}L^{2}$ , which is referred to as an effective stiffness, where $\overline{|\boldsymbol{F}^{n}|}$ is the time-averaged absolute normal force of the plates over one flapping period. Moreover, the characteristic velocity is proposed to be $U_{ref}^{\ast }=(2\overline{|\boldsymbol{F}^{n}|}/\unicode[STIX]{x1D70C}L)^{1/2}$ , and thus the normalized propulsive speed is $U^{\ast }=U/U_{ref}^{\ast }$ .

After the rescaling is performed, the propulsive speed $U^{\ast }$ and the efficiency $\unicode[STIX]{x1D702}$ of the antiphase flapping plates as functions of $K^{\ast }$ for different $H$ are shown in figure 21(a,b). From figure 21(a), it is seen that the plateaus with a maximum propulsive speed in figure 18(a) narrow into an identical peak at a critical value of $K^{\ast }\approx 0.25$ . Figure 21(b) shows that the curves of the propulsive efficiency also collapse together when $K^{\ast }$ is not too large. On the other hand, it is also noticed that the curves for different $H$ tend to deviate from each other slightly when $K^{\ast }$ is large. The possible reason is presented in the following. As discussed in the above, the thrust is primarily provided by the normal force for the plates with the moderate bending stiffness. While for the stiff plates, the thrust contributed by the normal force is not so significant due to the tiny deformation and the thrust is mainly induced by the leading-edge vortex. Therefore the scaling law based on $\overline{|\boldsymbol{F}^{n}|}$ is not so valid when $K^{\ast }$ is large. Hence, the rescaling works well only in the region of moderate $K^{\ast }$ .

4.7 Comparison with results in the literature

In the above, we concluded that under the circumstance of in-phase (IP) flapping, the SF mode is more efficient than the isolated case and the side-by-side case. However, for the antiphase (out-of-phase (OP)) flapping, the side-by-side is more efficient and there is no SF mode.

The optimal configurations for the two nearest neighbours in fish schooling are listed in table 2. It is seen that the SF mode in Katz et al. (Reference Katz, Tunstrøm, Ioannou, Huepe and Couzin2011) qualitatively agrees with our result about the SF mode. Besides, the SF cases have a very low distance $H<0.5$ (see figure 6) in terms of our simulation results. Movie S1 in Ashraf et al. (Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017) also shows that the distance is $H\approx 0.35$ , which is consistent with our results for the SF mode. It is also seen that both side-by-side and staggered-following modes are found in Partridge (Reference Partridge1982) and when the side-by-side configuration appears, the OP (out-of-phase) state is favoured (Ashraf et al. Reference Ashraf, Godoy-Diana, Halloy, Collignon and Thiria2016). Hence, some of our conclusions are basically consistent with the results in the literature.

Table 2.

Optimal configuration for two nearest neighbours. To quantify the orientation between a given fish and its nearest neighbour, the angle between the connected line of the two fish and the swimming direction is measured. Zero degrees is directly ahead of the fish; $180^{\circ }$ is directly behind.

It is also noticed that by using such a simple model it is impossible to fully reveal all hydrodynamic mechanisms of the fish schooling. There may exist some discrepancies between our results and those in the literature for the fish experiments. For example, the experimental study (Ashraf et al. Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017) showed that at the more energy-demanding, higher-velocity regimes, almost all schooling occurs in phalanx (side-by-side) formation. While in our simulations, not only the side-by-side formation but also the staggered formation (the SF mode) are found in the more energy-demanding, higher-velocity regimes. Another example is that Katz et al. (Reference Katz, Tunstrøm, Ioannou, Huepe and Couzin2011) and Ashraf et al. (Reference Ashraf, Bradshaw, Ha, Halloy, Godoy-Diana and Thiria2017) did not mention that the staggered arrangement leads to energy efficiency, while the SF mode in our study has higher efficiency. The inconsistency may be attributed to the limitation of our present model and the complexity of a natural fish school. As for the limitations of the present model, first our model consisting of two self-propulsion plates with simple actuation is a highly simplified model for a fish school. Second, it is only a two-dimensional study instead of three-dimensional. Third, the $Re$ is small with respect to the $Re$ regime of fish swimming.

It is worth mentioning that simulations/experiments from a limited number of bodies and at low $Re$ should be extended to a larger number and at high $Re$ with caution. In the present study, the wake interference of two self-propelled flapping plates at low $Re$ may not break the vortical structures into small structures and some coherent structures may still exist. At high $Re$ , although the wake of a single fish/object still shows coherent structures, wakes of many fish/object may become disorganized and the coherent structures may break down because of the interference of the wakes at high $Re$ (Daghooghi & Borazjani Reference Daghooghi and Borazjani2015).