2.1.1. Wing model and fluid model

As shown in figure 1, the pair of wings in tandem are each modelled as rigid plates of the same chord length $c$ . The up-and-down or heaving motions are prescribed by sinusoidal functions of the same frequency $f$ and (perhaps different) amplitudes $A_k$

(2.1) \begin{align} x_{k}(s,t) & = X_k(t)+s, \; s\in [-c/2,c/2], \end{align}

(2.2) \begin{align} y_{k}(t) & = A_k\cos (2\pi ft), \end{align}

where $k=1,2$ denotes the leader and follower wings, respectively, and $\boldsymbol{x}_k=(x_k,y_k)$ denotes their locations. Here, $y_k(t)$ is the same throughout the wing’s chord length (i.e. for any $s\in [-c/2,c/2]$ of each wing), and so it is a function of time only. The arc-length parameter $s$ denotes the signed distance from the wing centre $X_k(t)= x_k(0,t)$ . The leading edge of the wings (left end of the wing in figure 1) corresponds to $s=-c/2$ , and the trailing edge (right end) corresponds to $s=c/2$ .

The wings are individually free to swim or fly, that is, to move in the horizontal direction due to hydrodynamic interactions. Assuming a wing mass per unit span of $m$ , the horizontal dynamics of each individual satisfies

(2.3) \begin{align} m\ddot {X}_k = F_k = T_k+D_k,\quad k=1,2, \end{align}

where $F_k$ is the hydrodynamic force on the wing, composed of the hydrodynamic thrust $T_k$ and drag $D_k$ . Here, the propulsive thrust is negatively signed and the resistive drag is positively signed, given the convention that the wings swim or fly from right to left.

The surrounding inviscid fluid flow is described by the 2-D Euler equations

(2.4) \begin{align} \rho _{f}\frac {{\rm D}\boldsymbol{u}}{{\rm D}t}=-\boldsymbol{\nabla }p, \end{align}

where $\boldsymbol{u}=(u_x,u_y)$ is the fluid velocity, $p$ is the pressure and $\rho _{f}$ is the fluid density. We impose a no-penetration boundary condition for the flow velocity $\boldsymbol{u}(\boldsymbol{x}(s))$ on each wing, with the vertical component of the velocity being continuous and matching the wing flapping speed, in a reference frame where the fluid is at rest at infinity. That is, the vertical components $u_{y,\pm }$ above and below the wing are equal to the vertical component of the wing velocity

(2.5) \begin{eqnarray} u_{y,+}(\boldsymbol{x}_k(s,t))=u_{y,-}(\boldsymbol{x}_k(s,t))=\dot {y}_k(t),\quad k=1,2. \end{eqnarray}

In the tangential direction at the wing, i.e. the $x$ -direction, flow is allowed to slip freely along the surface. This horizontal component $u_{x,\pm }(\boldsymbol{x}_k(s,t))$ of the flow velocity can be viewed as a model for the flow just outside of the boundary layer, which approaches zero thickness in the limit of infinite Reynolds number. We later discuss a skin-friction model that estimates the drag due to shear in the viscous boundary layer.

Integration of the fluid pressure over a closed contour around each wing surface $C_k^w$ provides the thrust, i.e. the hydrodynamic force in the horizontal or tangential direction $\hat {\boldsymbol{s}}=(1,0)$ . This integral can be separated into parts associated with an infinitesimal circle (denoted by $r_k^{le}$ ) around the wing’s leading edge, the differential or jump in pressure $[p]$ across the length of the wing and an infinitesimal circle $r_k^{te}$ at the trailing edge

(2.6) \begin{align} \int _{C_k^w}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s \boldsymbol{\cdot }\hat {\boldsymbol{s}}=\left ( \int _{r_k^{le}}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s +\int _{-1}^1[p]\hat {\boldsymbol{n}}\,{\mathrm{d}} s + \int _{r_k^{te}}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s \right ) \boldsymbol{\cdot }\hat {\boldsymbol{s}}, \quad \hat {\boldsymbol{s}}=(1,0). \end{align}

In our model, we allow continuous shedding of vorticity at the trailing edge while keeping a flow singularity at the leading edge. The absence of leading-edge shedding is intended as a model applicable to the case of small-amplitude flapping and long wavelengths of motion, in which case the small angles of attack are associated with weak leading-edge separation. The singularity at the leading edge yields a leading-edge suction and thrust in the tangential direction to the wing (Saffman Reference Saffman1993). Note that the integral of the pressure jump $[p]$ along the wing is always normal to the wing and does not contribute to the thrust. The pressure integral around the wing trailing edge is also zero due to the imposition of the unsteady Kutta condition. Therefore, for a heaving wing that swims in a direction parallel to the body axis, the thrust results purely from leading-edge suction. Thus, (2.6) simplifies to

(2.7) \begin{align} T_k =\int _{C_k^w}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s \boldsymbol{\cdot }\hat {\boldsymbol{s}} =\int _{r_k^{le}}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s \boldsymbol{\cdot }\hat {\boldsymbol{s}}, \quad \hat {\boldsymbol{s}}=(1,0). \end{align}

The hydrodynamic drag in (2.3) results from skin friction along the upper and lower wing surfaces, which is modelled using Blasius laminar boundary-layer theory (Schlichting Reference Schlichting1968, Chap. VII). For one side of a plate in a stream of uniform far-field speed $U_f$ , this approximation gives

(2.8) \begin{align} D_f=0.664\sqrt {\rho _{f}\mu c}\,| U_f|^{3/2}, \end{align}

where $\mu$ is the dynamic viscosity of the fluid. In our model, we replace the free-stream velocity $U_f$ with the tangential flow speed averaged along the wing surface and relative to the swimming speed

(2.9) \begin{align} \overline {U}_{k,\pm }(t) = \frac {1}{c}\int _{-c/2}^{c/2}u_{x,\pm }(\boldsymbol{x}_k(s,t))\,{\mathrm{d}} s-U_k(t), \end{align}

where $U_k(t)=\dot {X}_{k}(t)$ is the horizontal swimming speed of the wing, which is assumed to move in the negative $x$ -direction (figure 1). Summing the contributions from the upper and lower surfaces of the wing, the total drag is

(2.10) \begin{align} D_k=0.664\sqrt {\rho _{f}\mu c}\,\left (|\overline {U}_{k,+}|^{3/2}+|\overline {U}_{k,-}|^{3/2}\right ). \end{align}

This form of the drag force was also implemented by Heydari & Kanso (Reference Heydari and Kanso2021). We note that the drag force scaling with $U^{3/2}$ reproduces the experimentally observed $(Af)^{4/3}$ dependence of the swimming speed for a single heaving wing (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016). This is consistent with the scaling identified by Gazzola, Argentina & Mahadevan (Reference Gazzola, Argentina and Mahadevan2014) for locomotion in the laminar-flow regime at Reynolds numbers up to approximately $10^4$ and where drag dominantly arises from skin friction due to laminar boundary-layer flows. At yet higher Reynolds numbers of the turbulent-flow regime, a quadratic drag law of the form $U^2$ is likely more appropriate and derives from pressure and separated flows (Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014).

We non-dimensionalise the system with the characteristic length $\tilde {L}=c/2,$ time $\tilde {T}=1/f$ , velocity $\tilde {U}=cf/2$ and pressure $P=\rho _f\tilde {U}^2=\rho _{f}c^{2}f^{2}/4.$ The dimensionless equations of motion are then

(2.11) \begin{align} M\ddot {X}_k & = T_k(t)+D_k(t), \end{align}

(2.12) \begin{eqnarray} y_{k}(s,t) & = \tilde {A}_k\cos (2\pi t),\ s\in [-1,1], \end{eqnarray}

where

(2.13) \begin{align} T_k = \int _{C_k^w}p\hat {\boldsymbol{n}}\,{\mathrm{d}} s\boldsymbol{\cdot }\hat {\boldsymbol{s}}\quad \text{and}\quad D_k = C_f\left (\overline {U}_{k,+}^{3/2}+\overline {U}_{k,-}^{3/2}\right ),\; \quad k=1,2. \end{align}

There are three dimensionless parameters that appear in the system: the wing–fluid mass ratio $M=4m/(\rho _fc^2)$ , the ratio of the heaving amplitude to wing length $\tilde {A}_k = 2A_k/c$ and the skin-friction drag coefficient

(2.14) \begin{align} C_f = 1.878\sqrt {\nu /(c^2f)}, \end{align}

where $\nu =\mu /\rho _f$ is the kinematic viscosity. This drag coefficient is related to the so-called frequency Reynolds number $Re_{fr}=\tilde {L}\tilde {U}/\nu =c^2f/4\nu$ used in the literature of flapping-wing locomotion (Alben & Shelley Reference Alben and Shelley2005): $C_f=0.939Re_{fr}^{-1/2}$ . When the kinematics, wing chord length and fluid parameters are appropriately matched to the experimental conditions of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), (2.14) yields $C_{f}\sim 0.02$ , reproducing the same trends in swimming speed observed in the experiments. In our study, we set the mass ratio to $M=2$ while exploring the effects of varying $C_{f}=0.01$ to $0.1$ and $\tilde {A}_{k}=0.1$ to $1$ . Here, the mass per unit span $m$ should be interpreted as the total solid mass of the locomotor (i.e. the free-swimming or flying animal), which sets the relevant inertia for the forward dynamics. We note that, even for an arbitrarily thin wing, a finite value of wing–fluid mass ratio $M$ remains meaningful because it represents the mass payload (e.g. the body of a fish or bird) that is towed by the propulsor (fin or wing). The choice of $M=2$ is not meant to match a specific biological system; rather, it lies within the range typical of swimming and flying animals (Greenewalt Reference Greenewalt1962; Rayner Reference Rayner1979; Viscor & Fuster Reference Viscor and Fuster1987; Dimitriadis et al. Reference Dimitriadis, Gardiner, Tickle, Codd and Nudds2015; Taylor et al. Reference Taylor, Taylor, Lambert, Walker, Biro and Portugal2019; Krishnan et al. Reference Krishnan2022) and is consistent with flapping-foil experiments and prior numerical studies (Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016; Newbolt et al. Reference Newbolt, Zhang and Ristroph2019, Reference Newbolt, Lewis, Bleu, Wu, Mavroyiakoumou, Ramananarivo and Ristroph2024; Nitsche et al. Reference Nitsche, Oza and Siegel2025). A smaller mass, and thus lower inertia, is expected to make the wing more susceptible to fluctuations. This leads to larger oscillations in the separation gap between the wings, which could undermine the positional stabilisation of the wing. The influence of mass will be examined further in future work.

2.1.2. Vortex-sheet shedding and wake model

In the vortex-sheet model, a flapping wing is a bound vortex sheet on which the no-penetration boundary condition (2.5) is imposed. The bound sheet encompasses the infinitesimally thin wing and its boundary layers that vanish in thickness at higher Reynolds numbers. At the trailing edge of each wing, a free vortex sheet is continuously shed into the fluid (Nitsche & Krasny Reference Nitsche and Krasny1994). In the free-sheet model the shed shear layer forms downstream eddies. In the limit of zero viscosity, the shed vorticity concentrates onto an infinitesimally thin layer that can be modelled as a 1-D vortex sheet that does not dissipate in time. The vortex shedding rate at the trailing edge is determined by the unsteady Kutta condition (Jones Reference Jones2003), which provides the direction of vortex shedding as well as the amount of circulation transmitted from the wing (bound sheet) to the wake (free sheet).

In our 2-D vortex-sheet method, we keep track of the sheet position $\zeta _k(s,t)$ as well as the true vortex-sheet strength $\bar {\gamma }_k(s,t)$ , defined as the tangential fluid velocity jump across the vortex sheet $\bar {\gamma }_k(s,t)=[\boldsymbol{u}]\hat {s}$ (Shelley Reference Shelley1992). Here, $k=1,2$ denotes the two wings, and $s\in [-1,s^k_{max }]$ is an arc-length parameterisation of the vortex sheet. The true vortex-sheet strength $\bar {\gamma }_k(s,t)=\gamma _k(\alpha )/s_{\alpha }(\alpha ,t)$ is related to the (unnormalised) vortex strength $\gamma _k(\alpha )$ which is a material variable, where $\alpha$ is a Lagrangian parameterisation of the sheet. Using complex variables, the vortex-sheet position $\zeta _k(s,t)$ and strength $\bar {\gamma }_k(s,t)$ are related through the Biot–Savart law for the mean velocity of the sheet $w_k(s,t)$

(2.15) \begin{align} \overline {w_k(s,t)} = \frac {1}{2\pi i}-\!\!\!\!\!\kern-1pt\int_{-1}^{{s}_{max}^{k}}\frac {\bar {\gamma }_{k}(s',t)\,{\mathrm{d}} s'}{\zeta _{k}(s,t)-\zeta _{k}(s',t)} +\sum _{j\neq k}\frac {1}{2\pi i}\int _{-1}^{s_{max }^j}\frac {\bar {\gamma }_{j}(s',t)\,{\mathrm{d}} s'}{\zeta _{k}(s,t)-\zeta _{j}(s',t)}, \end{align}

where the bar denotes the complex conjugate, and $ -\kern-1.5pt\!\!\!\!\int $ is the Cauchy principal value. Since the tangential component of fluid velocity is discontinuous at the vortex sheet, the mean velocity $w_k(s,t)$ is an average of the velocities on the two sides of the sheet boundary $\zeta _k(s,t)$ .

The description of the free-sheet boundary $s\in [1,s^k_{max }]$ dynamics is simplified if $w_{k}(s,t)$ is chosen as a velocity frame, in which case we arrive at the Birkhoff–Rott equation (Shelley Reference Shelley1992)

(2.16) \begin{align} D_t\zeta _{k}(s,t) =w_{k}(s,t),\quad 1\leqslant s\leqslant s_{max }^k, \end{align}

where $D_t$ denotes the material derivative, and $w_{k}(s,t)$ is given by (2.15). Following the Lagrangian path of the free vortex sheet, the circulation $\gamma (\alpha )$ on an infinitesimal segment is conserved

(2.17) \begin{align} D_t(\bar {\gamma }_{k}(s,t) s_{\alpha })=D_t\gamma _{k}(\alpha )=0,\quad k=1,2. \end{align}

The bound vortex-sheet position, i.e. the wing position

(2.18) \begin{align} \zeta _{k}(s,t) = X_k(t)+s+iy_k(t),\quad s\in [-1,1], \end{align}

is governed by (2.11). Using (2.15), the no-penetration boundary condition of (2.5) yields

(2.19) \begin{align} \textrm {Im}\left (w_{k}(s,t)\right )=\dot {y}_k(t),\quad s\in [-1,1], \end{align}

which is a Fredholm integral equation of the first kind for the bound-sheet strength $\bar {\gamma }_k(s,t),\ s\in [-1,1]$ .

The unsteady Kutta condition connects the bound and free sheets at the trailing edge (Jones Reference Jones2003). It requires the vortex shedding to be tangential, i.e. the interface of the bound and free vortex sheets is smooth, and also determines the amount of circulation shed from the bound sheet into the free sheet

(2.20) \begin{align} \dot {\varGamma }_k(t)+\left [\textrm {Re}(w_k(1,t))-U_k(t)\right ]\bar {\gamma }_k(1,t)=0, \end{align}

where $\varGamma _k(t) = \int _{-1}^1\bar {\gamma }_k(s,t)\,{\mathrm{d}} s$ denotes the total circulation on the bound sheet, and $U_k(t)=\dot {X}_{k}(t)$ is the wing horizontal swimming speed. Kelvin’s circulation conservation theorem ensures that $\varGamma _k(t) = -\int _1^{s_{max }^k}\bar {\gamma }_k(s,t)\,{\mathrm{d}} s$ .

An equation for the pressure jump $[p]_{k}(s,t)$ distributed on the bound sheet can be derived from the Euler equation (2.4) (Saffman Reference Saffman1993; Jones Reference Jones2003; Alben Reference Alben2009b ) and the Kutta condition (2.20)

(2.21) \begin{align} [p]_{k}(s,t)=\int _{1}^{s}\partial _{t}\bar {\gamma }_{k}(s',t)\,{\mathrm{d}} s'+\left [\textrm {Re}(w_k(s,t))-U_k(t)\right ]\bar {\gamma }_{k}(s,t)+\dot {\varGamma }_{k}(t). \end{align}

Note that the pressure jump across the wing contributes no force in the wing swimming direction, as the wing maintains a horizontal alignment at all times. The thrust arises from the suction force at the wing leading edge (Saffman Reference Saffman1993)

(2.22) \begin{align} T_k(t) = \frac {\pi }{8}\nu _k^{2}(-1,t), \end{align}

where $\nu _k(s,t)=\bar {\gamma }_k(s,t)\sqrt {1-s^2}$ has a finite value at $s=-1$ as $\bar {\gamma }_k$ has an inverse square root singularity at the wing leading edge. A full derivation of the leading-edge suction term can be found in Appendix A of Nitsche et al. (Reference Nitsche, Oza and Siegel2025). Note that (2.22) relates the hydrodynamics with the wing dynamics through (2.11).

On each wing, the instantaneous input power $P_{k,{in}}$ is defined as the product of flapping speed $\dot {y}_k$ and the pressure jump $\int _{-1}^{1}[p]_k\,{\mathrm{d}} s$ , which reflects the energy required to maintain the prescribed flapping motion. The instantaneous output power $P_{k,{out}}$ equals the rate of work done by the hydrodynamic thrust $T_k$ on the wing at swimming speed $U_k$ . Hence,

(2.23) \begin{align} P_{k,{in}}(t)&=-\dot {y}_k(t)\int _{-1}^{1}[p]_k(s,t)\,{\mathrm{d}} s, \end{align}

(2.24) \begin{align} P_{k,{out}}(t)=T_k(t)U_k(t)=\frac {\pi }{8}\nu _k^{2}(-1,t)U_k(t). \end{align}

The ratio of $P_{k,{in}}$ and $P_{k,{out}}$ gives the Froude efficiency on each wing $\eta _k,\ k=1,2$ and also the efficiency for the pair $\eta$ (Lighthill Reference Lighthill1960; Wu Reference Wu1961)

(2.25) \begin{align} \eta _k=\frac {\langle P_{k,{out}}\rangle }{\langle P_{k,{in}}\rangle },\quad \eta =\frac {\langle P_{1,{out}}\rangle +\langle P_{2,{out}}\rangle }{\langle P_{1,{in}}\rangle +\langle P_{2,{in}}\rangle }, \end{align}

where $\langle \boldsymbol{\cdot }\rangle =\int _t^{t+1}(\boldsymbol{\cdot }) \,{\mathrm{d}} t'$ denotes the time average over a flapping stroke.

In the end, we express the mean tangential flow speed along the wing $\overline {U}_{k,\pm }(t)$ (2.9) using vortex-sheet variables in order to calculate the skin-friction drag $D_k$ of (2.13). By the definitions of the mean velocity $w_k(s,t)$ and vortex-sheet strength $\bar {\gamma }_k(s,t)$ , the tangential speed along the two surfaces of each wing $u_{x,\pm }(\zeta _k(s,t))$ satisfies

(2.26) \begin{align} u_{x,+}(\zeta _k(s,t))+u_{x,-}(\zeta _k(s,t)) = 2\textrm {Re}(w_k(s,t)),\nonumber \\ u_{x,+}(\zeta _k(s,t))-u_{x,-}(\zeta _k(s,t)) = \bar {\gamma }_k(s,t), \end{align}

from which we obtain

(2.27) \begin{align} u_{x,\pm }(\zeta _k(s,t)) = \textrm {Re}(w_k(s,t))\pm \frac {1}{2}\bar {\gamma }_k(s,t). \end{align}

Substituting (2.27) into (2.9), $\overline {U}_{k,\pm }(t)$ is then expressed as

(2.28) \begin{align} \overline {U}_{k,\pm }(t)=\frac {1}{2}\int _{-1}^{1}\textrm {Re}(w_k(s,t))\,{\mathrm{d}} s-U_k(t)\pm \frac {1}{4}\int _{-1}^{1}\bar {\gamma }_k(s,t)\,{\mathrm{d}} s. \end{align}

4.1.1. Swimming of a single flapping wing

We begin by considering a single wing swimming with constant speed $U_{\infty }$ in a quiescent background. We use the same non-dimensionalisation as in § 2: the characteristic length $\tilde {L}=c/2,$ time $\tilde {T}=1/f,$ velocity $\tilde {U}=cf/2$ and pressure $P=\rho _f\tilde {U}^2=\rho _{f}c^{2}f^{2}/4$ . The heaving motion of the wing described in the body frame is

(4.1) \begin{align} h(x,t)=\tilde {A}\cos (2\pi t),\quad -1\leqslant x\leqslant 1, \end{align}

where the wing spans from $x=-1$ to $x=1$ .

In the frame of the flow, we assume the wing translates with constant velocity $-U_{\infty }$ in the negative $x$ -direction. Then in the wing frame, the wing is considered to be flapping in a uniform background stream of constant velocity $U_{\infty }$ in the positive $x$ -direction. Following Wu (Reference Wu1961), the flow velocity is denoted as $\boldsymbol{q}=(U_{\infty }+u,v),$ where velocity perturbations $(u,v)$ are assumed small compared with the uniform stream $U_{\infty }$ , i.e. $U_{\infty }\gg u$ and $U_{\infty }\gg v$ . The flow satisfies the 2-D Euler equation, and is linearised as

(4.2) \begin{align} \partial _{t}\boldsymbol{q}+U_{\infty }\partial _x \boldsymbol{q} = -\boldsymbol{\nabla }p,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} = 0. \end{align}

The kinematic boundary condition on the wing requires the velocity to be continuous in the normal direction. Assuming the heaving amplitude is small, $\tilde {A} \ll 1$ , the boundary condition is then linearised as

(4.3) \begin{align} v(x,y,t)\mid _{y=\pm 0}=\partial _t h(x,t) ,\quad -1\leqslant x\leqslant 1. \end{align}

Note that in dimensional quantities, (4.3) and condition $U_{\infty }\gg v$ imply that the Strouhal number ${\textit {St}}=Af/U_{\infty }$ is small, as $U_{\infty } \gg v\sim Af$ , where the Strouhal number denotes the ratio of the vertical flapping speed and the translational speed of the wing.

Instead of solving the fluid velocity directly, Wu (Reference Wu1961) solves for acceleration potential

(4.4) \begin{align} w(z,t) = \phi (x,y,t) + i\psi (x,y,t), \end{align}

in complex variables $z=x+iy$ . Here,

(4.5) \begin{align} \phi (x,y,t) = -p \end{align}

is Prandtl’s acceleration potential, which is a harmonic function of $(x,y)$ for every $t$ , and $\psi (x,y,t)$ is its harmonic conjugate. The pressure $p$ is continuous everywhere except across the wing boundary, $-1 \leqslant x \leqslant 1$ . It then follows that $\phi$ is regular inside the flow, so that

(4.6) \begin{align} \phi (x,+0,t)=\phi (x,-0,t), \quad \text{for} \; |x|\gt 1, \end{align}

and for all $t$ . The Kutta condition is required at the trailing edge of the wing

(4.7) \begin{align} |w(1,t)| \lt \infty , \quad \forall t \geqslant 0. \end{align}

Furthermore, far-field boundary conditions

(4.8) \begin{align} w(z,t)\rightarrow 0 \quad \text{as}\quad |z|\rightarrow \infty ; \quad u(z,t)-iv(z,t)\rightarrow 0 \quad \text{as}\quad z\rightarrow -\infty , \end{align}

are required to complete the problem.

Wu (Reference Wu1961) constructed the solution of $w(z,t)$ using Fourier series and conformal transformation. Adjusting this to the pure-heaving plate case considered in this work, simplifies the Fourier series expressions given in Wu (Reference Wu1961). Specifically, all higher Fourier modes vanish, leaving only the uniform mode, and any terms involving $\partial _x h$ are identically zero. Therefore, provided by the solution, the stroke-averaged hydrodynamic thrust on the wing is given as

(4.9) \begin{align} \langle T \rangle = -4\pi ^{3}\tilde {A}^{2}\left |C(\sigma )\right |^2, \end{align}

(4.10) \begin{align} C(\sigma )=\frac {K_{1}(i\sigma )}{K_{0}(i\sigma )+K_{1}(i\sigma )}=C_{1}(\sigma )+iC_{2}(\sigma ),\quad \sigma =2\pi /U_{\infty }, \end{align}

where $K_0$ and $K_1$ are the modified Bessel functions of the second kind of orders zero and one. The function $C(\sigma )$ is called Theodorsen’s function, and $C_1(\sigma )$ and $C_2(\sigma )$ are the real and imaginary parts of $C(\sigma )$ , respectively. Using dimensional quantities, $\sigma = \pi cf/U_{\infty }$ is usually called the reduced frequency. If we consider the problem in the flow frame, that the wing moves with velocity $U_{\infty }$ in a quiescent background, $\lambda =U_{\infty }/f$ denotes the wavelength of the wing motion and then $\sigma = \pi c/\lambda$ represents the ratio of the wing size to the wavelength. The thrust (4.9) is in the same direction as swimming (negative $x$ -direction) and is composed of the leading-edge suction only. As discussed in § 2.1 regarding the vortex-sheet model, the pressure distributed along the wing has zero horizontal component, since the body is always horizontally aligned with the flow.

The Froude efficiency (Lighthill Reference Lighthill1960), or the stroke-averaged efficiency of propulsion, is defined in the same way as in the vortex-sheet model (§ 2.1), which is provided by Wu (Reference Wu1961) as

(4.11) \begin{align} \eta = \frac {\langle P_{{out}} \rangle }{ \langle P_{{in}} \rangle } = \frac {\left |C(\sigma )\right |^2 }{C_{1}(\sigma )}. \end{align}

Here, the output power

(4.12) \begin{align} \langle P_{{out}}(t) \rangle =\langle T(t) \rangle U_{\infty } \end{align}

is computed from (4.9). The input power required to maintain the flapping motion is expressed as

(4.13) \begin{align} \langle P_{{in}} \rangle = \biggl \langle -\int _{-1}^{1}[p] \partial _t h \,{\mathrm{d}} x \biggr \rangle ; \end{align}

calculated by integrating the instantaneous pressure jump distribution $[p](x,t)$ along the wing, where

(4.14) \begin{align} [p](x,t)=p(x,+0,t)-p(x,-0,t),\quad -1\leqslant x\leqslant 1, \end{align}

computed from Wu’s solution of $w(z,t)$ .

Note that using dimensional quantities, the thrust formula (4.9) is expressed as

(4.15) \begin{align} \langle T \rangle = -2\pi ^{3}c \rho _{f} A^{2}f^{2}\left |C(\sigma )\right |^2, \quad \sigma =\pi cf/U_{\infty }, \end{align}

which is quadratic in flapping velocity $V_f=2\pi Af$ for fixed reduced frequency $\sigma$ . The limiting case of $\sigma \rightarrow 0$ is of particular interest. For small $\sigma$ , $K_0(i \sigma )$ and $K_1(i\sigma )$ can be expanded as

(4.16) \begin{align} K_{0}(i\sigma ) = -\ln \sigma +\ln 2-\gamma -i\pi /2 + O(\sigma ^{2}\ln \sigma ), \end{align}

(4.17) \begin{align} K_{1}(i\sigma ) =-\frac {\pi }{4}\sigma -\frac {i}{\sigma }+\frac {i\sigma }{4}\left (2\ln \sigma -2\ln 2+2\gamma -1\right )+ O(\sigma ^{3}\ln \sigma ), \end{align}

where $\gamma \approx 0.5771$ is Euler’s constant (Abramowitz & Stegun Reference Abramowitz and Stegun1964). Substituting (4.16)–(4.17) in Theodorsen’s function (4.10), we obtain $C(\sigma )\rightarrow 1$ , $C_1(\sigma )\rightarrow 1$ as $\sigma \rightarrow 0$ . Hence $\langle T \rangle \sim -2\pi ^{3}c \rho _{f} A^{2}f^{2}$ for small $\sigma$ , or long free-swimming wavelength $\lambda /c=\pi /\sigma$ . The Froude efficiency $\eta$ in (4.11) is determined only by the reduced frequency $\sigma$ , and has limit $\eta \rightarrow 1$ as $\sigma \rightarrow 0$ (see figure 13 b). In the other limiting case of $\sigma \rightarrow +\infty ,$ we could obtain $C(\sigma )=C_1(\sigma )\rightarrow 1/2$ , and thus $\langle T \rangle \sim -\pi ^{3}c \rho _{f} A^{2}f^{2}$ and $\eta \rightarrow 1/2$ for $\sigma \gg 1$ . The limit at large $\sigma$ , or small free-swimming wavelength $\lambda /c=\pi /\sigma$ , is not of much interest to free-swimming problems, as the free-swimming speed $U_{\infty } = \pi cf/\sigma$ , determined by the flapping amplitude $A$ and frequency $f$ , does not approach a small value for the parameter ranges of our interest.

Figure 13.

Comparison of linear theory (dashed lines) with simulation (symbols) of single self-propelled flapping wing, for various flapping amplitudes $\tilde {A}=2A/c$ and drag coefficients $C_f$ . Simulation data are the same as in figure 6. (a) Steady swimming speed $|U_{\infty }|$ , normalised by the characteristic speed $cf/2$ . The power law of steady speed $2|U_{\infty }|/(cf)$ in the amplitude $2A/c$ is between linear and quadratic, as shown in the log–log plot in the inset. (b) The Froude efficiency $\eta$ of the single wing at steady state. The linearised theory (4.11) predicts limits of $\eta \rightarrow 1$ as $\lambda /c\rightarrow +\infty$ , and $\eta \rightarrow 1/2$ as $\lambda /c\rightarrow 0$ .

For the freely swimming wing, the vertical flapping kinematics determines the free translational speed in the horizontal direction, that is, the free speed $U_{\infty }$ is a function of the flapping amplitude $A$ and frequency $f$ . The steady free speed $U_{\infty }$ can be solved from balancing the thrust and drag. The thrust $\langle T \rangle$ is provided by (4.9) (dimensionless), or (4.15) (dimensional). The drag $D$ is modelled using a simplification of the skin-friction model (2.13) used in simulations, where the boundary-layer velocity is approximated using the free-swimming velocity of the wing

(4.18) \begin{align} D=2C_f|U_{\infty }|^{3/2}, \end{align}

and the drag coefficient $C_f$ is provided by (2.14). The dimensionless free-swimming speed $U_{\infty }$ is then solved from balance of forces

(4.19) \begin{align} 0=\langle T \rangle + D = -4\pi ^{3}\tilde {A}^{2}\left |C(\sigma )\right |^2+2C_f|U_{\infty }|^{3/2},\quad \sigma =2\pi /U_{\infty }. \end{align}

The steady swimming speed $U_{\infty }$ , solved from (4.19), is shown in figure 13(a), compared with the steady speed of a single flapping wing calculated through the vortex-sheet simulations (see figure 6 a). Figure 13(b) compares the Froude efficiency of the single wing at steady swimming state, calculated using vortex-sheet simulation and the linearised theory (4.11). The results from linearised theory agree well with the simulation results. For the parameter values used in our simulations, the flapping amplitude ( $A/c = 0.05$ to $0.5$ ) and the Strouhal number ( ${\textit {St}}\approx 0.005$ to $0.1$ ) are within small ranges, so the linearisation provides a good approximation. The log–log plot of the free speed, in the inset of figure 13(a), shows that the free speed $U_{\infty }$ is superlinear in flapping amplitude, which agrees with the power law in experiments. This will be discussed further in § 5.

4.1.2. Fluid wake generated by single flapping wing

The vortex wake generated by a flapping wing is complex and self-developing, as shown by our simulation figure 11 and by experimental flow visualisations (Vandenberghe et al. Reference Vandenberghe, Childress and Zhang2006; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016). In the linear regime, the Strouhal number $St=Af /U_{\infty }$ , the ratio of the vortex wake width $A$ and wavelength $\lambda =U_{\infty }/f$ , is assumed to be small. The vortex sheet, shed tangentially from the trailing edge (by the Kutta condition; see Jones (Reference Jones2003)), is assumed to lie along the $x$ -axis downstream of the wing trailing edge (as in the schematic diagram in figure 14). We define the vortex-sheet strength $\gamma (x,t)$ at point $(x,0)$ as

(4.20) \begin{align} \gamma (x,t)=u(x,+0,t) - u(x,-0,t). \end{align}

We consider the vortex sheet in the wing frame; the sheet defined on $x\geqslant 1$ downstream the wing trailing edge ( $x=1$ ). According to (4.2) and (4.6), the vortex-sheet strength $\gamma$ satisfies

(4.21) \begin{align} \partial _t\gamma + U_{\infty }\partial _x\gamma = 0,\quad x\geqslant 1, \end{align}

which shows that $\gamma$ is dependent only on a single variable $(x-U_{\infty } t)$ , and

(4.22) \begin{align} \gamma (x,t) = \gamma \left (1,t-\frac {(x-1)}{U_{\infty }}\right ). \end{align}

The vortex-sheet strength at the trailing edge is calculated by Wu (Reference Wu1961) as

(4.23) \begin{align} \gamma (1,t)=\frac {4\pi ^{2}\tilde {A}}{G(\sigma )}\sin \left (2\pi \left (t-\frac {\sigma }{2\pi }-g(\sigma )\right )\right ),\quad \sigma =2\pi /U_{\infty }, \end{align}

(4.24) \begin{align} \mathcal{G}(\sigma )=K_{0}(i\sigma )+K_{1}(i\sigma )=G(\sigma )e^{2\pi ig(\sigma )}, \end{align}

where $G(\sigma )$ and $g(\sigma )$ are the modulus and angle of $\mathcal{G}(\sigma )$ , respectively. Substituting (4.23) in (4.22), the vortex-sheet strength is then provided as

(4.25) \begin{align} \gamma (x,t)=\frac {4\pi ^{2}\tilde {A}}{G(\sigma )}\sin \left (2\pi \left (t-\frac {\sigma x}{2\pi }-g(\sigma )\right )\right ),\quad \sigma =2\pi /U_{\infty }, \end{align}

where $x\geqslant 1$ denotes the downstream of the wing trailing edge.

Figure 14.

The vortex sheet generated by a single flapping wing with swimming velocity $-U_{\infty }$ (in the negative $x$ -direction) is assumed to lie along the $x$ -axis downstream of the wing. It has $\delta _1$ phase lag from the wing trailing-edge trajectory (dashed line). The wake induced by the flat vortex sheet can be approximated by a steady stationary wave $V_w$ , which has $\delta _2$ phase lag from the wing trailing-edge trajectory. See (4.29) and (4.37). The red denote positive vorticity and the blue $\ominus$ denote negative vorticity.

Now we change back to the flow frame, where the wing swims with speed $-U_{\infty }$ towards the negative $x$ -direction in a quiescent environment. The flow frame variables $(\tilde {x},\tilde {y},\tilde {t})$ and the wing frame variables $(x,y,t)$ are related through

(4.26) \begin{align} (\tilde {x},\tilde {y},\tilde {t})=(x-U_{\infty }t,y,t). \end{align}

The position of the wing trailing edge $(\tilde {x}_b,\tilde {y}_b)$ in the flow frame can be parameterised in $\tilde {t}$ as

(4.27) \begin{align} \tilde {x}_b = -U_{\infty } \tilde {t}+1,\quad \tilde {y}_b = \tilde {A}\cos \left (2\pi \tilde {t}\right ), \end{align}

from which the trajectory of the trailing edge is thus provided as

(4.28) \begin{align} \tilde {y}_b(\tilde {x})=\tilde {A}\cos \left (2\pi \frac {(1-\tilde {x})}{U_{\infty }}\right )=\tilde {A}\cos \left (2\pi \frac {(\tilde {x}-1)}{\lambda }\right ). \end{align}

Here, the wavelength of the motion is $\lambda =U_{\infty }=2\pi /\sigma$ using dimensionless variables. Substituting (4.26) in (4.25), we obtain that the vortex-sheet strength in the flow frame is given by

(4.29) \begin{align} \gamma (\tilde {x})=-\frac {4\pi ^{2}\tilde {A}}{G(\sigma )}\sin \left (2\pi \left (\frac {\tilde {x}}{\lambda }+g(\sigma )\right )\right ),\quad \tilde {x}\geqslant \tilde {x}_b, \end{align}

where $\tilde {x}_b=-U_{\infty } \tilde {t}+1$ denotes the location of the wing trailing edge.

We note that the vortex-sheet strength (4.29) has a phase lag of $\delta _1$ from the trailing-edge trajectory (4.28), where

(4.30) \begin{align} \delta _1(\sigma ) = -1/4-g(\sigma )-\sigma /(2\pi ). \end{align}

In the limiting case of $\sigma \rightarrow 0$ , we obtain

(4.31) \begin{align} G(\sigma )\rightarrow \infty ,\; g(\sigma )\rightarrow -1/4, \end{align}

using (4.16)–(4.17) and (4.24). Substituting (4.31) in (4.29)–(4.30), it then follows that the magnitude of the vortex-sheet strength $\gamma _0 = 4\pi ^2 \tilde {A}/G(\sigma )\rightarrow 0$ , and the phase lag $\delta _1(\sigma )\rightarrow 0$ as $\sigma \rightarrow 0$ . These limits indicate a weak vortex-sheet wake when the wavelength $\lambda$ , or the swimming speed $U_{\infty }$ is large. Moreover, the vorticity shed by the trailing edge is nearly in phase with the motion of the trailing edge, in the regime of large $\lambda$ or $U_{\infty }$ . This agrees with experimental flow visualisations in Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), which show that vortex cores are roughly located at the peak of the trailing-edge trajectory. The limit of $\sigma \rightarrow +\infty$ is of less interest as the steady speed of free swimming $U_{\infty }$ does not achieve small values for the parameter values of our interest. For the range $U_{\infty }\in [2,\infty )$ (i.e. $\sigma \in (0,\pi ]$ ), which covers the range of $U_{\infty }$ in our simulations, we compute numerical ranges $0 \lt 1/G(\sigma ) \lesssim 0.705$ , and $-0.12 \lesssim \delta _1(\sigma ) \lt 0$ which indicates that the phase-lag $\delta _1$ between the trailing-edge trajectory and the vorticity is always less than a quarter.

The fluid velocity induced by the vortex sheet (4.29) can be calculated through the Biot–Savart law. To calculate the flow velocity on the $x$ -axis downstream the wing (or the mean velocity at the vortex sheet), we ignore the circulation on the wing body, and approximate the flow field using that induced by the vortex sheet, provided as

(4.32) \begin{align} u(\tilde {x},0,t)-iv(\tilde {x},0,t)=\frac {1}{2\pi i}-\!\!\!\!\!\kern-1pt\int_{\tilde{x}_b}^{\infty} \frac {\gamma (\tilde {x}')}{\tilde {x}-\tilde {x}'}\,{\mathrm{d}}\tilde {x}', \quad \tilde {x}\geqslant \tilde {x}_b, \end{align}

where $\tilde {x}_b$ denotes the location of the wing trailing edge. Note that $u(\tilde {x},0,t)=0$ , as the right-hand side of (4.32) is purely imaginary. While the horizontal velocity on the $x$ -axis, downstream the wing, is zero, the vertical velocity is given by

(4.33) \begin{align} v(\tilde {x},0,t) & = \frac {1}{2\pi }-\!\!\!\!\!\kern-1pt\int_{\tilde{x}_b}^{+\infty} \frac {\gamma (\tilde {x}')}{\tilde {x}-\tilde {x}'}\,{\mathrm{d}}\tilde {x}'\nonumber \\ & = -\frac {2\pi \tilde {A}}{G(\sigma )}-\!\!\!\!\!\kern-1pt\int_{\tilde{x}_b}^{+\infty} \frac {\sin \left (2\pi (\tilde {x}/\lambda +g(\sigma ))\right )}{\tilde {x}-\tilde {x}'}\,{\mathrm{d}}\tilde {x}'\nonumber \\ & = \frac {2\pi ^2 \tilde {A}}{G(\sigma )}\left \{ \cos \left (2\pi \left (\frac {\tilde {x}}{\lambda }+g(\sigma )\right )\right )+J(\tilde {x},\tilde {x}_b)\right \}, \end{align}

(4.34) \begin{align} J(\tilde {x},\tilde {x}_b) & = \frac {1}{\pi }\left [\left (\mathrm{Si}(2\pi (\tilde {x}-\tilde {x}_b)/\lambda ) -\dfrac {1}{2}\pi \right )\cos \left (2\pi (\tilde {x}/\lambda +g(\sigma ))\right )\right .\nonumber \\ &\quad\left . -\mathrm{Ci}(2\pi (\tilde {x}-\tilde {x}_b)/\lambda )\sin \left (2\pi (\tilde {x}/\lambda +g(\sigma ))\right )\right ], \quad \tilde {x}\geqslant \tilde {x}_b, \end{align}

where $\mathrm{Si}(x)=\int _0^x(\sin t/t)\,{\mathrm{d}} t$ and $\mathrm{Ci}(x)=-\int _x^{+\infty }(\cos t/t)\,{\mathrm{d}} t$ are the sine integral and cosine integral of real arguments, respectively.

Note that (4.34) can be bounded by its magnitude function $I(\tilde {x}-\tilde {x}_b)$ , where

(4.35) \begin{align} I(x) = \sqrt {\left (\mathrm{Si}(2\pi x)-\dfrac {1}{2}\pi \right )^{2}+\mathrm{Ci}^{2}(2\pi x)}. \end{align}

It is easy to show that $I(x)$ is a monotonically decreasing function. For location $\tilde {x}$ that is at least one wavelength away from the wing trailing edge $\tilde {x}_b$ , i.e. $(\tilde {x}-\tilde {x}_b)/\lambda \geqslant 1$ , (4.34) is bounded numerically, as follows:

(4.36) \begin{align} |J(\tilde {x},\tilde {x}_b)| \leqslant \frac {1}{\pi } I((\tilde {x}-\tilde {x}_b)/\lambda )) \leqslant I(1) \approx 0.05. \end{align}

In (4.33), the magnitude of $J(\tilde {x},\tilde {x}_b)$ is small compared with $\cos (2\pi (\tilde {x}/\lambda +g) )$ , the other term in the bracket. Therefore, $J(\tilde {x},\tilde {x}_b)$ can be neglected and (4.33) can be approximated in a simple wave form (figure 14), as

(4.37) \begin{align} v(\tilde {x},0,t) \approx V_{w}(\tilde {x}) = \frac {2\pi ^2 \tilde {A}}{G(\sigma )} \cos \left (2\pi \left (\frac {\tilde {x}}{\lambda }+g(\sigma )\right )\right ). \end{align}

Moreover, for large $x$ , the asymptotic expansion $I(x)={1}/{(2\pi ^2 x)} + O(x^{-3})$ is obtained from expansions of $\mathrm{Si}(x)$ and $\mathrm{Ci}(x)$ . It then follows that $|J(\tilde {x},\tilde {x}_b)|\rightarrow 0$ for $(\tilde {x}-\tilde {x}_b)/\lambda \rightarrow +\infty$ , which indicates that the approximation $V_{w}(\tilde {x})$ in (4.37) is closer to $v(\tilde {x},0,t)$ for $\tilde {x}$ further away the wing trailing edge $\tilde {x}_b$ . Note that the wave $V_{w}(\tilde {x})$ (4.37) is a quarter lag in phase from the vorticity distribution $\gamma (\tilde {x})$ (4.29). Therefore, the phase difference between the wave $V_{w}(\tilde {x})$ and the trailing-edge trajectory $\tilde {y}_b(\tilde {x})$ (4.28) is then

(4.38) \begin{align} \delta _2(\sigma ) =\delta _1(\sigma )+1/4 = -g(\sigma )-\frac {\sigma }{2\pi }, \end{align}

and it has the limit $\delta _2(\sigma )\rightarrow 1/4$ as $\sigma \rightarrow 0$ .

The steady stationary wave $V_{w}(\tilde {x})$ in (4.37) provides a simple model of the flow at the centreline of a reverse von Kármán vortex street, in the limit of small Strouhal number $St=Af/U_{\infty }.$ For the parameters used in both the schooling experiments and our simulations, the flapping velocity $V_f=2\pi Af$ is small compared with the swimming velocity $U_{\infty }$ so that the Strouhal number is small. In the experiments of Ramananarivo et al. (Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016), ${\textit {St}}=0.1$ to $0.2$ , and in our simulations ${\textit {St}}=0.005$ to $0.1$ . From experimental flow visualisation of a single flapping wing (Becker et al. Reference Becker, Masoud, Newbolt, Shelley and Ristroph2015; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016) and our simulations (see the wake of the leader wing in figure 11), we find that the backward jet flow at the centreline in the reverse von Kármán street is weak, as the Strouhal number is small. Flow visualisation shows alternating upward and downward flow jets along the centreline of the wake generated by a single flapping wing. This allows us to use a waveform (4.37) to model the flow at the centreline of the wake produced by the leader wing in the school, and subsequently study its interaction with the follower wing.

4.1.3. Single flapping wing interacting with steady wave

Now we consider a single flapping wing interacting with a wavy stream in the background (figure 15). It serves as a model for the follower swimmer in the ‘school-of-two’, that is surrounded by the wake produced by the leader wing. In the frame of the wing, where the wing spans in the range $-1 \leqslant x \leqslant 1$ , we assume the wing is immersed in a wavy stream with velocity $(U_{\infty },V_{sw}(x,t)).$ The $x$ -component of this background flow is a uniform stream $U_{\infty },$ and the $y$ -component $V_{sw}$ is a travelling wave of the same velocity $U_{\infty }$ and the frequency is $f$ , the same frequency as the flapping wing. We denote the wavelength by $\lambda =U_{\infty }/f$ . Using dimensionless variables, where $f=1$ , the travelling wave is given by

(4.39) \begin{align} V_{sw}(x,t)=V_{0}\sin (2\pi (x/\lambda -t)),\quad \lambda =U_{\infty }/f, \end{align}

where $V_0$ is the amplitude of the wave, and the heaving motion of the wing has the form

(4.40) \begin{align} h(x,t)=\tilde {A}\cos (2\pi (t+\theta )),\quad -1\leqslant x \leqslant 1, \end{align}

where $\theta$ is a phase shift in time, and the flapping frequency ( $f=1$ ) is the same as the travelling wave.

The linearised problem of a rigid wing flapping in a wavy stream is solved by Wu & Chwang (Reference Wu and Chwang1975). There, the wing maintains a superposition motion of pitching and heaving, with the flapping frequency required to be the same as the frequency of the wavy stream. In this work, we adapt the calculation of Wu & Chwang (Reference Wu and Chwang1975) to our case of a purely heaving wing (4.40), flapping in a wavy stream $(U_{\infty },V_{sw})$ with the wave speed the same as the stream speed.

Figure 15.

A single flapping wing swimming with velocity $-U_{\infty }$ (negative $x$ -direction) interacts with a stationary wave $V_{sw}$ in the background. When the leading-edge velocity of the wing is in phase with the wave, the effective wing flapping is reduced and the wing thrust reaches a minimum value $T_{min }$ . When the wing leading edge is out of phase with the wave, the effective flapping speed is increased and the thrust reaches a maximum value $T_{max }$ .

We assume the heaving amplitude $\tilde {A}$ and the amplitude of the wave $V_0$ are small, and denote $\epsilon =V_{0}/(2\pi \tilde {A})=O(1)$ (using dimensionless variables; $\epsilon =V_{0}/(2\pi Af)$ if using dimensional variables) the ratio of the wave velocity to the flapping velocity. Similar to the linearisation of a single wing in a quiescent background (see § 4.1.1), we denote the flow velocity by $\boldsymbol{q} = (U_{\infty }+u,V_{sw}+v).$ Assuming $U_{\infty } \gg u$ and $U_{\infty } \gg v$ , the linearised 2-D Euler (4.2) has the following kinematic boundary condition:

(4.41) \begin{align} v_{1}(x,y,t)\mid _{y=\pm 0}=\partial _t h(x,t)-V_{sw}(x,t),\quad -1\leqslant x\leqslant 1. \end{align}

Using the same assumptions (4.6)–(4.8) in § 4.1.1, the acceleration potential $w(z,t)$ (defined in (4.4)) is solved by Wu & Chwang (Reference Wu and Chwang1975) through Fourier series analysis and conformal mapping techniques, a process similar to that in Wu (Reference Wu1961). The stroke-averaged thrust on the wing is provided by Wu & Chwang (Reference Wu and Chwang1975) as

(4.42) \begin{align} \langle T \rangle = -4\pi ^3 \tilde {A}^{2}\left |C(\sigma )\right |^{2}\left [ 1+\epsilon ^{2}P^{2}(\sigma )+2\epsilon P(\sigma )\sin (2\pi \left (\theta -p(\sigma )\right ))\right ], \end{align}

(4.43) \begin{align} \mathcal{P}(\sigma ) =\frac {\overline {K_{1}(i\sigma )}K_{0}(i\sigma )+\overline {K_{0}(i\sigma )}K_{1}(i\sigma )}{\pi K_{1}(i\sigma )}= P(\sigma )e^{2\pi ip(\sigma )}, \end{align}

where $\epsilon =V_{0}/(2\pi \tilde {A})$ , $\sigma =2\pi /U_{\infty },$ $P(\sigma )$ and $p(\sigma )$ are the modulus and angle of $\mathcal{P}(\sigma )$ , respectively. Note that the thrust (4.42) has a simple sinusoidal form of three variables $\theta$ , $\epsilon$ and $\sigma$ , where $\theta$ denotes the phase difference between the wave and the wing flapping, $\epsilon$ denotes the ratio of the flapping speed and the wave speed and $\sigma =2\pi /U_{\infty }=2\pi /\lambda$ ( $\sigma =\pi cf/U_{\infty }=\pi c/\lambda$ using dimensional quantities) represents the wavelength $\lambda$ of the wave, or stream speed $U_{\infty }$ .

In the flow frame, where the $x$ -component of the background flow is zero and the wing translates with velocity $-U_{\infty }$ in the negative $x$ -direction, by substituting (4.26) in (4.39), the background wave (4.39) now becomes a steady vertical wave

(4.44) \begin{align} \tilde {V}_{sw}(\tilde {x})=V_{0}\sin (2\pi \tilde {x}/\lambda ). \end{align}

The leading edge of the wing, parameterised in $\tilde {t}$ as

(4.45) \begin{align} \tilde {x} = -U_{\infty } \tilde {t}-1,\quad \tilde {y} = \tilde {A}\cos (2\pi (\tilde {t}+\theta )), \end{align}

follows the trajectory

(4.46) \begin{align} \tilde {y}(\tilde {x})=\tilde {A}\cos (2\pi \left ((\tilde {x}+1)/\lambda -\theta \right )). \end{align}

Using (4.46) and (4.44), we obtain the phase lag of the wing’s leading edge $\tilde {y}(\tilde {x})$ to the vertical wave $\tilde {V}_{sw}(\tilde {x})$ , that is,

(4.47) \begin{align} \tilde {\theta } = \theta - \sigma /(2\pi ) - 1/4, \end{align}

where $\sigma = 2\pi /\lambda ,$ and the thrust on the wing (4.42) can be expressed in terms of $\tilde {\theta }$ as

(4.48) \begin{align} \langle T \rangle = -4\pi ^3 \tilde {A}^{2}\left |C(\sigma )\right |^{2}\left [ 1+\epsilon ^{2}P^{2}(\sigma )+2\epsilon P(\sigma )\cos \left (2\pi \left (\tilde {\theta }-\delta _3(\sigma )\right )\right )\right ], \end{align}

(4.49) \begin{align} \delta _3(\sigma ) = p(\sigma )-\sigma /(2\pi ). \end{align}

It is shown by (4.48) that the extrema of $\langle T \rangle$ are determined by the phase difference between the wing and the wave $\tilde {\theta }$ , which is of particular interest. Extreme values of thrust $\langle T \rangle$ are

(4.50) \begin{align} \langle T \rangle _{max } & = -4\pi ^3 \tilde {A}^{2}\left |C(\sigma )\right |^{2}\left [ 1+\epsilon P(\sigma )\right ]^2, \quad \text{when} \quad \tilde {\theta }=\delta _3 + k, \end{align}

(4.51) \begin{align} \langle T \rangle _{min } & = -4\pi ^3 \tilde {A}^{2}\left |C(\sigma )\right |^{2}\left [ 1-\epsilon P(\sigma )\right ]^2, \quad \text{when} \quad \tilde {\theta }=\delta _3 + k + 1/2, \end{align}

where $k\in \mathbb{Z}^+$ . For $\lambda \in [2,\infty )$ , i.e. $\sigma = 2\pi /\lambda \in (0,\pi ]$ , which is the main range of interest, the numerical value is approximately $0.143 \lesssim \delta _3 \leqslant 1/4$ , which is approximately a quarter. Moreover, in the limit of $\sigma \rightarrow 0$ , we obtain $P(\sigma )\rightarrow 1$ and $p(\sigma )\rightarrow 1/4$ by substituting asymptotic expansions of the modified Bessel functions (4.16)–(4.17) in (4.43), and thus $\delta _3(\sigma ) \rightarrow 1/4$ . These values indicate that the maximum thrust is achieved when the wing leading-edge trajectory (4.46) has phase-lag $\tilde {\theta } \approx k+1/4$ to the stationary wave (4.44), and the minimum thrust is achieved when the leading-edge trajectory has about $\tilde {\theta }\approx k+3/4$ phase lag to the wave.

This could be explained by looking at the flapping velocity of the wing leading edge, which can be parameterised as

(4.52) \begin{align} \tilde {x} = -U_{\infty } \tilde {t}-1,\quad \dot {\tilde {y}} = -2\pi \tilde {A}\sin (2\pi (\tilde {t}+\theta )), \end{align}

by (4.26), and be expressed in $\tilde {x}$ as

(4.53) \begin{align} \dot {\tilde {y}}(\tilde {x})=2\pi \tilde {A}\sin \big (2\pi \big (\tilde {x}/\lambda -\tilde {\theta }-1/4 \big )\big ). \end{align}

We define an ‘effective flapping’ velocity at the wing leading edge

(4.54) \begin{align} V_{\textit{eff}}(\tilde {x})= \dot {\tilde {y}}(\tilde {x})-\tilde {V}_{sw}(\tilde {x})=2\pi \tilde {A}\sin \big (2\pi \big (\tilde {x}/\lambda -\tilde {\theta }-1/4 \big )\big )-V_{0}\sin (2\pi \tilde {x}/\lambda ), \end{align}

which is the relative velocity of the wing flapping (4.53) and the background wave (4.44). Note that the maximum thrust (4.50) is achieved when $\tilde {\theta } \approx k+1/4$ , at which time the flapping velocity of the leading edge (4.53) is out of phase with the wave, and the effective flapping has a maximum amplitude

(4.55) \begin{align} V_{{\textit{eff}},max }(\tilde {x}) = (2\pi \tilde {A}+ V_0)\sin \left (2\pi \tilde {x}/\lambda \right ), \end{align}

as shown in the schematic diagram in figure 15. Similarly, the minimum thrust corresponds to the heaving velocity at the leading edge in phase with the wave, and the amplitude of effective flapping achieves minimum

(4.56) \begin{align} V_{{\textit{eff}},min }(\tilde {x})=(2\pi \tilde {A}-V_0)\sin \left (2\pi \tilde {x}/\lambda \right ). \end{align}

These results are not surprising. As the thrust $\langle T \rangle$ consists only of the suction force at the leading edge of the wing, when the background wave is out of phase with the flapping velocity of the wing leading edge, the effective flapping is increased and the wing thrust is thus increased and reaches a maximum. Similarly when the wave is in phase with the flapping of the wing leading edge, the effective flapping is decreased and the thrust reaches a minimum.

4.2.1. Hydrodynamic thrust of the follower wing

We denote the distance between the follower’s head and the leader’s tail by $g$ (as in figure 16). The schooling number defined as $S=g/\lambda$ ( $\lambda =U_{\infty }/f$ ) denotes the phase lag of the follower’s head to leader’s tail, following which the phase lag of the follower’s leading edge to the stationary wave generated by the leader is

(4.57) \begin{align} \tilde {\theta } = S - \delta _2(\sigma ). \end{align}

Substituting (4.57) in (4.48), and by (4.24), (4.38), (4.43) and (4.49), we express the thrust on the follower in schooling number $S$ as

(4.58) \begin{align} \langle T_2 \rangle = \langle T_1 \rangle \left [ \varepsilon ^{2}+Q^{2}(\sigma )+2\varepsilon Q(\sigma )\cos \left (2\pi \left (S - \delta _4(\sigma )\right )\right )\right ],\quad \varepsilon =\tilde {A}_2/\tilde {A}_1, \end{align}

(4.59) \begin{align} \delta _4(\sigma ) = \delta _2(\sigma )+\delta _3(\sigma )=q(\sigma )-\sigma /\pi , \end{align}

where $\langle T_1 \rangle = -4\pi ^3 \tilde {A}_1^{2}\left |C(\sigma )\right |^{2}$ is the thrust of the leader wing (4.9), $Q(\sigma )$ and $q(\sigma )$ are the modulus and angle of

(4.60) \begin{align} \mathcal{Q}(\sigma ) = \frac {\pi \mathcal{P}(\sigma )}{\mathcal{G}(\sigma )} =\frac {\overline {K_{1}(i\sigma )}K_{0}(i\sigma )+\overline {K_{0}(i\sigma )}K_{1}(i\sigma )}{K_{1}(i\sigma )\left [K_{0}(i\sigma )+K_{1}(i\sigma )\right ]} =Q(\sigma )e^{2\pi iq(\sigma )}, \end{align}

respectively. In figure 17(a), the follower’s thrust (black solid curve) calculated from (4.58) is compared with the thrust computed by the vortex-sheet simulations (grey solid curve) in § 3, in which the thrust is measured on a ‘ghost follower’ maintaining a fixed separation to the leader. Here, the two wings have the same amplitude $\tilde {A}=0.2$ , the drag coefficient is fixed at $C_f=0.02$ and $\sigma$ in (4.58) is calculated from (4.19).

Figure 17.

Comparisons between the hydrodynamic forces and potential on the follower wing calculated from linearised theory (black) and from simulations (grey). (a) Hydrodynamic thrust $\langle T_2\rangle$ (solid lines) and drag $-\langle D_2\rangle$ (dashed lines) as functions of schooling number $S$ , where the two wings have the same flapping amplitude $\tilde{A}=0.2$ . The forces here are dimensionless. (b) Total hydrodynamic force on the follower $\langle F_2\rangle$ normalised by the single-wing thrust $|\langle T_0\rangle |$ . (c) Hydrodynamic potential on the follower. In (b)–(c), the leader’s amplitude is fixed at $\tilde {A}_1=0.2$ , while the follower’s amplitude $\tilde {A}_2$ is varied. Here, $\varepsilon =\tilde {A}_2/\tilde {A}_1=1-Q$ (red), $1$ (black and grey), and $1+Q$ (blue). In all cases the drag coefficient is $C_f=0.02$ .

Note that the follower’s thrust (4.58) is a sinusoidal function in schooling number $S$ , i.e. it depends on the phase difference between the two wings swimming. The coefficients of this sinusoidal function depend on the amplitude ratio of the wings $\varepsilon =\tilde {A}_2/\tilde {A}_1$ and the wavelength of the motion $\lambda /c=2 \pi /\sigma$ (dimensionless). For $\varepsilon =0$ , (4.58) simplifies to

(4.61) \begin{align} \langle T_2 \rangle |_{\varepsilon =0} = \langle T_1 \rangle Q^2(\sigma )=-4\pi ^3 \left |C(\sigma )\right |^{2}\left ( Q(\sigma )\tilde {A}_1\right )^2. \end{align}

This implies that when the follower wing is a non-flapping ‘dead wing’, if it somehow still maintains a swimming speed $U_{\infty }$ , it can generate a non-zero constant thrust due to the interaction with the wavy wake shed by the leader. This ‘dead-follower’ thrust $\langle T_2 \rangle |_{\varepsilon =0}\to 0$ at the long wavelength limit $\sigma \rightarrow 0$ , as $Q(\sigma )\rightarrow 0$ ((4.16)–(4.17) and (4.60)).

The follower’s thrust (4.58) achieves extreme values

(4.62) \begin{align} \langle T_2 \rangle _{max } & = \langle T_1 \rangle \left [ \varepsilon +Q(\sigma )\right ]^2=-4\pi ^3 \left |C(\sigma )\right |^{2}\left [ \tilde {A}_2+Q(\sigma )\tilde {A}_1\right ]^2,\; \text{when} \; S=\delta _4+k,\nonumber \\ \langle T_2 \rangle _{min } & = \langle T_1 \rangle \left [ \varepsilon -Q(\sigma )\right ]^2=-4\pi ^3 \left |C(\sigma )\right |^{2}\left [ \tilde {A}_2-Q(\sigma )\tilde {A}_1\right ]^2,\; \text{when} \; S=\delta _4+k+1/2, \end{align}

where $k\in \mathbb{Z}^+$ . In the limiting case of $\sigma \rightarrow 0$ , we have $q(\sigma )\rightarrow 1/2$ and thus $\delta _4\rightarrow 1/2$ (4.59). This implies that at the long wavelength limit, $\lambda /c\rightarrow \infty$ , the maximum thrust $\langle T_2 \rangle _{max }$ is achieved when the schooling number $S\approx 1/2+k$ and minimum thrust $\langle T_2 \rangle _{max }$ occurs at $S\approx k$ . Moreover, we note that the perturbation due to the schooling is weak as $Q(\sigma )\rightarrow 0$ when $\sigma \rightarrow 0$ . When $\sigma =0$ , the follower’s thrust $\langle T_2 \rangle = -4\pi ^3 \tilde {A}_2^2|C(\sigma )|^{2}$ is the thrust of a single wing maintaining the follower’s flapping amplitude $\tilde {A}_2$ .

4.2.2. Hydrodynamic potential, schooling numbers and ‘keep-up’ schooling condition

We approximate the free-stream velocity on wing surfaces, $U_{\infty }+u$ in (4.2), by the wing swimming speed $U_{\infty }$ , ignoring the small perturbation in horizontal velocity $u$ , as in the linearised calculation. The skin-friction $D_k(|U_{\infty }|),\ k=1,2$ is then a function of $U_{\infty }$ only. As the two wings swim at the same speed, they experience the same amount of drag $D_1=D_2$ . Note that the stroke-averaged thrust and drag are always balanced on the leader wing, $\langle T_{1} \rangle +D_1=0$ (4.19). It then follows that $0=\langle T_{1} \rangle +D_1=\langle T_{1} \rangle +D_2$ . Using this relation and (4.58), we obtain the following expression for the total hydrodynamic force on the follower:

(4.63) \begin{align} \langle F_2 \rangle = \langle T_2 \rangle + D_2 =\langle T_1 \rangle \left [ \varepsilon ^{2}+Q^{2}(\sigma )-1+2\varepsilon Q(\sigma )\cos \left (2\pi \left (S - \delta _4(\sigma )\right )\right )\right ]. \end{align}

In figure 17(a), the approximation of the follower’s drag $D_2=-\langle T_{1} \rangle$ , calculated by (4.9) (black dashed curve), is compared with the drag computed from the simulations (grey dashed curve), which shows the theoretical approximation is close to the simulation when the schooling number $S$ is not too small, i.e. the follower is not too close to the leader. It then follows that the total hydrodynamic force formula (4.63) serves as a good approximation to the simulations, as shown in figure 17(b).

The hydrodynamic potential $\phi (S)$ is the negative integral of the the hydrodynamic force $\langle F_2 \rangle$ , given as

(4.64) \begin{align} \phi (S) = -\int \langle F_2 \rangle \,{\mathrm{d}} S=-\langle T_1 \rangle \left [\left ( \varepsilon ^{2}+Q^{2}(\sigma )-1\right )S+\frac {\varepsilon Q(\sigma )}{\pi }\sin \left (2\pi \left (S - \delta _4(\sigma )\right )\right )\right ]. \end{align}

In figures 17(b) and 17(c), the hydrodynamic force calculated from (4.63) and hydrodynamic potential from (4.64) are shown, and compared with simulations of two wings having the same flapping amplitude (the same simulation data as in figure 17 b). Here, $\sigma$ in (4.63)–(4.64) is solved from (4.19) given a flapping amplitude $\tilde {A}=0.2$ and a drag coefficient $C_f=0.02$ . For other values of $C_f$ (ranging from 0.01 to 0.1), the system also exhibits multiple stable equilibria, producing a corrugated potential landscape $\phi (S)$ characterised by an array of stable wells, similar to that shown in figure 17(c) (for $\varepsilon =1$ ). The positions of the local minima shift only moderately with $C_f$ , moving toward slightly smaller $S$ as $C_f$ increases. Moreover, the widths of these wells, which corresponds to their basins of attraction, show minimal variation with $C_f$ .

The steady schooling numbers are roots of $\langle F_2 \rangle$ , or local extrema of $\phi (S)$ , which can be solved from (4.63)

(4.65) \begin{align} S_{k}^{\pm }=\pm \frac {1}{2\pi }\arccos \left (\frac {1-\varepsilon ^2-Q^2(\sigma )}{2\varepsilon Q(\sigma )}\right )+q(\sigma )-\frac {\sigma }{\pi } + k,\quad k\in \mathbb{Z}^+. \end{align}

Moreover, it can be shown from (4.64) that $S_{k}^{-}$ are stable schooling numbers (stable critical points) and $S_{k}^{+}$ are unstable (see figure 18). If the two wings have the same flapping amplitude, $\varepsilon = \tilde {A}_2/\tilde {A}_1=1$ , (4.65) simplifies to

(4.66) \begin{align} S_{k}^{\pm }=\pm \frac {1}{2\pi }\arccos \left (-Q(\sigma )/2\right )+q(\sigma )-\frac {\sigma }{\pi }+ k,\quad k\in \mathbb{Z}^+. \end{align}

This is compared in figure 18(a,b) with simulations of two wings having the same flapping amplitude, where $\sigma$ is solved from (4.19) given a flapping amplitude $\tilde {A}$ and a drag coefficient $C_f$ . At the long wavelength limit $\sigma =\pi c/\lambda \rightarrow 0$ , we have $Q(\sigma )\rightarrow 0$ and $q(\sigma )\rightarrow 1/2$ , so that the schooling numbers achieve limits

(4.67) \begin{align} S_{k}^{-}\rightarrow k+\frac {1}{4},\quad S_{k}^{+}\rightarrow k+\frac {3}{4},\quad k\in \mathbb{Z}^+. \end{align}

Figure 18.

Schooling numbers computed from simulations (symbols) are compared with the linearised theory calculations (lines). Stable schooling numbers $S_k^-$ are black and unstable schooling numbers $S_k^+$ are grey. (a)–(b) Flapping amplitudes are the same on both wings, the amplitude $\tilde {A}$ and drag coefficient $C_f$ are varied. In (b), integers $k=1,2,3$ are subtracted from stable schooling numbers $S_k^-$ . The limit $S_k^{-}-k\rightarrow 1/4$ as $\lambda /c\rightarrow \infty$ . The simulation data are the same as in figure 7. (c) Fix leader’s amplitude $\tilde {A}_1=0.2$ , and vary follower’s amplitude $\tilde {A}_2$ . Stable and unstable schooling numbers exist when $1-Q\leqslant \tilde {A}_2/\tilde {A}_1 \leqslant 1+Q$ . The simulation data are the same as in figure 12(a). The drag coefficient is $C_f=0.02$ .

Note that if the two wings have different flapping amplitudes, (4.65) indicates a condition for steady schooling states, i.e.

(4.68) \begin{align} -1\leqslant \frac {1-\varepsilon ^2-Q^2(\sigma )}{2\varepsilon Q(\sigma )} \leqslant 1. \end{align}

This provides a ‘keep-up’ condition for the ratio of the follower’s amplitude $\tilde {A}_2$ to the leader’s amplitude $\tilde {A}_1$

(4.69) \begin{align} 1-Q(\sigma )\leqslant \tilde {A}_2/\tilde {A}_1 \leqslant 1+Q(\sigma ). \end{align}

As figure 18(c) shows, the schooling number exists only for $1-Q(\sigma )\leqslant \varepsilon \leqslant 1+Q(\sigma )$ . Stable and unstable branches merge at $\varepsilon =1\pm Q(\sigma ).$ It can also be seen in figure 17(b,c) that for $\varepsilon =1\pm Q(\sigma ),$ the hydrodynamic force $\langle F_2 \rangle$ in (4.63) has only one root in each period, where the stable and unstable solutions collapse. Beyond the range of (4.69), $\langle F_2 \rangle$ does not intersect with the $x$ -axis. When $\tilde {A}_2/\tilde {A}_1 \lt 1-Q(\sigma )$ , the hydrodynamic force on the follower wing is always strictly positive (always in the opposite direction of swimming). In this case, the follower cannot generate enough force maintaining the same swimming speed $U_{\infty }$ as the leader, so that the follower is not able to ‘keep up’ in the school and will drift away downstream. It can also be explained by the hydrodynamic potential profile $\phi (S)$ in figure 17(c), that $\phi (S)$ is a monotonically decreasing function in $S$ , and wherever the follower locates in this potential the schooling number $S$ could not be trapped and will keep increasing. For $\tilde {A}_2/\tilde {A}_1 \gt 1+Q(\sigma )$ , the hydrodynamic force is always strictly negative (always in the swimming direction), therefore the follower will keep accelerating, eventually leading to ‘rear-ending’ collisions within the school. Tests with other values of skin-friction drag coefficient $C_f$ show that the ‘separation’ and ‘collision’ regions (as in figure 18 c) shrink as $C_f$ increases. This indicates that the follower is more likely to maintain a stable formation with the leader at higher $C_f$ , and that the schooling number exists over a wider range of $\varepsilon$ , i.e. the stability of the schooling configuration improves with increasing $C_f$ .