2 Formulation of the problem. Kinematics

We consider a two-dimensional, incompressible and nearly inviscid flow over a flexible thin foil of chord length $c$ that moves with constant speed $U$ along the negative $x$-axis. Superimposed to heaving and pitching motions (see figure 1a), the plate undergoes an undulating motion with a given wavelength (defined below), and a quadratic flexure (or deflection) motion of trailing-edge amplitude $\unicode[STIX]{x1D6FF}_{m}$ (see figure 1b). All the amplitudes are small compared with the chord $c$, so that the airfoil, and every point of the trail of vortices that it leaves behind, may be considered to be on the horizontal plane $z=0$ to a first approximation.

Figure 1. Schematic of the oscillating foil for heaving and pitching rigid motion (a) and for heaving and flexural motion (b), both in the absence of undulatory motion (i.e. for $\mathsf{b}=0$). See (2.1)–(2.4) for the meaning of the various non-dimensional quantities.

We use non-dimensional variables scaled with the half-chord length $c/2$ and the velocity $U$, so that the plate extends from $x=-1$ to $x=1$ in a reference frame translating with it at speed $U$ along the $x$-axis (see figure 1), and the non-dimensional time $t$ is scaled with $c/(2U)$. In this reference frame, the vertical motion of the foil, given by the vertical displacement of its mean-chamber line $z_{s}$, can be written as

(2.1)$$\begin{eqnarray}z_{s}(x,t)=h(x,t)-(x-a)\unicode[STIX]{x1D6FC}(x,t)+(x-p)^{2}\unicode[STIX]{x1D6FF}(x,t),\quad -1\leqslant x\leqslant 1,\end{eqnarray}$$

with

(2.2a-c)$$\begin{eqnarray}h(x,t)=\text{Re}[h_{0}\text{e}^{\text{i}kt+\mathsf{b}(x+1)}],\quad \unicode[STIX]{x1D6FC}(x,t)=\text{Re}[\unicode[STIX]{x1D6FC}_{0}\text{e}^{\text{i}kt+\mathsf{b}(x+1)}],\quad \unicode[STIX]{x1D6FF}(x,t)=\text{Re}[\unicode[STIX]{x1D6FF}_{0}\text{e}^{\text{i}kt+\mathsf{b}(x+1)}],\end{eqnarray}$$

(2.3a,b)$$\begin{eqnarray}\mathsf{b}=b_{1}-\text{i}b_{2},\quad k=\frac{\unicode[STIX]{x1D714}c}{2U}.\end{eqnarray}$$

In these expressions $\text{Re}$ means real part and $k$ is the reduced frequency associated with the frequency $\unicode[STIX]{x1D714}$ of the harmonic motion, with non-dimensional period $T=2\unicode[STIX]{x03C0}/k$, which consists of a heaving displacement $h(x,t)$, a pitching rotation $\unicode[STIX]{x1D6FC}(x,t)$ pivoting at $x=a$ and a quadratic deflection motion $\unicode[STIX]{x1D6FF}(x,t)$ pivoting at $x=p$. Superimposed to each of these three displacements, we have included an undulatory motion of the foil with non-dimensional wavenumber $b_{2}$ whose amplitude grows exponentially to the trailing edge with a non-dimensional factor $b_{1}$, both quantities grouped together in the complex number $\mathsf{b}$. This form, motivated by the analysis of the carangiform fish dynamics by Lighthill (Reference Lighthill1970), produces a travelling wave of growing amplitude that propagates from the leading edge to the trailing edge with non-dimensional phase speed $k/b_{2}$. Figure 2 shows some examples when the pitching and flexure motions are absent.

Figure 2. Non-dimensional vertical displacement $z_{s}(x,t)$ for an undulatory heaving motion $(\unicode[STIX]{x1D6FC}_{0}=\unicode[STIX]{x1D6FF}_{0}=0)$ at different instants of the period’s oscillation $T$ for $h_{0}=0.05$ and $b_{1}=0.8$, with $b_{2}=1$ (a) and with $b_{2}=\unicode[STIX]{x03C0}$ (b) (see (2.1)–(2.4) for the non-dimensional notation).

The amplitudes $h_{0},\unicode[STIX]{x1D6FC}_{0}$ and $\unicode[STIX]{x1D6FF}_{0}$ are, in general, complex constants satisfying $|h_{0}|\ll 1$, $|\unicode[STIX]{x1D6FC}_{0}|\ll 1$ and $|\unicode[STIX]{x1D6FF}_{0}|\ll 1$. For simplicity we select $h_{0}$ real and

(2.4a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{0}=a_{0}\text{e}^{\text{i}\unicode[STIX]{x1D719}},\quad \unicode[STIX]{x1D6FF}_{0}=\frac{\unicode[STIX]{x1D6FF}_{m}}{(1-p)^{2}}\text{e}^{\text{i}\unicode[STIX]{x1D713}},\end{eqnarray}$$

with $\unicode[STIX]{x1D719}$ the phase shift between the heaving and pitching motions of the foil, $\unicode[STIX]{x1D713}$ the phase shift between the heaving and deflection motions, $a_{0}$ the maximum pitch amplitude of the plate and $\unicode[STIX]{x1D6FF}_{m}$ the maximum amplitude of the flexure component of the motion at the trailing edge ($x=1$). In what follows we shall work with the complex expressions knowing that we have to take the real part of the results. Equations (2.1)–(2.4) describe a fairly broad class of the flapping undulatory motion of a flexible or compliant foil, with nine non-dimensional kinematic parameters (plus the reduced frequency), for which we shall derive analytical expressions for the force and moment using the impulse theory in the linear potential limit. General undulatory kinematics, with infinitely many wavenumbers, have been considered numerically by several authors both in the present small-amplitude inviscid limit (e.g. Alben Reference Alben2008; Moore Reference Moore2017; Tzezana & Breuer Reference Tzezana and Breuer2019) and using direct numerical simulation of the Navier–Stokes equations (e.g. Hoover et al. Reference Hoover, Cortez, Tytell and Fauci2018).

To facilitate the computations, the vertical displacement of the foil will be written as

(2.5)$$\begin{eqnarray}\displaystyle z_{s}(x,t)=[{\mathcal{F}}(t)+{\mathcal{E}}(t)x+{\mathcal{D}}(t)x^{2}]\text{e}^{\mathsf{b}x},\quad -1\leqslant x\leqslant 1, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{F}}(t),{\mathcal{E}}(t)$ and ${\mathcal{D}}(t)$ are given by

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{F}}(t)\equiv (h_{0}+a\unicode[STIX]{x1D6FC}_{0}+p^{2}\unicode[STIX]{x1D6FF}_{0})\text{e}^{\text{i}kt+\mathsf{b}}\equiv {\mathcal{F}}_{0}\text{e}^{\text{i}kt+\mathsf{b}}, & \displaystyle\end{eqnarray}$$

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}(t)\equiv -(\unicode[STIX]{x1D6FC}_{0}+2p\unicode[STIX]{x1D6FF}_{0})\text{e}^{\text{i}kt+\mathsf{b}}\equiv {\mathcal{E}}_{0}\text{e}^{\text{i}kt+\mathsf{b}}, & \displaystyle\end{eqnarray}$$

(2.8)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{D}}(t)\equiv \unicode[STIX]{x1D6FF}_{0}\text{e}^{\text{i}kt+\mathsf{b}}. & \displaystyle\end{eqnarray}$$

The corresponding non-dimensional vertical velocity of the foil’s mean-chamber line can be written as

(2.9)$$\begin{eqnarray}\displaystyle & \displaystyle v_{0}(x,t)=\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}x}=[{\mathcal{A}}(t)+{\mathcal{B}}(t)x+{\mathcal{C}}(t)x^{2}]\text{e}^{\mathsf{b}x},\quad -1\leqslant x\leqslant 1, & \displaystyle\end{eqnarray}$$

where ${\mathcal{A}}(t),{\mathcal{B}}(t)$ and ${\mathcal{C}}(t)$ are given by

(2.10)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}(t)\equiv \{(\text{i}k+\mathsf{b}){\mathcal{F}}_{0}+{\mathcal{E}}_{0}\}\text{e}^{\text{i}kt+\mathsf{b}}\equiv {\mathcal{A}}_{0}\text{e}^{\text{i}kt+\mathsf{b}}, & \displaystyle\end{eqnarray}$$

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}(t)\equiv \{(\text{i}k+\mathsf{b}){\mathcal{E}}_{0}+2\unicode[STIX]{x1D6FF}_{0}\}\text{e}^{\text{i}kt+\mathsf{b}}\equiv {\mathcal{B}}_{0}\text{e}^{\text{i}kt+\mathsf{b}}, & \displaystyle\end{eqnarray}$$

(2.12)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{C}}(t)\equiv (\text{i}k+\mathsf{b})\unicode[STIX]{x1D6FF}_{0}\text{e}^{\text{i}kt+\mathsf{b}}\equiv {\mathcal{C}}_{0}\text{e}^{\text{i}kt+\mathsf{b}}. & \displaystyle\end{eqnarray}$$

3 General expressions for the lift, thrust, moment and input power

The vortical impulse theory for an incompressible and unbounded flow (Wu Reference Wu1981) is used to obtain the force and moment on the airfoil. Neglecting the volume (section) of the airfoil, the total force $\boldsymbol{F}$ can be written as

(3.1)$$\begin{eqnarray}\boldsymbol{F}\equiv D\boldsymbol{e}_{x}+L\boldsymbol{e}_{z}=-\unicode[STIX]{x1D70C}\frac{\text{d}\boldsymbol{I}}{\text{d}t},\end{eqnarray}$$

where the $x$-component $D$ is the drag (or minus the thrust), the $z$-component $L$ is the lift, $\unicode[STIX]{x1D70C}$ is the fluid density and the vortical impulse (or vorticity moment) $\boldsymbol{I}$ is defined as

(3.2)$$\begin{eqnarray}\boldsymbol{I}=\int _{{\mathcal{V}}}\boldsymbol{x}\wedge \unicode[STIX]{x1D74E}\,\text{d}{\mathcal{V}},\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\wedge \boldsymbol{v}$ is the vorticity field and ${\mathcal{V}}$ is the entire volume (plane surface in this case) occupied by the fluid plus the airfoil. In writing (3.1) it is assumed that ${\mathcal{V}}$ is unbounded and that the flow is irrotational far from the airfoil. In fact, we shall assume that the vorticity, which is directed along the normal $\boldsymbol{e}_{y}$ to the plane of the fluid motion, is concentrated on the airfoil surface and in its trailing wake, both of them considered as planar (bound and free, respectively) vortex sheets. Thus,

(3.3)$$\begin{eqnarray}\displaystyle \boldsymbol{I}\simeq \int _{-1}^{1}(-z_{s}\unicode[STIX]{x1D71B}_{s}\boldsymbol{e}_{x}+x\unicode[STIX]{x1D71B}_{s}\boldsymbol{e}_{z})\,\text{d}x+\int _{1}^{\infty }(-z_{e}\unicode[STIX]{x1D71B}_{e}\boldsymbol{e}_{x}+x\unicode[STIX]{x1D71B}_{e}\boldsymbol{e}_{z})\,\text{d}x, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D71B}_{s}(x,t)$ is the vorticity density distribution on the airfoil, $\unicode[STIX]{x1D71B}_{e}(x,t)$ is the vorticity density distribution in the trailing wake and $z_{e}(x,t)$ is the vertical position of each point in this vortex wake (see figure 1). We consider the large-time behaviour in which the vortex wake sheet extends many chord lengths downstream of the airfoil, so that, to a first approximation, $1\leqslant x\leqslant \infty$ for both $\unicode[STIX]{x1D71B}_{e}(x,t)$ and $z_{e}(x,t)$, with $|z_{e}|\ll 1$ in the present linear approximation. Consequently, under the assumptions made, the total drag and lift force components can be written as

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle D=\unicode[STIX]{x1D70C}\frac{\text{d}}{\text{d}t}\left[\int _{-1}^{1}z_{s}\unicode[STIX]{x1D71B}_{s}\,\text{d}x+\int _{1}^{\infty }z_{e}\unicode[STIX]{x1D71B}_{e}\,\text{d}x\right], & \displaystyle\end{eqnarray}$$

(3.5)$$\begin{eqnarray}\displaystyle & \displaystyle L=-\unicode[STIX]{x1D70C}\frac{\text{d}}{\text{d}t}\left[\int _{-1}^{1}x\unicode[STIX]{x1D71B}_{s}\,\text{d}x+\int _{1}^{\infty }x\unicode[STIX]{x1D71B}_{e}\,\text{d}x\right]. & \displaystyle\end{eqnarray}$$

Similarly, the vortical impulse theory also provides the total moment on the airfoil (Wu Reference Wu1981):

(3.6)$$\begin{eqnarray}\boldsymbol{M}=-M\boldsymbol{e}_{y}=-\unicode[STIX]{x1D70C}\frac{\text{d}\boldsymbol{A}}{\text{d}t},\end{eqnarray}$$

where

(3.7)$$\begin{eqnarray}\boldsymbol{A}=-\frac{1}{2}\int _{{\mathcal{V}}}|\boldsymbol{x}-a\boldsymbol{e}_{x}|^{2}\,\text{d}{\mathcal{V}}\end{eqnarray}$$

is the angular impulse in relation to the pitching axis $x=a$ moving with speed $U$ along the $x$-axis (note that the distance $Ut$ is also scaled with $c/2$). Thus, on using the same approximations made in (3.3),

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle M\simeq \frac{1}{2}\unicode[STIX]{x1D70C}\frac{\text{d}}{\text{d}t}\left[\int _{-1}^{1}(x-a)^{2}\unicode[STIX]{x1D71B}_{s}\,\text{d}x+\int _{1}^{\infty }(x-a)^{2}\unicode[STIX]{x1D71B}_{e}\,\text{d}x\right]. & \displaystyle\end{eqnarray}$$

Finally, the input power is given by

(3.9)$$\begin{eqnarray}P=\int _{-1}^{1}\unicode[STIX]{x0394}p\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x,\end{eqnarray}$$

where $\unicode[STIX]{x0394}p(x,t)\equiv p_{+}(x,t)-p_{-}(x,t)$ is the local pressure difference between both sides of the foil, with subscripts $+$ and $-$ denoting its upper and lower surfaces, respectively. Pressure difference $\unicode[STIX]{x0394}p$ can be obtained from the unsteady Bernoulli equation on the foil surface as

(3.10)$$\begin{eqnarray}\unicode[STIX]{x0394}p=-\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7})-\frac{1}{2}\unicode[STIX]{x1D70C}(u_{+}^{2}-u_{-}^{2}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}$ is the velocity potential and $u$ the tangential velocity component. Taking into account that

(3.11a-c)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\int _{-1}^{x}(u_{+}-u_{-})\,\text{d}x,\quad \unicode[STIX]{x1D71B}_{s}=u_{+}-u_{-},\quad U={\textstyle \frac{1}{2}}(u_{+}+u_{-}),\end{eqnarray}$$

and substituting into (3.9), the input power can be written in terms of the vorticity distribution as

(3.12)$$\begin{eqnarray}P=-\unicode[STIX]{x1D70C}U\int _{-1}^{1}\unicode[STIX]{x1D71B}_{s}\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x-\unicode[STIX]{x1D70C}\int _{-1}^{1}\left(\int _{-1}^{x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D71B}_{s}}{\unicode[STIX]{x2202}t}\,\text{d}\unicode[STIX]{x1D709}\right)\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x,\end{eqnarray}$$

which, after integrating by parts the double integral, can be more conveniently written as

(3.13)$$\begin{eqnarray}P=-\unicode[STIX]{x1D70C}U\int _{-1}^{1}\unicode[STIX]{x1D71B}_{s}\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x-\unicode[STIX]{x1D70C}\int _{-1}^{1}\left(\int _{x}^{1}\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}\unicode[STIX]{x1D709}\right)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D71B}_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x.\end{eqnarray}$$

3.1 Vorticity distribution

Following von Kármán & Sears (Reference von Kármán and Sears1938), who used the above linearized version of the impulse theory in the inviscid limit for the lift and for the moment on a rigid airfoil (but not for the thrust) much earlier than Wu’s (Reference Wu1981) general impulse formulation was derived, and invoking the linearity of the problem, we separate the different contributions to the vorticity on the airfoil as

(3.14)$$\begin{eqnarray}\unicode[STIX]{x1D71B}_{s}(x,t)=\unicode[STIX]{x1D71B}_{0}(x,t)+\unicode[STIX]{x1D71B}_{se}(x,t),\end{eqnarray}$$

where the first term $\unicode[STIX]{x1D71B}_{0}$, with associated circulation

(3.15)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{0}(t)=\int _{-1}^{1}\unicode[STIX]{x1D71B}_{0}(x,t)\,\text{d}x,\end{eqnarray}$$

is the quasi-steady contribution, i.e. that for an airfoil without considering the effect of its unsteady wake, such that the corresponding lift would be $\unicode[STIX]{x1D70C}U\unicode[STIX]{x1D6E4}_{0}$. The last term in (3.14), $\unicode[STIX]{x1D71B}_{se}$, is the contribution to $\unicode[STIX]{x1D71B}_{s}$ induced by its wake vortex sheet, of strength $\unicode[STIX]{x1D71B}_{e}$ (remember from (3.3) that subscript $s$ stands for the surface of the foil and $e$ for the trailing wake).

Kelvin’s total-circulation conservation theorem requires that

(3.16)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{0}+\unicode[STIX]{x1D6E4}_{se}+\int _{1}^{\infty }\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)\,\text{d}\unicode[STIX]{x1D709}=0,\end{eqnarray}$$

with

(3.17)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{se}(t)=\int _{-1}^{1}\unicode[STIX]{x1D71B}_{se}(x,t)\,\text{d}x.\end{eqnarray}$$

To obtain $\unicode[STIX]{x1D71B}_{0}$ and $\unicode[STIX]{x1D71B}_{se}$ one has to apply the boundary condition of the vertical velocity (2.9) at $z=0$, for $-1\leqslant x\leqslant 1$, which is induced by the whole vorticity distribution. From the linearity of the problem, one may separate the two different sources (see e.g. Newman Reference Newman1977; Fernandez-Feria & Alaminos-Quesada Reference Fernandez-Feria and Alaminos-Quesada2018) to obtain the following two integral equations for $\unicode[STIX]{x1D71B}_{0}$ and $\unicode[STIX]{x1D71B}_{se}$:

(3.18)$$\begin{eqnarray}\displaystyle & \displaystyle v_{0}(x,t)=\frac{1}{2\unicode[STIX]{x03C0}}\unicode[STIX]{x2A0D}_{-1}^{1}\frac{\unicode[STIX]{x1D71B}_{0}(\unicode[STIX]{x1D709},t)}{\unicode[STIX]{x1D709}-x}\,\text{d}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

(3.19)$$\begin{eqnarray}\displaystyle & \displaystyle -\frac{1}{2\unicode[STIX]{x03C0}}\int _{1}^{\infty }\frac{\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)}{\unicode[STIX]{x1D709}-x}\,\text{d}\unicode[STIX]{x1D709}=\frac{1}{2\unicode[STIX]{x03C0}}\unicode[STIX]{x2A0D}_{-1}^{1}\frac{\unicode[STIX]{x1D71B}_{se}(\unicode[STIX]{x1D709},t)}{\unicode[STIX]{x1D709}-x}\,\text{d}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x2A0D}$ denotes Cauchy’s principal value of the integral (e.g. Butkov Reference Butkov1968, chap. 2) and $v_{0}$ is given by (2.9). The solutions of these singular, linear integral equations of the first kind with constant integration limits are given by (Polyanin & Manzhirov Reference Polyanin and Manzhirov1998)

(3.20)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D71B}_{0}(x,t)=\frac{1}{\sqrt{1-x^{2}}}\left\{\frac{\unicode[STIX]{x1D6E4}_{0}(t)}{\unicode[STIX]{x03C0}}-\frac{2}{\unicode[STIX]{x03C0}}\unicode[STIX]{x2A0D}_{-1}^{1}\frac{\sqrt{1-\unicode[STIX]{x1D709}^{2}}}{\unicode[STIX]{x1D709}-x}\left[{\mathcal{A}}(t)+{\mathcal{B}}(t)\unicode[STIX]{x1D709}+{\mathcal{C}}(t)\unicode[STIX]{x1D709}^{2}\right]\text{e}^{\mathsf{b}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}\right\}\!,\qquad & \displaystyle\end{eqnarray}$$

(3.21)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D71B}_{se}(x,t)=\frac{1}{\unicode[STIX]{x03C0}}\sqrt{\frac{1-x}{1+x}}\int _{1}^{\infty }\sqrt{\frac{\unicode[STIX]{x1D709}+1}{\unicode[STIX]{x1D709}-1}}\frac{\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)}{\unicode[STIX]{x1D709}-x}\,\text{d}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$

where the regularity of $\unicode[STIX]{x1D71B}_{s}$ at its corresponding trailing edge $x=1$, or Kutta condition, has been applied. In addition, substituting into (3.15) and (3.17) one obtains

(3.22)$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6E4}_{0}(t) & = & \displaystyle -2\unicode[STIX]{x03C0}\left\{{\mathcal{A}}(t)\text{I}_{0}(\mathsf{b})+\left[{\mathcal{A}}(t)+\left(1+\frac{1}{\mathsf{b}}\right)({\mathcal{B}}(t)+{\mathcal{C}}(t))\right]\text{I}_{1}(\mathsf{b})\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left[{\mathcal{B}}(t)+{\mathcal{C}}(t)\left(1-\frac{1}{\mathsf{b}}\right)\right]\text{I}_{2}(\mathsf{b})\right\},\end{eqnarray}$$

(3.23)$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6E4}_{se}(t)=\int _{1}^{\infty }\left(\sqrt{\frac{\unicode[STIX]{x1D709}+1}{\unicode[STIX]{x1D709}-1}}-1\right)\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)\,\text{d}\unicode[STIX]{x1D709}, & & \displaystyle\end{eqnarray}$$

where $\text{I}_{n}(\mathsf{b})$, $n=0,1,2$, is the modified Bessel function of the first kind and order $n$ (Olver & Maximon Reference Olver, Maximon, Olver, Lozier, Boisvert and Clark2010) applied to the complex number $\mathsf{b}$.

Finally, from Kelvin’s theorem (3.16), and considering that the vorticity in the wake is convected downstream with velocity $U$, so it remains constant in a reference frame moving with the fluid,

(3.24a-c)$$\begin{eqnarray}\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)=\unicode[STIX]{x1D71B}_{e}(X),\quad z_{e}(\unicode[STIX]{x1D709},t)=z_{e}(X),\quad X=\unicode[STIX]{x1D709}-Ut,\end{eqnarray}$$

the wake vorticity distribution $\unicode[STIX]{x1D71B}_{e}$ is given by the well-known solution in terms of $\unicode[STIX]{x1D6E4}_{0}$ (Theodorsen Reference Theodorsen1935; von Kármán & Sears Reference von Kármán and Sears1938)

(3.25)$$\begin{eqnarray}\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)=\frac{2\unicode[STIX]{x1D6E4}_{0}(t)}{\unicode[STIX]{x03C0}}\frac{\text{e}^{-\text{i}k\unicode[STIX]{x1D709}}}{\text{i}\text{H}_{0}^{(2)}(k)+\text{H}_{1}^{(2)}(k)},\end{eqnarray}$$

where $\text{H}_{n}^{(2)}(z)=\text{J}_{n}(z)-\text{i}\text{Y}_{n}(z)$, $n=0,1$, is the Hankel function of the second kind and order $n$, related to the Bessel functions of the first and second kind $\text{J}_{n}(z)$ and $\text{Y}_{n}(z)$ (Olver & Maximon Reference Olver, Maximon, Olver, Lozier, Boisvert and Clark2010).

3.2 Lift, thrust, moment and input power

Taking into account (3.24), the temporal derivatives of $\unicode[STIX]{x1D71B}_{e}$ and $z_{e}$ in the general expressions for $D$, $L$ and $M$ can be easily computed using Leibniz’ rule. For the lift one obtains, in non-dimensional form,

(3.26)$$\begin{eqnarray}C_{L}=\frac{L}{\frac{1}{2}\unicode[STIX]{x1D70C}U^{2}c}=C_{L0}+C_{L1}+C_{L2},\end{eqnarray}$$

where

(3.27a-c)$$\begin{eqnarray}C_{L0}=\unicode[STIX]{x1D6E4}_{0},\quad C_{L1}=-\frac{\text{d}}{\text{d}t}\int _{-1}^{1}x\unicode[STIX]{x1D71B}_{0}(x,t)\,\text{d}x,\quad C_{L2}=\int _{1}^{\infty }\frac{\unicode[STIX]{x1D71B}_{e}(\unicode[STIX]{x1D709},t)}{\sqrt{\unicode[STIX]{x1D709}^{2}-1}}\,\text{d}\unicode[STIX]{x1D709}\end{eqnarray}$$

are the quasi-steady lift, the apparent-mass lift and the lift induced by its own unsteady wake. Solving the integrals,

(3.28)$$\begin{eqnarray}C_{L}=-2\unicode[STIX]{x03C0}\frac{\text{i}k}{\mathsf{b}}\left[({\mathcal{A}}+{\mathcal{C}})\text{I}_{1}(\mathsf{b})+\left({\mathcal{B}}-\frac{3}{\mathsf{b}}{\mathcal{C}}\right)\text{I}_{2}(\mathsf{b})\right]+\unicode[STIX]{x1D6E4}_{0}C(k),\end{eqnarray}$$

where $C(k)$ is the Theodorsen function given by

(3.29)$$\begin{eqnarray}C(k)=\frac{\text{H}_{1}^{(2)}(k)}{\text{i}\text{H}_{0}^{(2)}(k)+\text{H}_{1}^{(2)}(k)}=F(k)+\text{i}G(k).\end{eqnarray}$$

To compute $D$ from (3.4) one assumes that the vertical displacement of the wake coincides with the trailing-edge location at the time $t^{\prime }=t+(1-\unicode[STIX]{x1D709})/U$ when it was shed from the airfoil, $z_{s}(1,t^{\prime })$; i.e.

(3.30)$$\begin{eqnarray}z_{e}(X)=h\left(\frac{1-X}{U}\right)-(1-a)\unicode[STIX]{x1D6FC}\left(\frac{1-X}{U}\right)+(1-p)^{2}\unicode[STIX]{x1D6FF}\left(\frac{1-X}{U}\right).\end{eqnarray}$$

Thus,

(3.31)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int _{1}^{\infty }z_{e}\unicode[STIX]{x1D71B}_{e}\,\text{d}\unicode[STIX]{x1D709}=z_{s}(1,t)\unicode[STIX]{x1D71B}_{e}(1,t).\end{eqnarray}$$

Therefore, the thrust, or minus the drag (3.4), in non-dimensional form can be written as

(3.32)$$\begin{eqnarray}\displaystyle C_{T} & = & \displaystyle \frac{T}{\frac{1}{2}\unicode[STIX]{x1D70C}U^{2}c}=-\frac{\text{d}}{\text{d}t}\{\!\text{Re}[\unicode[STIX]{x1D6E4}_{0}]\times \text{Re}[{\mathcal{Q}}({\mathcal{F}},{\mathcal{E}},{\mathcal{D}})]+\text{Re}[{\mathcal{A}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{0}({\mathcal{F}},{\mathcal{E}},{\mathcal{D}})]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[{\mathcal{B}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{1}({\mathcal{F}},{\mathcal{E}},{\mathcal{D}})]+\text{Re}[{\mathcal{C}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{2}({\mathcal{F}},{\mathcal{E}},{\mathcal{D}})]\!\}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{2k}{\unicode[STIX]{x03C0}}\text{Re}[\unicode[STIX]{x1D6E4}_{0}C_{1}(k)]\times \text{Re}\{[{\mathcal{F}}+{\mathcal{E}}+{\mathcal{D}}]\text{e}^{\mathsf{b}}\},\end{eqnarray}$$

where the functions ${\mathcal{Q}}$, $\unicode[STIX]{x1D6FA}_{n}$ and $C_{1}(k)$ are defined as

(3.33)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{Q}}(\tilde{a},\tilde{b},\tilde{c})=\tilde{a}(\text{I}_{0}(\mathsf{b})+{\mathcal{J}}_{0})+\tilde{b}(\text{I}_{1}(\mathsf{b})+{\mathcal{J}}_{1})+\tilde{c}\left(\frac{\text{I}_{1}(\mathsf{b})}{\mathsf{b}}+I_{2}(\mathsf{b})+{\mathcal{J}}_{2}\right), & \displaystyle\end{eqnarray}$$

(3.34)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{n}(\tilde{a},\tilde{b},\tilde{c})=\tilde{a}{\mathcal{I}}_{0,n}+\tilde{b}{\mathcal{I}}_{1,n}+\tilde{c}{\mathcal{I}}_{2,n},\quad n=0,1,2, & \displaystyle\end{eqnarray}$$

(3.35)$$\begin{eqnarray}\displaystyle & \displaystyle C_{1}(k)=\frac{\frac{1}{k}\text{e}^{-\text{i}k}}{\text{i}\text{H}_{0}^{(2)}(k)+\text{H}_{1}^{(2)}(k)}=F_{1}(k)+\text{i}G_{1}(k), & \displaystyle\end{eqnarray}$$

and the functions ${\mathcal{J}}_{n}$ and ${\mathcal{I}}_{n,m}$ are given by

(3.36)$$\begin{eqnarray}\displaystyle \displaystyle {\mathcal{J}}_{n}(\mathsf{b},k) & \equiv & \displaystyle \frac{2/\unicode[STIX]{x03C0}^{2}}{\text{i}\text{H}_{0}^{(2)}(k)+\text{H}_{1}^{(2)}(k)}\int _{-1}^{1}x^{n}\text{e}^{\mathsf{b}x}\sqrt{\frac{1-x}{1+x}}\nonumber\\ \displaystyle & & \displaystyle \times \left(\int _{1}^{\infty }\sqrt{\frac{\unicode[STIX]{x1D709}+1}{\unicode[STIX]{x1D709}-1}}\frac{\text{e}^{-\text{i}k\unicode[STIX]{x1D709}}}{\unicode[STIX]{x1D709}-x}\,\text{d}\unicode[STIX]{x1D709}\right)\,\text{d}x={\mathcal{J}}_{n}^{\text{R}}+\text{i}{\mathcal{J}}_{n}^{\text{I}},\end{eqnarray}$$

(3.37)$$\begin{eqnarray}\displaystyle \displaystyle {\mathcal{I}}_{n,m}(\mathsf{b}) & \equiv & \displaystyle -\frac{2}{\unicode[STIX]{x03C0}}\int _{-1}^{1}\frac{x^{n}\text{e}^{\mathsf{b}x}}{\sqrt{1-x^{2}}}\left(\int _{-1}^{1}\frac{\sqrt{1-\unicode[STIX]{x1D709}^{2}}}{\unicode[STIX]{x1D709}-x}\unicode[STIX]{x1D709}^{m}\text{e}^{\mathsf{b}\unicode[STIX]{x1D709}}\,\text{d}\unicode[STIX]{x1D709}\right)\text{d}x\nonumber\\ \displaystyle & = & \displaystyle {\mathcal{I}}_{n,m}^{\text{R}}+\text{i}{\mathcal{I}}_{n,m}^{\text{I}},\quad n,m=0,1,2,\end{eqnarray}$$

where, in general, these last integrals have to be solved numerically. It should be noted that ${\mathcal{I}}_{1,1}={\mathcal{I}}_{2,2}=0$ and ${\mathcal{I}}_{0,0}={\mathcal{I}}_{2,0}$. In addition, when the growth factor of the undulatory motion vanishes ($b_{1}=0$) then

(3.38a-c)$$\begin{eqnarray}{\mathcal{I}}_{0,0}^{\text{R}}={\mathcal{I}}_{0,1}^{\text{I}}={\mathcal{I}}_{0,2}^{\text{R}}=0,\quad {\mathcal{I}}_{1,0}^{\text{I}}={\mathcal{I}}_{1,1}={\mathcal{I}}_{1,2}^{\text{I}}=0,\quad {\mathcal{I}}_{2,0}^{\text{R}}={\mathcal{I}}_{2,1}^{\text{I}}={\mathcal{I}}_{2,2}=0,\end{eqnarray}$$

where the superscripts R and I denote real and imaginary parts, respectively. As discussed in Fernandez-Feria (Reference Fernandez-Feria2016) for a rigid foil, this impulse formulation includes the leading-edge suction force together with all the other contributions to the thrust force coming from the vorticity distributions on the airfoil and the vortex wake.

The moment (3.8) can be written, in non-dimensional form, as

(3.39)$$\begin{eqnarray}C_{M}=\frac{M}{\frac{1}{2}\unicode[STIX]{x1D70C}U^{2}c^{2}}=C_{M0}+C_{M1}+C_{M2}+\frac{a}{2}C_{L},\end{eqnarray}$$

where

(3.40a-c)$$\begin{eqnarray}C_{M0}=-\frac{1}{2}\int _{-1}^{1}x\unicode[STIX]{x1D71B}_{0}(x,t)\,\text{d}x,\quad C_{M1}=\frac{1}{4}\frac{\text{d}}{\text{d}t}\int _{-1}^{1}\left(x^{2}-\frac{1}{2}\right)\unicode[STIX]{x1D71B}_{0}(x,t)\,\text{d}x,\quad C_{M2}=\frac{C_{L2}}{4}.\end{eqnarray}$$

Solving the integrals,

(3.41)$$\begin{eqnarray}\displaystyle C_{M} & = & \displaystyle \frac{\unicode[STIX]{x1D6E4}_{0}}{2}\left[C(k)\left(\frac{1}{2}+a\right)-1\right]-\frac{\unicode[STIX]{x03C0}}{\mathsf{b}}(\text{i}ka+1)\left[({\mathcal{A}}+{\mathcal{C}})\text{I}_{1}(\mathsf{b})+\left({\mathcal{B}}-\frac{3}{\mathsf{b}}{\mathcal{C}}\right)\text{I}_{2}(\mathsf{b})\right]\nonumber\\ \displaystyle & & \displaystyle +\frac{\unicode[STIX]{x03C0}\text{i}k}{2\mathsf{b}^{2}}\{[{\mathcal{B}}+\mathsf{b}({\mathcal{A}}+{\mathcal{C}})]\text{I}_{2}(\mathsf{b})+[\mathsf{b}{\mathcal{B}}-3{\mathcal{C}}]\text{I}_{3}(\mathsf{b})\}.\end{eqnarray}$$

Finally, the input power in non-dimensional form is given by

(3.42)$$\begin{eqnarray}\displaystyle C_{P} & = & \displaystyle \frac{P}{\frac{1}{2}\unicode[STIX]{x1D70C}U^{3}c}=-\text{Re}[\unicode[STIX]{x1D6E4}_{0}]\times \text{Re}[{\mathcal{Q}}(\dot{{\mathcal{F}}},\dot{{\mathcal{E}}},\dot{{\mathcal{D}}})]-\text{Re}[{\mathcal{H}}]\times \text{Re}[\unicode[STIX]{x1D6E4}_{0}C_{1}(k)]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[\dot{\unicode[STIX]{x1D6E4}}_{0}]\times \text{Re}[{\mathcal{Q}}(g_{1},g_{2},g_{3})]-\text{Re}[{\mathcal{A}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{0}(\dot{{\mathcal{F}}},\dot{{\mathcal{E}}},\dot{{\mathcal{D}}})]\nonumber\\ \displaystyle & & \displaystyle -\,\text{Re}[{\mathcal{B}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{1}(\dot{{\mathcal{F}}},\dot{{\mathcal{E}}},\dot{{\mathcal{D}}})]-\text{Re}[{\mathcal{C}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{2}(\dot{{\mathcal{F}}},\dot{{\mathcal{E}}},\dot{{\mathcal{D}}})]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[\dot{{\mathcal{A}}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{0}(g_{1},g_{2},g_{3})]+\text{Re}[\dot{{\mathcal{B}}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{1}(g_{1},g_{2},g_{3})]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[\dot{{\mathcal{C}}}]\times \text{Re}[\unicode[STIX]{x1D6FA}_{2}(g_{1},g_{2},g_{3})],\end{eqnarray}$$

where a dot denotes the time derivative, and the functions ${\mathcal{H}}$, $g_{1},g_{2}$ and $g_{3}$ are defined as

(3.43)$$\begin{eqnarray}{\mathcal{H}}=\frac{\text{e}^{\mathsf{b}}}{\mathsf{b}^{2}}\left\{\left[2\left(1-\frac{1}{\mathsf{b}}\right)-\mathsf{b}\right]{\mathcal{D}}+(1-\mathsf{b}){\mathcal{E}}-\mathsf{b}{\mathcal{F}}\right\},\end{eqnarray}$$

(3.44a-c)$$\begin{eqnarray}g_{1}=\frac{1}{\mathsf{b}^{3}}\left(2{\mathcal{D}}+\mathsf{b}(\mathsf{b}{\mathcal{F}}-{\mathcal{E}})\right),\quad g_{2}=\frac{1}{\mathsf{b}^{2}}(\mathsf{b}{\mathcal{E}}-2{\mathcal{D}}),\quad g_{3}=\frac{{\mathcal{D}}}{\mathsf{b}}.\end{eqnarray}$$

3.3 Propulsive efficiency

Propulsive efficiency is defined as the ratio between the time-averaged power output of the airfoil (thrust force multiplied by the forward speed $U$) and the time-averaged input power required to drive the airfoil:

(3.45)$$\begin{eqnarray}\unicode[STIX]{x1D702}=\frac{\bar{C}_{T}}{\bar{C}_{P}},\end{eqnarray}$$

where the time-averaged quantities are defined as

(3.46a,b)$$\begin{eqnarray}\bar{C}_{T}\equiv \frac{1}{T}\int _{t}^{t+T}C_{T}(t)\,\text{d}t,\quad \bar{C}_{P}\equiv \frac{1}{T}\int _{t}^{t+T}C_{P}(t)\,\text{d}t,\end{eqnarray}$$

with $T=2\unicode[STIX]{x03C0}/k$ being the non-dimensional period of the oscillation. This quantity will be computed explicitly for the particular case discussed below.

4 First oscillation mode of a flexible plate

The special case with $\mathsf{b}=0$ corresponds to a quadratic flexural motion of a heaving and pitching foil, which is of particular interest because it is equivalent to the first oscillation mode of a flexible foil. The corresponding expressions for the lift, thrust, moment and input power are obtained by making the limit $\mathsf{b}\rightarrow 0$ in the above expressions.

For the lift one obtains

(4.1)$$\begin{eqnarray}C_{L}=\unicode[STIX]{x03C0}\left[\dot{\unicode[STIX]{x1D6FC}}+2p\dot{\unicode[STIX]{x1D6FF}}-\ddot{h}-a\ddot{\unicode[STIX]{x1D6FC}}-\left(p^{2}+{\textstyle \frac{1}{4}}\right)\ddot{\unicode[STIX]{x1D6FF}}\right]+\unicode[STIX]{x1D6E4}_{0}C(k).\end{eqnarray}$$

When $\unicode[STIX]{x1D6FF}_{m}=0$ this expression coincides with the lift coefficient obtained by von Kármán & Sears (Reference von Kármán and Sears1938) using the present impulse theory, and previously by Theodorsen (Reference Theodorsen1935) from a more standard potential flow approach.

For the thrust, the functions ${\mathcal{Q}}$ and $\unicode[STIX]{x1D6FA}_{n}$ become

(4.2)$$\begin{eqnarray}\lim _{\mathsf{b}\rightarrow 0}{\mathcal{Q}}(\tilde{a},\tilde{b},\tilde{c})=\frac{\text{i}\tilde{b}}{k}C(k)+\left[\left(\frac{1+\text{i}k}{k}\right)\tilde{b}+\text{i}\tilde{a}+\left(\frac{2+\text{i}k}{k}-\frac{2\text{i}}{k^{2}}\right)\tilde{c}\right]\frac{2}{\unicode[STIX]{x03C0}}C_{1}(k)+\frac{\tilde{c}}{k}C_{2}(k),\end{eqnarray}$$

(4.3a-c)$$\begin{eqnarray}\lim _{\mathsf{b}\rightarrow 0}\unicode[STIX]{x1D6FA}_{0}(\tilde{a},\tilde{b},\tilde{c})=\tilde{b}\unicode[STIX]{x03C0},\quad \lim _{\mathsf{b}\rightarrow 0}\unicode[STIX]{x1D6FA}_{1}(\tilde{a},\tilde{b},\tilde{c})=\tilde{c}\frac{\unicode[STIX]{x03C0}}{4},\quad \lim _{\mathsf{b}\rightarrow 0}\unicode[STIX]{x1D6FA}_{2}(\tilde{a},\tilde{b},\tilde{c})=\tilde{b}\frac{\unicode[STIX]{x03C0}}{4},\end{eqnarray}$$

where $C_{2}(k)$ is given by

(4.4)$$\begin{eqnarray}C_{2}(k)\equiv \frac{\text{H}_{2}^{(2)}(k)}{\text{i}\text{H}_{0}^{(2)}(k)+\text{H}_{1}^{(2)}(k)}=F_{2}(k)+\text{i}G_{2}(k).\end{eqnarray}$$

Consequently, the thrust coefficient can be written as

(4.5)$$\begin{eqnarray}\displaystyle C_{T} & = & \displaystyle -\text{Re}\left[\unicode[STIX]{x1D6FC}+2p\unicode[STIX]{x1D6FF}\right]\times C_{L}+\frac{\unicode[STIX]{x03C0}}{2}\left\{\text{Re}[\unicode[STIX]{x1D6FF}]\times \text{Re}\left[\frac{\ddot{\unicode[STIX]{x1D6FC}}}{2}+p\ddot{\unicode[STIX]{x1D6FF}}-\dot{\unicode[STIX]{x1D6FF}}\right]\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\text{Re}[\dot{\unicode[STIX]{x1D6FF}}]\times \text{Re}\left[\frac{\dot{\unicode[STIX]{x1D6FC}}}{2}+p\dot{\unicode[STIX]{x1D6FF}}-\unicode[STIX]{x1D6FF}\right]\right\}+\text{Re}\left[\dot{\unicode[STIX]{x1D6FC}}+2p\dot{\unicode[STIX]{x1D6FF}}\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\text{Re}\left\{\unicode[STIX]{x03C0}\left[{\dot{h}}+a\dot{\unicode[STIX]{x1D6FC}}+\left(p^{2}+\frac{1}{4}\right)\dot{\unicode[STIX]{x1D6FF}}-\left(\unicode[STIX]{x1D6FC}+2p\unicode[STIX]{x1D6FF}\right)\right]+\unicode[STIX]{x1D6E4}_{0}\left[\frac{\text{i}}{k}C+\left(\frac{1+\text{i}k}{k}\right)\frac{2}{\unicode[STIX]{x03C0}}C_{1}\right]\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}\left[{\dot{h}}+a\dot{\unicode[STIX]{x1D6FC}}-\left(\unicode[STIX]{x1D6FC}+2p\unicode[STIX]{x1D6FF}\right)+p^{2}\dot{\unicode[STIX]{x1D6FF}}\right]\times \text{Re}\left[\unicode[STIX]{x1D6E4}_{0}\frac{-2\text{i}}{\unicode[STIX]{x03C0}}C_{1}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[\unicode[STIX]{x1D6FF}]\times \text{Re}\left[\unicode[STIX]{x1D6E4}_{0}\left(-\text{i}C_{2}-2\frac{1+\text{i}k}{k}\frac{2}{\unicode[STIX]{x03C0}}C_{1}\right)\right]\nonumber\\ \displaystyle & & \displaystyle +\,\text{Re}[\dot{\unicode[STIX]{x1D6FF}}]\times \text{Re}\left\{\unicode[STIX]{x1D6E4}_{0}\left[\left(\frac{2\text{i}}{k^{2}}-\frac{2+\text{i}k}{k}\right)\frac{2}{\unicode[STIX]{x03C0}}C_{1}-\frac{C_{2}}{k}\right]\right\}.\end{eqnarray}$$

The complex form of this expression coincides with the thrust obtained by Fernandez-Feria (Reference Fernandez-Feria2016) when $\unicode[STIX]{x1D6FF}_{m}=0$, though here the real part of the nonlinear expression is separated in a slightly different form.

The moment coefficient becomes

(4.6)$$\begin{eqnarray}\displaystyle \displaystyle C_{M} & = & \displaystyle \frac{\unicode[STIX]{x1D6E4}_{0}}{2}\left(\frac{1}{2}+a\right)C(k)-\frac{\unicode[STIX]{x03C0}}{2}\left\{\left(\frac{1}{2}-a\right)\dot{\unicode[STIX]{x1D6FC}}+a\ddot{h}+\left(\frac{1}{8}+a^{2}\right)\ddot{\unicode[STIX]{x1D6FC}}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\displaystyle \left[\frac{p}{4}+a\left(p^{2}+\frac{1}{4}\right)\right]\ddot{\unicode[STIX]{x1D6FF}}+\left[\left(\frac{1}{2}-a\right)p-1\right]\frac{\dot{\unicode[STIX]{x1D6FF}}}{2}-\unicode[STIX]{x1D6FF}\right\}.\end{eqnarray}$$

As in the case of the lift, the above expression coincides with the moment obtained by von Kármán & Sears (Reference von Kármán and Sears1938) when $\unicode[STIX]{x1D6FF}_{m}=0$.

Finally, the input power coefficient can be written in the physically relevant form

(4.7)$$\begin{eqnarray}C_{P}=-{\dot{h}}(t)C_{L}-2\dot{\unicode[STIX]{x1D6FC}}(t)C_{M}-\dot{\unicode[STIX]{x1D6FF}}(t)\unicode[STIX]{x1D6F6},\end{eqnarray}$$

where

(4.8)$$\begin{eqnarray}\unicode[STIX]{x1D6F6}=p(4C_{Mp}-pC_{L})+\unicode[STIX]{x1D6E4}_{0}C(k)-\frac{\unicode[STIX]{x03C0}}{4}\left[\ddot{h}+a\ddot{\unicode[STIX]{x1D6FC}}+\left(p^{2}+\frac{1}{3}\right)\ddot{\unicode[STIX]{x1D6FF}}-2\unicode[STIX]{x1D6FF}\right],\end{eqnarray}$$

and $C_{Mp}\equiv C_{M}(a=p)$ is the moment coefficient but in relation to the point $x=p$, where the flexure component of the motion vanishes. It must be noted that, to obtain the above input power coefficient, the limit $\mathsf{b}\rightarrow 0$ in the double integral appearing in (3.13), which now can be written as

(4.9)$$\begin{eqnarray}\int _{-1}^{1}\left(\int _{x}^{1}\frac{\unicode[STIX]{x2202}z_{s}}{\unicode[STIX]{x2202}t}\,\text{d}\unicode[STIX]{x1D709}\right)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D71B}_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x=\int _{-1}^{1}[{\mathcal{H}}+(g_{1}+g_{2}x+g_{3}x^{2})\text{e}^{\mathsf{b}x}]\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D71B}_{s}}{\unicode[STIX]{x2202}t}\,\text{d}x,\end{eqnarray}$$

has to be done collectively to all the terms inside the integrand, because the functions ${\mathcal{H}},g_{1},g_{2}$ and $g_{3}$ diverge separately as $\mathsf{b}\rightarrow 0$, but together

(4.10)$$\begin{eqnarray}\lim _{\mathsf{b}\rightarrow 0}[{\mathcal{H}}+(g_{1}+g_{2}x+g_{3}x^{2})\text{e}^{\mathsf{b}x}]={\mathcal{F}}(x-1)+\frac{{\mathcal{E}}}{2}(x^{2}-1)+\frac{{\mathcal{D}}}{3}(x^{3}-1).\end{eqnarray}$$

Thus, the general expression (3.42) is valid when $\mathsf{b}\neq 0$, and one has to use expression (4.7) when $\mathsf{b}=0$.

4.1 Time-averaged coefficients and propulsive efficiency

The time-averaged coefficients (3.46) can now be written in simple closed forms. To that end it is convenient to define the following non-dimensional parameters to separate the different components of the foil’s motion:

(4.11a-c)$$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{a_{0}}{kh_{0}},\quad \unicode[STIX]{x1D703}_{hd}=\frac{\unicode[STIX]{x1D6FF}_{m}}{kh_{0}},\quad \unicode[STIX]{x1D703}_{pd}=\frac{\unicode[STIX]{x1D6FF}_{m}}{a_{0}},\end{eqnarray}$$

where the first one, $\unicode[STIX]{x1D703}$, is the well-known feathering parameter (Lighthill Reference Lighthill1969) for a pitching and heaving motion of a rigid foil. Thus, the time-averaged thrust coefficient can be written either normalized in relation to a pure heaving motion,

(4.12)$$\begin{eqnarray}\displaystyle {\hat{C}}_{Th}\equiv \frac{\bar{C}_{T}}{(kh_{0})^{2}} & = & \displaystyle t_{h}(k)+t_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}+t_{p}(k,a)\unicode[STIX]{x1D703}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,t_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{hd}+t_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{hd}\unicode[STIX]{x1D703}+t_{d}(k,p)\unicode[STIX]{x1D703}_{hd}^{2},\end{eqnarray}$$

or in relation to a pure pitching motion,

(4.13)$$\begin{eqnarray}\displaystyle \displaystyle {\hat{C}}_{Tp}\equiv \frac{\bar{C}_{T}}{a_{0}^{2}} & = & \displaystyle t_{h}(k)\unicode[STIX]{x1D703}^{-2}+t_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}^{-1}+t_{p}(k,a)\nonumber\\ \displaystyle & & \displaystyle +\,t_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{pd}\unicode[STIX]{x1D703}^{-1}+t_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{pd}+t_{d}(k,p)\unicode[STIX]{x1D703}_{pd}^{2},\end{eqnarray}$$

where the functions $t_{h}(k)$, $t_{hp}(k,a,\unicode[STIX]{x1D719})$, $t_{p}(k,a)$, $t_{hd}(k,p,\unicode[STIX]{x1D713})$, $t_{pd}(k,a,p,\unicode[STIX]{x1D719},\unicode[STIX]{x1D713})$ and $t_{d}(k,p)$ are given by (A 1)–(A 6) in appendix A. Note that for $\unicode[STIX]{x1D703}_{hd}=0$ or $\unicode[STIX]{x1D703}_{pd}=0$, these expressions do not coincide exactly with those for an oscillating rigid foil given in Fernandez-Feria (Reference Fernandez-Feria2017) because here the real part of the complex nonlinear expression for the thrust force is separated in a slightly different way, but the results are practically indistinguishable.

The time-averaged input power coefficient is given by

(4.14)$$\begin{eqnarray}\bar{C}_{P}=-\overline{{\dot{h}}(t)C_{L}}-2\overline{\dot{\unicode[STIX]{x1D6FC}}(t)C_{M}}-\overline{\dot{\unicode[STIX]{x1D6FF}}(t)\unicode[STIX]{x1D6F6}},\end{eqnarray}$$

which, similarly to the mean thrust coefficient, can be written in terms of the parameters (4.11) either as

(4.15)$$\begin{eqnarray}\displaystyle {\hat{C}}_{Ph}\equiv \frac{\bar{C}_{P}}{(kh_{0})^{2}} & = & \displaystyle p_{h}(k)+p_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}+p_{p}(k,a)\unicode[STIX]{x1D703}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,p_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{hd}+p_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{hd}\unicode[STIX]{x1D703}+p_{d}(k,p)\unicode[STIX]{x1D703}_{hd}^{2}\end{eqnarray}$$

or as

(4.16)$$\begin{eqnarray}\displaystyle {\hat{C}}_{Pp}\equiv \frac{\bar{C}_{P}}{a_{0}^{2}} & = & \displaystyle p_{h}(k)\unicode[STIX]{x1D703}^{-2}+p_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}^{-1}+p_{p}(k,a)\nonumber\\ \displaystyle & & \displaystyle +\,p_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{pd}\unicode[STIX]{x1D703}^{-1}+p_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{pd}+p_{d}(k,p)\unicode[STIX]{x1D703}_{pd}^{2},\end{eqnarray}$$

where the functions $p_{h}(k)$, $p_{hp}(k,a,\unicode[STIX]{x1D719})$, $p_{p}(k,a)$, $p_{hd}(k,p,\unicode[STIX]{x1D713})$, $p_{pd}(k,a,p,\unicode[STIX]{x1D719},\unicode[STIX]{x1D713})$ and $p_{d}(k,p)$ are given by (A 7)–(A 12) in appendix A.

Finally, the propulsive efficiency (3.45) can be computed using any of the above alternative forms of $\bar{C}_{T}$ and $\bar{C}_{P}$. It is convenient to use a propulsive efficiency relative to that of a rigid oscillating foil, $\unicode[STIX]{x1D702}_{0}$:

(4.17)$$\begin{eqnarray}\hat{\unicode[STIX]{x1D702}}\equiv \unicode[STIX]{x1D702}-\unicode[STIX]{x1D702}_{0},\quad \text{with }\unicode[STIX]{x1D702}_{0}=\frac{t_{h}(k)+t_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}+t_{p}(k,a)\unicode[STIX]{x1D703}^{2}}{p_{h}(k)+p_{hp}(k,a,\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}+p_{p}(k,a)\unicode[STIX]{x1D703}^{2}}.\end{eqnarray}$$

5 Model validation

Before analysing in detail the results corresponding to the first oscillation mode of a flexible foil ($\mathsf{b}=0$) considered in the previous section, it is convenient to validate the present theoretical results by comparing them with available numerical data for both $\mathsf{b}=0$ and $\mathsf{b}\neq 0$.

We first consider the numerical results of Le et al. (Reference Le, Ko, Byun, Park and Park2010), who investigated the performance of a flapping wing with different chord flexures, corresponding, in the notation of the present model, to $\mathsf{b}=0$ and $a_{0}=0$ (no pitching motion). Those authors reported results for different values of the flexure amplitude $\unicode[STIX]{x1D6FF}_{m}$ and the heave amplitude $h_{0}$, with varying phase shift $\unicode[STIX]{x1D713}$ between the two components of the foil’s motion. In particular, they considered the case of a plunging NACA0012 foil clamped at the leading edge (i.e. with $p=-1$ in our notation and without rotation at the leading edge) for $Re=3\times 10^{4}$.

Figure 3. Comparison of the present results for the mean thrust, input power and efficiency with numerical results from Le et al. (Reference Le, Ko, Byun, Park and Park2010) for a pure heaving motion with deflection $(\mathsf{b}=a_{0}=0)$, with $h_{0}=0.35,k=1.82,p=-1,C_{T0}=-0.055$: (a) $\unicode[STIX]{x1D6FF}_{m}=0.1$, $\unicode[STIX]{x0394}t/T=-0.06$ for $\bar{C}_{P}$; (b) $\unicode[STIX]{x1D6FF}_{m}=0.4$, $\unicode[STIX]{x0394}t/T=-0.025$ for $\bar{C}_{P}$.

Figure 3 compares the numerical results of Le et al. (Reference Le, Ko, Byun, Park and Park2010) for the mean values of the thrust and input power coefficients, and for the propulsive efficiency, as functions of the phase shift $\unicode[STIX]{x1D713}$ between the heaving and deflection motions, with the present theoretical results for two values of $\unicode[STIX]{x1D6FF}_{m}$. For reference, the results for a rigid foil are also included, which are obviously independent of $\unicode[STIX]{x1D713}$. The agreement is quite good for the three mean magnitudes represented in all the range of values of $\unicode[STIX]{x1D713}$. But it should be noted that, since the computation of the input power is very sensitive to any temporal shift between the lift coefficient $C_{L}(t)$ obtained numerically and the vertical motion ${\dot{h}}(t)$, to fit $\bar{C}_{P}$ we have included a very small temporal shift $\unicode[STIX]{x0394}t$ in the lift when computing the mean value of $C_{P}$ from our theoretical formulation as follows:

(5.1)$$\begin{eqnarray}\bar{C}_{P}=-\frac{1}{T}\int _{t}^{t+T}[{\dot{h}}(t)C_{L}(t+\unicode[STIX]{x0394}t)+\dot{\unicode[STIX]{x1D6FF}}(t)\unicode[STIX]{x1D6F6}(t)]\,\text{d}t.\end{eqnarray}$$

The values of $\unicode[STIX]{x0394}t$ are given in the caption of figure 3, together with all the kinematic parameters. In figure 3(b), $\unicode[STIX]{x0394}t$ is taken directly from $C_{L}(t)$ data in Le et al. (Reference Le, Ko, Byun, Park and Park2010), which are not available for the case of figure 3(a), and we take the values that best fit $C_{P}$. These $\unicode[STIX]{x0394}t$ are very small compared with the period $T$, justifying their use to compensate any small temporal shift in the computation of $C_{L}$ from the numerical simulations. Note also the excellent agreement between the theoretical results for $\bar{C_{T}}$ and the numerical results obtained by Le et al. (Reference Le, Ko, Byun, Park and Park2010) for any value of $\unicode[STIX]{x1D713}$, particularly in figure 3(a) where the non-dimensional flexure amplitude $\unicode[STIX]{x1D6FF}_{m}$ is smaller.

Next we consider the numerical results of Zhang et al. (Reference Zhang, Pan, Chao and Zhang2018) and of Dong & Lu (Reference Dong and Lu2007), who analysed two different undulatory motions, i.e. with $\mathsf{b}\neq 0$ in the present notation. As commented on in the Introduction, this kind of motion is commonly used to model fish locomotion. To put them into a biological context, the body or caudal fin of fishes can be classified into five groups that differ in the fraction of the body that is displaced laterally (coordinate $z$ in the present notation) (Breder Reference Breder1926): anguilliform (e.g. eels), sub-carangiform (e.g. trout), carangiform (e.g. jack mackerels), thunniform (e.g. tunas) and ostraciiform (e.g. boxfishes). Examples of three of these motions are plotted in figure 4, with the corresponding values of the parameters in the present model given in table 1. Also included are the values of the parameters corresponding to the selected numerical results of Zhang et al. (Reference Zhang, Pan, Chao and Zhang2018) and Dong & Lu (Reference Dong and Lu2007), which are compared below with the present theoretical results, together with other references and some additional information about the numerical simulations.

Figure 4. (a) Anguilliform, carangiform and thunniform motion at the mid-line. (b) Superimposed body outlines separated by $\unicode[STIX]{x0394}t/T=0.25$. The corresponding values of the parameters of the present model and references are given in table 1.

Table 1. Values of the kinematic parameters for the different motions considered in figures 4–6. Anguilliform from Tytell (Reference Tytell2004) and Zhang et al. (Reference Zhang, Pan, Chao and Zhang2018), respectively; carangiform from Dong & Lu (Reference Dong and Lu2007) and Zhang & Eldredge (Reference Zhang and Eldredge2009), respectively; and thunniform from Chang, Zhang & He (Reference Chang, Zhang and He2012).

Figure 5. Anguilliform motion. Time-dependent lift coefficient during one cycle for$k=\unicode[STIX]{x03C0}$. Numerical results from figure 9(d) in Zhang et al. (Reference Zhang, Pan, Chao and Zhang2018) for $Re=5\times 10^{4}$. Note that those authors non-dimensionalize the lift as $L/(\unicode[STIX]{x1D70C}U^{2}c^{2})$, so that the present $C_{L}/4$ is plotted. The corresponding values of the kinematic parameters in the present model are given in table 1.

Figure 6. Carangiform motion. (a) Time-averaged thrust coefficient, $\bar{C}_{T}$, as a function of $k$, with two different values of the offset drag: $\bar{C}_{T0}=-0.055$ and $-0.075$. (b) Time-dependent thrust coefficient during one cycle for $k=1.5\unicode[STIX]{x03C0}$. Numerical results from Dong & Lu (Reference Dong and Lu2007). The corresponding values of the kinematic parameters in the present model are given in table 1.

Figure 5 compares the lift coefficient from the numerical results of Zhang et al. (Reference Zhang, Pan, Chao and Zhang2018) for several NACA profiles with the present theoretical results. The numerical results for the different foil thicknesses are all quite close to each other, and with a good agreement with the present zero-thickness theoretical results, with a slightly better agreement for the NACA0016 profile, but for no particular reason.

Figure 6 compares the thrust coefficient computed numerically by Dong & Lu (Reference Dong and Lu2007) with the present theoretical results, both the time-averaged values as a function of the reduced frequency $k$, and its temporal evolution during one cycle for $k=1.5\unicode[STIX]{x03C0}$. As discussed in Fernandez-Feria (Reference Fernandez-Feria2017), the mean values are corrected to account for the viscous drag by adding to the present theoretical results a quasi-static thrust $\bar{C}_{T0}<0$ (i.e. by subtracting an offset quasi-static drag $-\bar{C}_{T0}$), corresponding to the numerical results for $k=0$. Since those authors do not provide numerical results for $k=0$, we have selected two values of $\bar{C}_{T0}$ in figure 6(a), one that yields a better fit for small values of $k$ and the other one for higher values of $k$, though both values of $\bar{C}_{T0}$ are quite close to each other. Notice that no mean drag correction is needed in figure 6(b) since we plot the instantaneous thrust coefficient minus its numerical value at $t=0$, $C_{T}(t=0)$.

Figure 7. Kinematics schemes for pure flexural motion with different values of $p$ (a), pure heaving motion with chordwise deflection for different values of $\unicode[STIX]{x1D713}$ and $p=-1$ (b) and pure pitching motion with chordwise deflection for different values of $\unicode[STIX]{x1D713}$ and $p=-1$ (c). In all cases the same three instants of time are plotted, as indicated.

6 Results and discussion for $\mathsf{b}=0$

Results with no undulatory motion ($\mathsf{b}=0$) are presented and discussed in this section. The kinematics of the three different cases considered are sketched in figure 7.

Figure 8. Contours of the normalized thrust ${\hat{C}}_{Td}$ (a) and efficiency $\unicode[STIX]{x1D702}$ (b) on the $(p,k)$-plane in a pure flexural motion $(\mathsf{b}=h_{0}=a_{0}=0)$.

6.1 Pure flexural motion

We consider first the simplest case of a pure flexural motion of the first oscillating mode, without heaving or pitching motions ($h_{0}=0$ and $a_{0}=0$). Note from (4.12) and (4.15) that in this case only the last terms of the thrust and power coefficients are different from zero, which are quadratic in the trailing-edge deflection amplitude $\unicode[STIX]{x1D6FF}_{m}$. Thus, it is convenient to use the following normalized thrust and power coefficients:

(6.1a,b)$$\begin{eqnarray}{\hat{C}}_{Td}\equiv \frac{\bar{C}_{T}}{\unicode[STIX]{x1D6FF}_{m}^{2}}=t_{d}(k,p),\quad {\hat{C}}_{Pd}\equiv \frac{\bar{C}_{P}}{\unicode[STIX]{x1D6FF}_{m}^{2}}=p_{d}(k,p),\end{eqnarray}$$

where the functions $t_{d}$ and $p_{d}$ are given in appendix A. The efficiency is $\unicode[STIX]{x1D702}=t_{d}(k,p)/p_{d}(k,p)$.

Figure 8 shows the contours of the normalized thrust ${\hat{C}}_{Td}$ and efficiency $\unicode[STIX]{x1D702}$ as $k$ and the location of the deflection pivot point $p$ are varied. It is clear from this figure that no thrust is generated when $k\lesssim 1$ (i.e. ${\hat{C}}_{Td}<0$). For $k\gtrsim 1$, the maximum propulsive efficiency is always reached for a clamped leading edge ($p=-1$). This can be justified by the consideration that in this case all the bending of the foil is useful to generate thrust (see figure 7a). For this reason we only consider the cases with $p=-1$ in the results reported below when a heaving or a pitching motion is added to the flexural motion. The optimum value of the reduced frequency $k$ in terms of propulsive efficiency for $p=-1$ is about $3$.

6.2 Pure heaving motion with chordwise deflection

Next in complexity is the case of a pure heaving motion combined with a deflection motion of the foil, i.e. the first oscillation mode of a flexible airfoil for a pure heaving motion, for which (4.12) and (4.15) become

(6.2)$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{C}}_{Th}=t_{h}(k)+t_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{hd}+t_{d}(k,p)\unicode[STIX]{x1D703}_{hd}^{2}, & \displaystyle\end{eqnarray}$$

(6.3)$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{C}}_{Ph}=p_{h}(k)+p_{hd}(k,p,\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}_{hd}+p_{d}(k,p)\unicode[STIX]{x1D703}_{hd}^{2}, & \displaystyle\end{eqnarray}$$

with the ratio $\unicode[STIX]{x1D703}_{hd}$ between the deflection and heaving amplitudes defined by (4.11) and the different functions given in appendix A. To compare the thrust of the flexible foil with the thrust generated by its rigid counterpart, we define a normalized thrust as

(6.4)$$\begin{eqnarray}{\hat{C}}_{Th}^{\ast }=\frac{{\hat{C}}_{Th}}{t_{h}(k)},\end{eqnarray}$$

which only depends on the non-dimensional parameters $k,p,\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D703}_{hd}$, as does the normalized efficiency $\hat{\unicode[STIX]{x1D702}}$ defined in (4.17). In order to reduce further the number of parameters we focus on the case of a clamped leading edge, i.e. $p=-1$, which was shown in § 6.1 to yield the best propulsion efficiency.

Figure 9. Contours of the normalized thrust ${\hat{C}}_{Th}^{\ast }$ (a) and efficiency $\hat{\unicode[STIX]{x1D702}}$ (b) on the $(\unicode[STIX]{x1D713},k)$-plane for $\unicode[STIX]{x1D703}_{hd}=0.75$ and $p=-1$ in a pure heaving and deflection motion $(\mathsf{b}=\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{pd}=0)$.

Figure 9 shows the contours of the normalized thrust ${\hat{C}}_{Th}^{\ast }$ and efficiency $\hat{\unicode[STIX]{x1D702}}$ as $k$ and $\unicode[STIX]{x1D713}$ are varied with $\unicode[STIX]{x1D703}_{hd}=0.75$. It is observed that, for $1\lesssim k\lesssim 4$, the flexible foil generates more thrust than the rigid airfoil (${\hat{C}}_{Th}^{\ast }>1$) when $250^{\circ }\lesssim \unicode[STIX]{x1D713}\lesssim 100^{\circ }$, and less when $100^{\circ }\lesssim \unicode[STIX]{x1D713}\lesssim 250^{\circ }$. In fact, a phase shift $\unicode[STIX]{x1D713}$ between $100^{\circ }$ and $250^{\circ }$, approximately, corresponds to a motion of the foil with the trailing edge pointing downwards at the beginning of the downstroke, and opposite at the beginning of the upstroke (see figure 7b), generating less thrust than its rigid counterpart. In relation to the efficiency, figure 9(b) shows that the model yields a singularity for a particular combination of $k$ and $\unicode[STIX]{x1D713}$ due to the vanishing of the power coefficient in the present model, where $\unicode[STIX]{x1D702}$ changes from $+\infty$ to $-\infty$, so that it is expected that close to this curve on the $(\unicode[STIX]{x1D713},k)$-plane the propulsive efficiency reaches a local maximum, greatly enhancing the propulsive efficiency of a rigid heaving foil. For $k$ of order unity this corresponds, approximately, to $\unicode[STIX]{x1D713}$ between $250^{\circ }$ and $270^{\circ }$.

Figure 10. Curves ${\hat{C}}_{Th}^{\ast }=1$ (a), and both $\hat{\unicode[STIX]{x1D702}}=0$ and $\hat{\unicode[STIX]{x1D702}}=1$ (b), on the $(\unicode[STIX]{x1D713},k)$-plane for pure heaving and deflection motion $(\mathsf{b}=\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{pd}=0)$ with $0<\unicode[STIX]{x1D703}_{hd}\leqslant 1$ and $p=-1$.

To see the effect of $\unicode[STIX]{x1D703}_{hd}$ we plot the curves on the $(\unicode[STIX]{x1D713},k)$-plane corresponding to ${\hat{C}}_{Th}^{\ast }=1$ in figure 10(a), and both the curves corresponding to $\hat{\unicode[STIX]{x1D702}}=0$ and $\hat{\unicode[STIX]{x1D702}}=1$ (i.e. close to the singularity) in figure 10(b), for $0<\unicode[STIX]{x1D703}_{hd}\leqslant 1$. Figure 10(a) shows that the regions of enhanced and reduced thrust in relation to a heaving rigid foil remain almost independent of $\unicode[STIX]{x1D703}_{hd}$ for $k$ of order unity. For the efficiency (figure 10b) the situation is somewhat similar, but the singularity moves towards smaller values of $k$ as $\unicode[STIX]{x1D703}_{hd}$ increases. Thus, the present model predicts a maximum enhancement of the propulsive efficiency in relation to a rigid heaving foil when $k\approx 1$ for $\unicode[STIX]{x1D703}_{hd}$ approaching unity and $\unicode[STIX]{x1D713}\approx 260^{\circ }$. This result is in qualitative agreement with the experimental results of Ramananarivo et al. (Reference Ramananarivo, Godoy-Diana and Thiria2011) for a similar foil motion, where the maximum efficiency is found to be reached when the trailing-edge deflection angle is approximately equal to the effective angle of attack at the mid-downstroke, which is the situation for $\unicode[STIX]{x1D713}$ approaching $270^{\circ }$ in figure 7(b).

6.3 Pure pitching motion with chordwise deflection

As the third simple case, we consider a pure pitching motion combined with a deflection motion of the foil, i.e. the first oscillation mode of a flexible airfoil for pure pitching motion, for which (4.13) and (4.16) become

(6.5)$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{C}}_{Tp}=t_{p}(k,a)+t_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{pd}+t_{d}(k,p)\unicode[STIX]{x1D703}_{pd}^{2}, & \displaystyle\end{eqnarray}$$

(6.6)$$\begin{eqnarray}\displaystyle & \displaystyle {\hat{C}}_{Pp}=p_{p}(k,a)+p_{pd}(k,p,a,\unicode[STIX]{x1D713},\unicode[STIX]{x1D719})\unicode[STIX]{x1D703}_{pd}+p_{d}(k,p)\unicode[STIX]{x1D703}_{pd}^{2}, & \displaystyle\end{eqnarray}$$

with the ratio $\unicode[STIX]{x1D703}_{pd}$ between the deflection and pitching amplitudes defined in (4.11). Since the thrust of a rigid plate may change its sign for a pure pitching motion, to avoid singularities it is convenient to redefine the normalized thrust coefficient in this case as

(6.7)$$\begin{eqnarray}{\hat{C}}_{Tp}^{\ast }={\hat{C}}_{Tp}-t_{p}(k,a).\end{eqnarray}$$

Now ${\hat{C}}_{Tp}^{\ast }=0$ means the same thrust coefficient as that of the equivalent pitching rigid foil. The definition of the normalized efficiency (4.17) does not change. Both ${\hat{C}}_{Tp}^{\ast }$ and $\hat{\unicode[STIX]{x1D702}}$ depend on the non-dimensional parameters $k,a,p,\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D703}_{pd}$ (note that one may set $\unicode[STIX]{x1D719}=0$ since there is no heaving motion, and $\unicode[STIX]{x1D713}$ is now the phase shift of the deflection in relation to the pitch). Similarly to the previous case, to reduce the number of parameters we focus on a deflection pivoting at the leading edge ($p=-1$), which in the present case implies also a pitching motion about the leading edge, i.e. $a=p=-1$.

Figure 11. Contours of the normalized thrust ${\hat{C}}_{Tp}^{\ast }$ (a) and efficiency $\hat{\unicode[STIX]{x1D702}}$ (b) on the $(\unicode[STIX]{x1D713},k)$-plane for $\unicode[STIX]{x1D703}_{pd}=0.75$ and $a=p=-1$ in a pure pitching and deflection motion $(\mathsf{b}=\unicode[STIX]{x1D703}^{-1}=\unicode[STIX]{x1D703}_{hd}=0)$.

Figure 12. Curves ${\hat{C}}_{Tp}^{\ast }=0$ (a) and $\hat{\unicode[STIX]{x1D702}}=0$ (b) on the $(\unicode[STIX]{x1D713},k)$-plane for pure pitching and deflection motion $(\mathsf{b}=\unicode[STIX]{x1D703}^{-1}=\unicode[STIX]{x1D703}_{hd}=0)$ with $0<\unicode[STIX]{x1D703}_{pd}\leqslant 1$ and $a=p=-1$.

Figure 11 shows the contours of the normalized thrust ${\hat{C}}_{Tp}^{\ast }$ and efficiency $\hat{\unicode[STIX]{x1D702}}$ on the $(\unicode[STIX]{x1D713},k)$-plane for $\unicode[STIX]{x1D703}_{pd}=0.75$. For the thrust (figure 11a) we can distinguish two regions: one for small $k$ with two maxima of the relative thrust, around $\unicode[STIX]{x1D713}=0^{\circ }$ and for $\unicode[STIX]{x1D713}\gtrsim 270^{\circ }$, and another region for $k\gtrsim 1$ where the enhanced thrust is for $90^{\circ }\lesssim \unicode[STIX]{x1D713}\lesssim 270^{\circ }$. Figure 11(b) shows that these regions approximately coincide with those of enhanced efficiency.

The effect of $\unicode[STIX]{x1D703}_{pd}$ is plotted in figure 12 as curves corresponding to ${\hat{C}}_{Tp}^{\ast }=0$ and $\hat{\unicode[STIX]{x1D702}}=0$ on the $(\unicode[STIX]{x1D713},k)$-plane for $0<\unicode[STIX]{x1D703}_{pd}\leqslant 1$. The regions of enhanced and reduced thrust in relation to a pitching rigid foil are nearly independent of $\unicode[STIX]{x1D703}_{pd}$, with larger variations as $k$ increases. However, in contrast to the pure heaving motion, the regions of positive and negative relative efficiency $\hat{\unicode[STIX]{x1D702}}$ remain practically unchanged as $\unicode[STIX]{x1D703}_{pd}$ varies from zero to unity (figure 12b). Therefore, the above discussed region of maximum enhancement in the propulsive efficiency of a flexible foil in relation to a pitching rigid foil about its leading edge remains practically independent of the trailing-edge deflection amplitude in relation to the pitch amplitude, within the present linearized potential theory valid for small amplitudes.

7 Concluding remarks

The closed expressions obtained in the present work for the aerodynamic force components and moment on a two-dimensional flexible foil undergoing a quite general undulatory motion, which are validated against available numerical results for sufficiently high Reynolds numbers and small amplitude of the oscillations, constitute a convenient tool for predicting and evaluating the optimal conditions for propulsion in terms of thrust generation and efficiency in a wide range of animal and bioinspired robotic locomotion.

Relatively simple analytic expressions are obtained for the interesting cases of pitching and heaving motions superimposed to a chordwise flexibility of the foil, with the additional force and moment terms in relation to an oscillating rigid foil neatly separated, and characterized by the ratio between the amplitudes of the deflection and the rigid motions. When this non-dimensional parameter vanishes one recovers previous results for a pitching and heaving rigid foil. A detailed evaluation of the propulsion performance is made for these particular cases, mapping the regions of thrust and of propulsive efficiency enhancement in relation to the rigid foil counterpart in the parameter space of the reduced frequency and the relevant deflection parameters.

The present results are limited to small amplitudes of the oscillations and sufficiently high Reynolds number for which the linearized potential theory applies. In addition, no analysis is made of the fluid–structure interaction that may produce the particular deflection or undulatory motion of the foil. Previous works that have analysed this fluid–structure interaction have shown that, in some particular pitching and heaving motions of the foil, the regions of optimal propulsion are related to structural resonant frequencies of the foil (Michelin & Llewellyn Smith Reference Michelin and Llewellyn Smith2009; Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013; Moore Reference Moore2014; Paraz et al. Reference Paraz, Schouvelier and Eloy2016; Moore Reference Moore2017). The present analytical results are independent of the material properties of the foil, depending only on the prescribed kinematics of the flexible foil. Although to analyse the fluid–structure interaction one needs to model a truly general motion of the deformable foil, with infinite kinematic parameters, as recently done numerically in this limit of linear potential theory by Tzezana & Breuer (Reference Tzezana and Breuer2019), we believe that with the present fairly broad class of flapping undulatory motion, with nine non-dimensional parameters, one may undertake a similar, but more limited, analysis of the fluid–structure interaction to obtain approximately these kinematic parameters for given properties of the flexible or compliant material and the boundary conditions at the leading and trailing edges, with the advantage that now the aerodynamic force and moment are obtained analytically in terms of these parameters from a general impulse theory. This is left for future research.