2.1. Non-dimensional equations

The equations are non-dimensionalized using the fluid density $\rho$ , the free stream velocity $U$ and the half-chord length $c/2$ . To simplify the notation, the same symbols are kept for the coordinates $(x,z)$ , the time $t$ (scaled with $2U/c$ ) and the foil displacement $z_s$ . For $z_s(x,t)$ to take into account the pitch and heave motions of the foil about the leading edge, as well as two flexural deformation degrees of freedom, at least a fifth-order polynomial must be used in $x$ , instead of the quartic one used in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021). Thus, imposing a pitch and heave motion at the leading edge ( $x=-1$ ) and a free trailing edge ( $\partial ^2 z_s/\partial x^2 = \partial ^3 z_s/\partial x^3 =0$ at $x=1$ ), this fifth-order polynomial approximation can be written as

(2.8) \begin{align} z_s(x,t)&= h(t) - \alpha (t) (x+1) + d_1(t) [ 24(x+1)^2-8(x+1)^3+(x+1)^4] \nonumber \\ &\quad + d_2(t)[160(x+1)^2 -40(x+1)^3+(x+1)^5] \,, \quad \quad -1 \leqslant x \leqslant 1 , \end{align}

where $h(t)$ and $\alpha (t)$ characterize the heave and pitch motions, respectively, and the unknown functions $d_1(t)$ and $d_2(t)$ characterize the flexural deformation. Note that the pitch angle $\alpha$ is positive when it is clockwise, as it is usual in airfoil aerodynamics, and that ${d}(t)$ in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021) is $24 d_1(t)$ here. The corresponding non-dimensional velocity of the foil’s centreline can be written as

(2.9) \begin{align} v_0(x,t)&= \frac {\partial z_s}{\partial t}+ \frac {\partial z_s}{\partial x}\nonumber \\ &= \dot {h} - \alpha + (x+1) [- \dot {\alpha } + 48 \dot {d}_1+320 \dot {d}_2] + (x+1)^2 [24 \dot {d}_1 \nonumber \\ &\quad +160 \dot {d}_2-24 d_1-120 d_2]+(x+1)^3 [-8 \dot {d}_1-40 \dot {d}_2+4 d_1] \\ &\quad +(x+1)^4 [\dot {d}_1+5 d_2] + (x+1)^5 \dot {d}_2 , \nonumber \end{align}

where dots are used for the derivatives with respect to the non-dimensional time $t$ .

Assuming that $\rho _s$ , $\varepsilon$ and $E$ are constant along the foil’s chord length, the non-dimensional counterparts of (2.2)–(2.5) can be written as

(2.10) \begin{equation} R \left ( \ddot {h} - \ddot {\alpha } + \frac {96}{5} \ddot {d}_1 + \frac {416}{3} \ddot {d}_2 \right) =C_L+C_{L_p} , \end{equation}

(2.11) \begin{equation} R \left ( \frac {1}{2} \ddot {h} - \frac {2}{3} \ddot {\alpha } + \frac {208}{15} \ddot {d}_1 + \frac {704}{7} \ddot {d}_2 \right) =C_M+C_{M_p} , \end{equation}

(2.12) \begin{equation} R \left ( \frac {4}{3} \ddot {h} - 2 \ddot {\alpha } + \frac {4544}{105} \ddot {d}_1 + \frac {944}{3} \ddot {d}_2 \right) + S \left ( \frac {32}{3} d_1 + 80 d_2 \right) = C_{F_1} , \end{equation}

(2.13) \begin{equation} R \left ( 2 \ddot {h} - \frac {16}{5} \ddot {\alpha } + \frac {496}{7} \ddot {d}_1 + \frac {32 512}{63} \ddot {d}_2 \right) + S \left ( 16 d_1 + 128 d_2 \right) = C_{F_2} , \end{equation}

where the non-dimensional parameters

(2.14) \begin{equation} R = \frac {4 \rho _s \varepsilon }{\rho c} \quad \quad \text{and} \quad \quad S= \frac {4 E \epsilon ^3}{\rho U^2 c^3} \end{equation}

are the mass ratio (or inertia parameter) and the bending stiffness parameter, respectively, and the following fluid force and moment coefficients have been defined:

(2.15) \begin{equation} C_L = \frac {F_z}{\frac {1}{2} \rho U^2 c} = \int _{-1}^1 (\Delta P) {\textrm d} x \,, \quad \quad \text{with} \quad \quad \Delta P = \frac {\Delta p}{\rho U^2} , \end{equation}

(2.16) \begin{equation} C_M = \frac {M}{\frac {1}{2} \rho U^2 c^2} = \frac {1}{2} \int _{-1}^1 (x+1) (\Delta P) {\textrm d} x = \frac {1}{2} C_L + \frac {1}{2} \int _{-1}^1 x (\Delta P) {\textrm d} x , \end{equation}

(2.17) \begin{equation} C_{F_1} = \frac {D_1}{\rho U^2 \left ( \frac {c}{2} \right)^3} = \int _{-1}^1 (x+1)^2 (\Delta P) {\textrm d} x = - C_L + 4 C_M + \int _{-1}^1 x^2 (\Delta P) {\textrm d} x , \end{equation}

(2.18) \begin{equation} C_{F_2} = \frac {D_2}{\rho U^2 \left ( \frac {c}{2} \right)^4} = \int _{-1}^1 (x+1)^3 (\Delta P) {\textrm d} x = C_L - 6 C_M + 3 C_{F_1} + \int _{-1}^1 x^3 (\Delta P) {\textrm d} x \,; \end{equation}

together with the point force and moment coefficients,

(2.19) \begin{equation} C_{L_p} = \frac {F_{pz}}{\frac {1}{2} \rho U^2 c} , \quad C_{M_p} = \frac {M_p}{\frac {1}{2} \rho U^2 c^2} . \end{equation}

The above expressions of the fluid coefficients are written in terms of the non-dimensional pressure difference $\Delta P=(p^-p^+)/(\rho U^2)$ to facilitate their computation in terms of $h(t)$ , $\alpha (t)$ , $d_1(t)$ and $d_2(t)$ , and their derivatives (see § 3). Although the expressions for the coefficients will be obtained for any temporal variation of these functions, a general harmonic motion of the plate with frequency $\omega$ will be considered in the computations; i.e. the real parts (say) of

(2.20) \begin{equation} h(t)= h_0 e^{i k t} \,, \quad \alpha (t)= \alpha _0 e^{ikt} \, \quad {\textrm{with}} \quad \alpha _0= a_0 e^{i \phi } , \end{equation}

(2.21) \begin{equation} d_1(t)= d_{10} e^{i k t}, \quad d_2(t)= d_{20} e^{ikt}, \quad {\textrm{with}} \quad d_{10}= d_{1m} e^{i \psi _1}, \quad d_{20}= d_{2m} e^{i \psi _2} , \end{equation}

where $h_0$ and $a_0$ are the heave and pitch amplitudes, $d_{1m}$ and $d_{2m}$ the amplitudes of the first and second flexural deflection modes, and the heave phase is taken zero as usual, so that $\phi$ , $\psi _1$ and $\psi _2$ are the phase shifts of pitch and the two deflection modes, respectively, in relation to the heave phase. Finally,

(2.22) \begin{equation} k= \frac {\omega c}{2 U} = \frac {\pi f c}{U} \end{equation}

is the non-dimensional, or reduced, frequency.

2.2. Validation of the FSI model with the first two natural frequencies in vacuum

Before deriving the fluid coefficients $C_L, C_M, C_{F_1}$ and $C_{F_2}$ , it is shown here that the above formulation with the fifth-order polynomial approximation (2.8) for $z_s$ captures very accurately the first two natural frequencies of the foil in vacuum.

Substituting (2.20)–(2.21) into the moment equations (2.12) and (2.13) with $C_{F_1}=C_{F_2}=0$ , i.e. in absence of fluid–foil interaction, one obtains the following linear equation for the two deflection amplitudes $d_{10}$ and $d_{20}$ in terms of the pitch and heave amplitudes, $\alpha _0$ and $h_0$ , the structural parameters and the frequency $k$ :

(2.23) \begin{equation} { \mathsf{A}}_0 \cdot {\boldsymbol {d}}_0 = {\boldsymbol {b}}_0 , \end{equation}

with

(2.24) \begin{equation} { \mathsf{A}}_0 = \left ( \begin{array}{cc} \displaystyle \frac {32}{3} \, S - \displaystyle \frac {4544}{105} R k^2 & 80 \, S - \displaystyle \frac {944}{3} R k^2 \\ \\ 16 \, S - \displaystyle \frac {496}{7} R k^2 & 128 \, S - \displaystyle \frac {32512}{63} R k^2 \end{array} \right)\!, \end{equation}

(2.25) \begin{equation} {\boldsymbol {d}}_0 = \left ( \begin{array}{c} d_{10} \\ \\ d_{20} \end{array} \right)\!, \quad {\boldsymbol {b}}_0 = \left ( \begin{array}{c} \displaystyle \frac {4}{3} R k^2 h_0 - 2 R k^2 \alpha _0 \\ \\ 2 R k^2 h_0 - \displaystyle \frac {16}{5} R k^2 \alpha _{0} \end{array} \right) . \end{equation}

The resonant frequencies, or frequencies where the flexural deformation amplitudes $d_{10}$ and $d_{20}$ diverge, are obtained from det $({ \mathsf{A}}_0)=0$ , yielding the first two resonant frequencies in vacuum as the two real positive roots of that equation,

(2.26) \begin{equation} k_{r10}= \sqrt {\frac {7}{953} \left ( 629 -2 \sqrt {88\,189} \right) \frac {S}{R} } \simeq 0.507521 \sqrt {\frac {S}{R} } , \end{equation}

(2.27) \begin{equation} k_{r20}= \sqrt {\frac {7}{953} \left ( 629 +2 \sqrt {88\,189} \right) \frac {S}{R} } \simeq 2.997118 \sqrt {\frac {S}{R} } . \end{equation}

These values are in very good agreement with the exact results for the first two resonant frequencies of a beam clamped at the leading edge (e.g. Timoshenko, Young & Weaver Reference Timoshenko, Young and Weaver1974, § 5.11): $\omega _{i} = \lambda _i^2 \sqrt {E \varepsilon ^2 /(12 \rho _s \varepsilon c^4)}$ , with $\lambda _1=1.87510$ and $\lambda _2=4.69409$ , which in the present non-dimensional notation are $k_1=0.507491 \sqrt {S/R}$ and $k_2=3.180403 \sqrt {S/R}$ (relative errors with (2.26) and (2.27) of $0.006 \, \%$ and $5.7 \, \%$ , respectively). Thus, (2.8) is the simplest model for the chordwise deformation of a pitching and heaving flexible plate covering the first two resonant natural modes.

3.1. Lift coefficient

According to (2.15) and (3.2), the lift coefficient can be obtained from

(3.9) \begin{equation} C_L(t)= \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 (1-x) \varpi _s(x,t) {\textrm d}x + \int _{-1}^1 \varpi _s(x,t) {\textrm d} x . \end{equation}

After using Kelvin’s theorem (3.8), integrating by parts and performing the integral $\int _{-1}^1 x \varpi _{se} {\textrm d}x$ using (3.5), this can be written as (von Kármán & Sears Reference von Kármán and Sears1938)

(3.10) \begin{equation} C_L(t) = - \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x \varpi _0(x,t) {\textrm d}x - \int _{1}^\infty \frac {x \, \varpi _e(x,t)}{\sqrt {x^2-1}} {\textrm d} x \equiv C_{L_0}(t) + C_{L_e} (t). \end{equation}

The first term, $C_{L_0}$ , is the added-mass or non-circulatory contribution to the lift, which can be readily obtained once $\varpi _0(x,t)$ is computed from (3.6) using (2.9):

(3.11) \begin{align} \varpi _0(x,t) &= \sqrt {\frac {1 - x}{1 + x}} \, \left [ - 2 \dot {h} + \dot {\alpha } (4+2x)+2 \alpha + \dot {d}_1 \left ( - \frac {371}{4} -65 x - 5 x^2 + 6 x^3 - 2 x^4 \right) \right.\nonumber \\ &\quad + d_1\left ( -72 - 4 x + 16 x^2 - 8 x^3 \right) \nonumber \\ &\quad + \dot {d}_2 \left ( - \frac {1353}{2} - \frac {1943}{4} x - 46 x^2 + 49 x^3-12x^4-2 x^5 \right)\nonumber \\ &\quad \left. + \,d_2 \left ( - \frac {2175}{4} -45 x +135 x^2-50 x^3- 10 x^4 \right) \right] ; \end{align}

(3.12) \begin{align} C_{L_0}(t)& = - \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x \varpi _0(x,t) {\textrm d}x \nonumber \\ &= \pi \left ( \dot {\alpha } + \ddot {\alpha } - \ddot {h} - 25 \dot {d}_1 - \frac {149}{8} \ddot {d}_1 - \frac {1465}{8} \dot {d}_2 - \frac {1073}{8} \ddot {d}_2 \right) . \end{align}

Obviously, the contributions from the pitch and heave, $\alpha (t)$ and $h(t)$ , coincide with those obtained by Theodorsen (Reference Theodorsen1935) and von Kármán & Sears (Reference von Kármán and Sears1938), while the contribution from $d_1(t)$ coincide with that obtained in Alaminos-Quesada & Fernandez-Feria (Reference Alaminos-Quesada and Fernandez-Feria2020) with ${d}(t)=24 d_1(t)$ . Only the contributions from $d_2(t)$ , which captures the second natural frequency of the foil, is new here. Similarly it happens for all the other coefficients calculated in the rest of this section, except for $C_{F_2}$ , which is completely new, but their derivations are summarized here for completeness, and because the new contributions from $d_2(t)$ have to be obtained together with the rest to be integrated into the new expressions derived here.

To compute the circulatory contribution $C_{L_c}$ , the wake vorticity distribution $\varpi _e$ is obtained from Kelvin’s theorem (3.8) assuming the harmonic motion (2.20)–(2.21) and taking into account that $\varpi _e(x,t)=\varpi _e(X)$ , with $X=x-t$ , so that, writing $\Gamma _0(t)= G_0 e^{ikt}$ , it results that $\varpi _e(x,t)= g e^{i k (t-x)}$ with $g= -G_0 / \int _1^\infty \sqrt {(x+1)/(x-1)} \, e^{- i k x} {\textrm d}x$ . Solving this integral in terms of the Hankel functions $H_0^{(2)} (k)$ and $H_1^{(2)} (k)$ , and writing the result again in terms of $\Gamma _0(t)$ , one gets (von Kármán & Sears Reference von Kármán and Sears1938)

(3.13) \begin{equation} C_{L_c}(t) = - \int _{1}^\infty \frac {x}{\sqrt {x^2-1}} \, \varpi _e(x,t) {\textrm d} x = \Gamma _0(t) \mathcal{C} (k) , \end{equation}

where

(3.14) \begin{equation} \mathcal{C}(k)= \frac {H_1^{(2)} (k)}{H_1^{(2)} (k)+ i H_0^{(2)} (k)} \end{equation}

is Theodorsen’s function. This result coincides formally with Theodorsen (Reference Theodorsen1935) and von Kármán & Sears (Reference von Kármán and Sears1938), but now $\Gamma _0(t)$ contains more terms. It is obtained from (3.7) using (2.9):

(3.15) \begin{equation} \Gamma _0(t)= \pi \left [ - 2 \dot {h} + 3 \dot {\alpha } + 2 \alpha - \frac {263}{4} \dot {d}_1 - 59 d_1 - \frac {3831}{8} \dot {d}_2 - \frac {1755}{4} d_2 \right] . \end{equation}

Note that although (3.13) has been derived for a harmonic motion, it can be extrapolated to more general foil motions using this expression for $\Gamma _0(t)$ , not restricted to a harmonic motion of the foil.

3.2. Moment coefficient

Once $C_L$ is obtained, to compute $C_M$ from (2.16) it suffices to obtain

(3.16) \begin{equation} \int _{-1}^1 x (\Delta P) {\textrm d}x = \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 \frac {1-x^2}{2} \varpi _s(x,t) {\textrm d}x + \int _{-1}^1 x \varpi _s(x,t) {\textrm d} x . \end{equation}

Again, after some manipulations using Kelvin’s theorem (3.8), integration by parts, performing the integrals $\int _{-1}^1 x^2 \varpi _{se} {\textrm d}x$ and $\int _{-1}^1 x \varpi _{se} {\textrm d}x$ using (3.5) and taking into account that $d (\int _1^\infty A(x) \varpi _e(x-t) {\textrm d}x)/{\textrm d}t = \int _1^\infty A^{\prime}(x) \varpi _e(x-t) {\textrm d}x$ when $A(x=1)=0$ , one gets

(3.17) \begin{align} \int _{-1}^1 x (\Delta P) {\textrm d}x & = \frac {1}{4} \frac {{\textrm d} \Gamma _0}{{\textrm d}t} - \frac {1}{2} \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x^2 \varpi _0(x,t) {\textrm d}x\nonumber\\ &\quad + \int _{-1}^1 x \varpi _0(x,t) {\textrm d}x - \frac {1}{2} \int _1^\infty \frac {\varpi _e(x,t)}{\sqrt {x^2-1}} {\textrm d} x . \end{align}

Only the last term contributes to the circulatory moment. Combining this expression with $C_L$ in (2.16) and performing the integrals in a similar way to how it has been done for $C_L$ , the moment coefficient can be conveniently written as

(3.18) \begin{equation} C_M(t)=C_{M_0}(t)+ C_{M_c}(t) , \end{equation}

with the added-mass term

(3.19) \begin{align} C_{M_0}(t)&= \frac {\pi }{256} \left ( 192 \dot {\alpha } + 144 \ddot {\alpha } - 128 \ddot {h} - 576 d_1 -5248 \dot {d}_1 - 2800 \ddot {d}_1 \right.\nonumber \\ & \quad \left. - 4640 d_2 - 38680 \dot {d}_2-20213 \ddot {d}_2 \right)\!, \end{align}

and

(3.20) \begin{equation} C_{M_c}(t)= \tfrac {1}{4} \Gamma _0 (t) \mathcal{C}(k) , \end{equation}

where $\Gamma _0(t)$ is given by (3.15) and $\mathcal{C}(k)$ by (3.14).

3.3. First flexural moment coefficient

According to (2.17), to obtain $C_{F_1}$ one needs the additional computation of

(3.21) \begin{equation} \int _{-1}^1 x^2 (\Delta P) {\textrm d}x = \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 \frac {1-x^3}{3} \varpi _s(x,t) {\textrm d}x + \int _{-1}^1 x^2 \varpi _s(x,t) {\textrm d} x , \end{equation}

which after some lengthy manipulations can be written as

(3.22) \begin{align} \int _{-1}^1 x^2 (\Delta P) {\textrm d}x & = - \frac {1}{2} \Gamma _0 - \frac {1}{3} \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x^3 \varpi _0(x,t) {\textrm d}x \nonumber\\ & \quad + \int _{-1}^1 x^2 \varpi _0(x,t) {\textrm d}x - \frac {1}{2} \int _1^\infty \frac {x \varpi _e(x,t)}{\sqrt {x^2-1}} {\textrm d} x . \end{align}

Thus, performing the integrals and using the above expressions of $C_L$ and $C_M$ in (2.17), one gets

(3.23) \begin{equation} C_{F_1}(t)=C_{F_{10}}(t)+ C_{F_{1c}}(t) , \end{equation}

with

(3.24) \begin{align} C_{F_{10}}(t)&= \frac {\pi }{192} \left ( 384 \dot {\alpha } + 288 \ddot {\alpha } - 240 \ddot {h} - 1056 d_1 -10848 \dot {d}_1 - 5745 \ddot {d}_1 \right.\nonumber \\ & \quad \left. - 8640 d_2 - 80190 \dot {d}_2-41540 \ddot {d}_2 \right)\!, \end{align}

(3.25) \begin{equation} C_{F_{1c}}(t)= \tfrac {1}{2} \Gamma _0 (t) \mathcal{C}(k) . \end{equation}

Except for the new terms with $d_2(t)$ , this expression coincides with $C_F$ obtained in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021) by a slightly different approach, with $d(t)=24 d_1(t)$ in that reference.

3.4. Second flexural moment coefficient

Finally, to obtain $C_{F_2}$ one needs to compute

(3.26) \begin{equation} \int _{-1}^1 x^3 (\Delta P) {\textrm d}x = \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 \frac {1-x^4}{4} \varpi _s(x,t) {\textrm d}x + \int _{-1}^1 x^4 \varpi _s(x,t) {\textrm d} x , \end{equation}

which can be written as

(3.27) \begin{align} \int _{-1}^1 x^3 (\Delta P) {\textrm d}x & = \frac {3}{32} \frac {{\textrm d} \Gamma _0}{{\textrm d}t} - \frac {1}{4} \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x^4 \varpi _0(x,t) {\textrm d}x \nonumber\\ &\quad + \int _{-1}^1 x^3 \varpi _0(x,t) {\textrm d}x - \frac {3}{8} \int _1^\infty \frac {\varpi _e(x,t)}{\sqrt {x^2-1}} {\textrm d} x . \end{align}

Substituting into (2.18) and performing the integrals, the second flexural moment coefficient can be written as

(3.28) \begin{equation} C_{F_2}(t)=C_{F_{20}}(t)+ C_{F_{2c}}(t) , \end{equation}

with the added-mass term

(3.29) \begin{align} C_{F_{20}}(t)&= \frac {\pi }{128} \left ( 368 \dot {\alpha } + 280 \ddot {\alpha } - 224 \ddot {h} - 912 d_1 -10580 \dot {d}_1 - 5680 \ddot {d}_1 \right.\nonumber \\ & \quad \left. - 7530 d_2 - 78350 \dot {d}_2-41117 \ddot {d}_2 \right)\!, \end{align}

and the circulatory term

(3.30) \begin{equation} C_{F_{2c}}(t)= \frac {5}{8} \Gamma _0 (t) \mathcal{C}(k) , \end{equation}

where $\Gamma _0(t)$ is given by (3.15) and $\mathcal{C}(k)$ by (3.14).

4.1. Comparison with experimental data. Deflection correction with a nonlinear fluid damping

Although the natural frequencies around which the flexural deflection reaches its maximum are quite well predicted by the present formulation, it is well known that the amplitudes of these resonant deflections are greatly overestimated by the linear theory. Actually, they reach infinity in vacuum (see § 2.2), and, as shown in figures 2 and 3, $A_T/A_0$ may become much larger than unity near the natural frequencies, violating the linear hypothesis. This behaviour can be corrected by taking into account an external resistive term in the E–B equation modelling the lateral (or transverse) fluid dynamic drag (Taylor Reference Taylor1952; Eloy et al. Reference Eloy, Kofman and Schouveiler2012). Since the correction has to be mainly done around the natural frequencies, for outside these regions the flexural deflection amplitude remain within the limits of the linear theory, the use of a nonlinear damping term on the left-hand side of (2.1) of the form

(4.12) \begin{equation} \frac {1}{2} C_{Dz} \left | \frac {\partial z_s}{\partial t} \right | \frac {\partial z_s}{\partial t} , \end{equation}

where $C_{Dz}$ is an empirically fitted transverse drag coefficient, would be enough (Paraz et al. Reference Paraz, Schouveiler and Eloy2016; Piñeirua et al. Reference Piñeirua, Thiria and Godoy-Diana2017). In fact, only terms involving products of the flexural deflections $d_1$ and $d_2$ would be relevant in (4.12), because any other quadratic term would be negligible for small pitching and heaving amplitudes. But, for completeness, all the terms are retained.

In non-dimensional form, the nonlinear damping (4.12) would introduce, for the harmonic motion (2.20)–(2.21), the following terms into the left-hand side of the flexural deflection moment equations (2.12) and (2.13), respectively:

(4.13) \begin{align} & i k^2 C_{Dz} \left ( \left | \frac {4}{3} h_0 -2 \alpha _0 +k_{1h1} d_{10} + k_{1h2} d_{20} \right | h_0 + \left | -2 h_0 + \frac {16}{5} \alpha _0 +k_{1a1} d_{10} + k_{1a2} d_{20} \right | \alpha _0 \right .\nonumber \\ &\qquad \qquad + \left | k_{1h1} h_0+ k_{1a1} \alpha _0 + k_{111} d_{10}+k_{112} d_{20} \right | d_{10} \nonumber \\ &\qquad\! \qquad \left. + \,\left | k_{1h2} h_0+ k_{1a2} \alpha _0 +k_{112} d_{10}+k_{122} d_{20} \right | d_{20} \vphantom{\frac{1}{1}}\right)\!, \end{align}

and

(4.14) \begin{align} & i k^2 C_{Dz} \left ( \left | 2 h_0 - \frac {16}{5} \alpha _0 +k_{2h1} d_{10} + k_{2h2} d_{20} \right | h_0 \right. \nonumber\\ &\qquad \qquad+ \left | - \frac {16}{5} h_0 + \frac {16}{3} \alpha _0 +k_{2a1} d_{10} + k_{2a2} d_{20} \right | \alpha _0 \nonumber\\ & \qquad \qquad + \left | k_{2h1} h_0+ k_{2a1} \alpha _0 + k_{211} d_{10}+k_{212} d_{20} \right | d_{10} \nonumber\\ & \qquad \qquad \left. + \,\left | k_{2h2} h_0+ k_{2a2} \alpha _0 + k_{212} d_{10}+k_{222} d_{20} \right | d_{20} \vphantom{\frac{1}{1}}\right)\!, \end{align}

where the $k_{\ldots }\hbox{-}$ constants are given in Appendix B. Thus, instead of (4.1), the equation for determining the flexural deflection amplitude and phase becomes

(4.15) \begin{equation} \left ( { \mathsf{A}} + { \mathsf{A}}_c \right) \boldsymbol{\cdot} {\boldsymbol {d}}_0 = {\boldsymbol {b}} + {\boldsymbol {b}}_c , \end{equation}

where the new matrix ${ \mathsf{A}}_c$ and vector ${\boldsymbol {b}}_c$ , which depend on ${\boldsymbol {d}}_0$ , are also given in Appendix B.

This is a nonlinear equation for ${\boldsymbol {d}}_0$ which can easily be solved iteratively starting from the solution for $C_{Dz}=0$ given by (4.1). Since both ${ \mathsf{A}}_c$ and ${\boldsymbol {b}}_c$ are proportional to $k^2$ , the correction to the linear solution will be more important as the frequency $k$ increases, which is consistent with the fact that the largest deflection amplitude coming from the linear theory is commonly observed at the second natural frequency $k_{r2}$ , so that, as we shall see below, the nonlinear damping correction must be the greatest around this frequency.

Figure 4. Comparison of experimental results by Paraz, Eloy & Schouveiler (Reference Paraz, Eloy and Schouveiler2014) for $A_T/A_0=A_{TE}/A_{LE}$ versus the frequency (triangles, case $B=0.053$ N m in their figure 5a) with the present linear results ( $C_{Dz}=0$ , dashed line), and damped with $C_{Dz}=1, \,10$ and $12$ (solid lines). Also shown with a dotted line is the linear result from Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021), which only captures the first natural frequency.

Figure 4 shows a comparison of the flexural deflection amplitude obtained from the present theoretical formulation with experimental results reported by Paraz et al. (Reference Paraz, Eloy and Schouveiler2014) for a flexible plate with $c=12$ cm undergoing a pure heaving motion in a water current with velocity $U=0.05$ m s^–1 ( $Re=U c/\nu = 6000$ ). The mass ratio of the foil is $R=0.16$ and its bending rigidity $0.053$ N m, corresponding to the non-dimensional stiffness $S=589$ in the present notation. This example also serves to illustrate how the nonlinear fluid damping described above corrects the flexural deflection amplitude as the drag coefficient $C_{Dz}$ and the forcing frequency $f$ are varied. The figure plots the trailing-edge amplitude $A_{TE}$ normalized with the heaving amplitude, $A_{LE}=0.004$ m ( $h_0 \simeq 0.067$ ), a ratio that coincides with $A_T/A_0$ defined in (4.10) for a pure heaving motion, versus the dimensional frequency in hertz. The present linear theory with no damping ( $C_{Dz}=0$ , dashed line) presents two marked maxima at the first two natural frequencies, $f_{r1}=1.03$ Hz and $f_{r2}=7.90$ Hz. Also shown for comparison sake is the deflection obtained with a quartic polynomial approximation for $z_s$ in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021) (dotted line with $C_{Dz}=0$ ), which captures only the first resonant frequency, obviously coinciding with the present results around $f_{r1}$ , but decaying beyond that frequency without presenting any other peak. The trailing edge deflection ratio $A_T/A_0$ is also plotted for three values of $C_{Dz} \ne 0$ to show how the linear resonant peaks are damped as $C_{Dz}$ increases, but barely affecting the rest of the deflection amplitude. It is found that the present theoretical results with $C_{Dz}=10$ best fit the experimental results around the first natural frequency, which is close to the value $C_{Dz}=12$ found by Paraz et al. (Reference Paraz, Schouveiler and Eloy2016) to best fit their analytical model with the experimental data. The present results with $C_{Dz}=12$ are also plotted to show that their differences with $C_{Dz}=10$ are quite small. Around the second natural frequency, the agreement of the present results (using either $C_{Dz}=10$ or $12$ ) with the experimental data is poorer than around the first natural frequency, but they qualitatively capture the magnitude of the experimental deflection. Figure 5 illustrates the shape of the flexible plate along a flapping cycle at the two natural frequencies.

Figure 5. Plate shape at different times for the same case plotted in figure 4 ( $h_0=0.067$ , $\alpha _0=0$ , $R=0.16$ , $S=589$ ) at the first (a) and second (b) resonant frequencies.

Further comparison of the present flexural deflection theoretical results with other experimental data is provided in § 5.4.

5.1. Thrust coefficient

Before deriving the expression for the thrust force generated by the flexible foil’s kinematics, (2.8)–(2.9), in the present framework of linear potential theory, it is instructive to write the lift coefficient (2.15) in terms of the vorticity distribution, but not in the form (3.9) used in § 3.1, but according to the general vortical impulse theory (Wu Reference Wu1981; Saffman Reference Saffman1992).

Starting from (3.9) and using Kelvin’s theorem, $\int _{-1}^1 \varpi _s(x,t) {\textrm d}x + \int _1^\infty \varpi _e (x,t) {\textrm d}x =0$ , it can be written as

(5.1) \begin{equation} C_L(t)= - \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x \varpi _s(x,t) {\textrm d}x - \frac {\textrm d}{{\textrm d}t} \int _1^\infty \varpi _e(x,t) {\textrm d}x - \int _1^\infty \varpi _e(x,t) {\textrm d} x . \end{equation}

Now, since $\varpi _e(x,t)=\varpi _e(X)$ , with $X=x-t$ (von Kármán & Sears Reference von Kármán and Sears1938),

(5.2) \begin{align} \frac {\textrm d}{{\textrm d}t} \int _1^\infty (x-1) \varpi _e(X) {\textrm d}x = \frac {\textrm d}{{\textrm d}t} \int _{1-t}^\infty (X+t-1) \varpi _e(X) {\textrm d} X = \int _1^\infty \varpi _e(x,t) {\textrm d}x . \end{align}

Substituting this expression into the last term of (5.1) one gets the lift coefficient as the $z\hbox{-}$ component of the 2-D force $\boldsymbol {F}$ in terms of vortical impulse, ${\boldsymbol {F}} = - \rho d (\int _{\mathcal{V}_\infty } ({\boldsymbol {x}} \wedge {\boldsymbol{\omega }}) {\textrm d}\mathcal{V})/{\textrm d} t$ , where $\boldsymbol{\omega }$ is the vorticity field and $\mathcal{V}_\infty$ the entire fluid domain. In non-dimensional form and in the small-amplitude linearized limit,

(5.3) \begin{equation} C_L(t)= - \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 x \varpi _s(x,t) {\textrm d}x - \frac {\textrm d}{{\textrm d}t} \int _1^\infty x \varpi _e(x,t) {\textrm d}x . \end{equation}

This expression was already used by von Kármán & Sears (Reference von Kármán and Sears1938) in their pioneering work on unsteady airfoil theory, recovering the previous result by Theodorsen (Reference Theodorsen1935) from the standard expression (2.15) for a rigid foil undergoing a general harmonic pitching and heaving motion of small amplitude. Likewise, the non-dimensional $x\hbox{-}$ component of the force $C_x$ , or better, the thrust coefficient $C_T=-C_x$ , can be written as

(5.4) \begin{equation} C_T(t)= - \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 z_s(x,t) \varpi _s(x,t) {\textrm d}x - \frac {\textrm d}{{\textrm d}t} \int _1^\infty z_e(x,t) \varpi _e(x,t) {\textrm d}x , \end{equation}

where $z_e(x,t)$ is the displacement of the vortex wake in the $z\hbox{-}$ direction. This expression, valid for $|z_s| \ll 1$ and $|z_e| \ll 1$ as it is (5.3), is now used for the flexible foil’s kinematics (2.8)–(2.9).

Assuming that $z_e(x,t)$ is, like $\varpi _e(x,t)$ , a function only of $X=x-t$ , i.e. that the wake is stationary relative to the fluid (von Kármán & Sears Reference von Kármán and Sears1938), the last term can be written as

(5.5) \begin{equation} - \frac {\textrm d}{{\textrm d}t} \int _1^\infty z_e(x,t) \varpi _e(x,t) {\textrm d}x = - z_e(1,t) \varpi _e(1,t) = - z_s(1,t) \varpi _e(1,t) . \end{equation}

The thrust force obtained with this assumption in the simplest case of a rigid pitching and heaving foil has shown to compare very well with several sets of experimental data for small-amplitude oscillations (Fernandez-Feria Reference Fernandez-Feria2016, Reference Fernandez-Feria2017; Alaminos-Quesada Reference Alaminos-Quesada2021).

Taking into account (3.3), the first term in (5.4) contains terms associated with $\varpi _0(x,t)$ and to $\varpi _{se}(x,t)$ . The terms corresponding to $\varpi _0$ are easily computed using (2.8) and (3.11), to yield

(5.6) \begin{equation}\! - \frac {\textrm d}{{\textrm d}t}\! \int _{-1}^1 \! z_s(x,t) \varpi _0(x,t) {\textrm d}x = - \frac {\textrm d}{{\textrm d}t} \left [ h(t) \Gamma _0(t) + \alpha (t) \Gamma _a(t) + d_1(t) \Gamma _{d_1}(t) + d_2(t) \Gamma _{d_2}(t) \right]\!, \end{equation}

where $\Gamma _0(t)$ , defined in (3.7), is given by (3.15), whereas the other functions $\Gamma _{\ldots } (t)$ are given in Appendix C. It is understood that the terms inside brackets are the product of the real parts of each factor; e.g. $\textrm{Re} [h(t)] \textrm{Re} [\Gamma _0(t)]$ .

For the terms coming from $\varpi _{se}$ , one can first use Kelvin’s theorem (3.8), and then make use the expression (3.5) for $\varpi _{se}$ to integrate in $x$ (except for the term corresponding to $h(t)$ in $z_s(x,t)$ for which this last step is not needed):

(5.7) \begin{align} &- \frac {\textrm d}{{\textrm d}t} \int _{-1}^1 z_s(x,t) \varpi _{se}(x,t) {\textrm d}x = - \frac {\textrm d}{{\textrm d}t} \left [ h(t) \left ( - \Gamma _0(t) - \int _1^\infty \varpi _e(\xi,t) {\textrm d}\xi \right) \right. \nonumber \\ &\quad + \alpha (t) \left (\Gamma _0(t) + \int _1^\infty \varpi _e(\xi,t) {\textrm d}\xi - \int _1^\infty I_a(\xi) \varpi _e(\xi,t) {\textrm d}\xi \right) \nonumber \\ &\quad + d_1(t) \left ( - \frac {163}{8} \Gamma _0(t) - \frac {163}{8} \int _1^\infty \varpi _e(\xi,t) {\textrm d}\xi + \int _1^\infty I_{d1}(\xi) \varpi _e(\xi,t) {\textrm d}\xi \right) \nonumber \\ &\quad \left. +\, d_2(t) \left ( - \frac {1183}{8} \Gamma _0(t) - \frac {1183}{8} \int _1^\infty \varpi _e(\xi,t) {\textrm d}\xi + \int _1^\infty I_{d2}(\xi) \varpi _e(\xi,t) {\textrm d}\xi \right) \right]\!, \end{align}

where

(5.8) \begin{equation} I_a(\xi) = \sqrt {\xi ^2-1}-\xi , \end{equation}

(5.9) \begin{equation} I_{d1}(\xi) = \frac {1}{8} \left [27-224 \xi -48 \xi ^2+32 \xi ^3-8 \xi ^4+ \sqrt {\xi ^2-1} \left (8 \xi ^3-32 \xi ^2+52 \xi +208\right)\right]\!, \end{equation}

(5.10) \begin{align} I_{d2}(\xi) & = \frac {1}{8} \left [215-1640 \xi -400 \xi ^2+240 \xi ^3-40 \xi ^4-8 \xi ^5 \vphantom{\frac{1}{1}}\right.\nonumber \\ &\quad \left. +\, \sqrt {\xi ^2-1} \left (8 \xi ^4 + 40\xi ^3-236 \xi ^2+420 \xi +1523\right)\right] . \end{align}

The factors $163/8$ and $1183/8$ appear when applying Kelvin’s theorem to the terms corresponding to $d_1$ and $d_2$ , respectively, because these factors are subtracted to the functions of $x$ multiplying $d_1$ and $d_2$ in $z_s$ , so that the functions $I_{d1}(\xi)$ and $I_{d2}(\xi)$ appearing in the integrals in $\xi$ are regular at the leading edge $\xi =1$ ; particularly, their values at $\xi =1$ are

(5.11) \begin{equation} I_{d1}(1) = - \frac {221}{8} \,, \quad I_{d2}(1) = - \frac {1633}{8} . \end{equation}

This regularity is important to compute the time derivatives of these integrals taking into account that $\varpi _e(\xi,t) = \varpi _e(X)$ , with $X=\xi -t$ , so that, for any function $G(\xi)$ regular at $\xi =1$ , the time derivative can be written as (von Kármán & Sears Reference von Kármán and Sears1938)

(5.12) \begin{align} \frac {\textrm d}{{\textrm d}t} \int _1^\infty G(\xi) \varpi _e(\xi,t) {\textrm d}x = \int _1^\infty \frac {{\textrm d}G(\xi)}{{\textrm d}\xi } \varpi _e(\xi,t) {\textrm d}\xi + G(1) \varpi _e(1,t) . \end{align}

Performing these time derivatives in (5.7) and then the resulting integrals in $\xi$ , taking into account that $\varpi _e(\xi,t)$ can be written as (see text above (3.13))

(5.13) \begin{equation} \varpi _e(\xi,t)= \Gamma _0(t) \mathcal{H}(k) e^{-i k \xi } \,, \quad \mathcal{H} (k) = \frac {2}{\pi } \frac {1}{i H_0^{(2)} (k) + H_1^{(0)}(k) } , \end{equation}

and substituting (5.5), (5.6) and (5.7) into (5.4), the thrust coefficient can finally be expressed as

(5.14) \begin{equation} C_T(t)= C_{T0}(t) + C_{Tc}(t) , \end{equation}

with

(5.15) \begin{equation} C_{T0}= - \frac {\textrm d}{{\textrm d}t} \left [ \alpha \Gamma _a +d_1 \Gamma _{0d1} + d_2 \Gamma _{0d2} \right]\!, \end{equation}

and

(5.16) \begin{align} C_{Tc} &= \left ( \dot {h}- \dot {\alpha } + \frac {163}{8} \dot {d}_1+ \frac {1183}{8} \dot {d}_2 \right) [\Gamma _0 \mathcal{C}_h] + \dot {\alpha } [\Gamma _0 \mathcal{C}_{a1}] + \alpha [\Gamma _0 \mathcal{C}_{a0}] \nonumber \\ &\quad + \dot {d}_1 [\Gamma _0 \mathcal{C}_{d11}] + d_1 [\Gamma _0 \mathcal{C}_{d10}] + \dot {d}_2 [\Gamma _0 \mathcal{C}_{d21}] + d_2 [\Gamma _0 \mathcal{C}_{d20}] , \end{align}

where the different functions $\Gamma _{0 \ldots }(t)$ and complex functions $\mathcal{C}_{\ldots }(k)$ are given in the Appendices C and D, respectively. Note that all the terms involving $\varpi _e(1,t)$ in (5.7), some of them affected by the constants (5.11) after time derivation, cancel exactly out, this being an indication of the coherence of the method and the correctness of the computations. All the terms in (5.15) and (5.16) are actually the product of the real parts of two factors, the second one always written between brackets in (5.16) for clarity sake. The new contributions to $C_T$ due to the second flexural mode $d_2(t)$ are, in addition to those which are explicit in (5.15)–(5.16), scattered inside the different functions $\Gamma _{0\ldots }(t)$ (see Appendix C). With the use of these $\Gamma _{0\ldots }$ functions, the present notation (5.14)–(5.16) is more compact, and with a clearer separation between added-mas and circulatory terms, than in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021).

5.2. Cycle-averaged thrust

Time-averaged quantities over the flapping cycle $2 \pi /k$ are more useful for practical applications than the instantaneous ones. Thus, one may define the cycle-averaged (non-dimensional) thrust as

(5.17) \begin{equation} \overline {C}_T= \frac {k}{2 \pi } \int _t^{t+ 2 \pi /k} C_T(t) {\textrm d} t . \end{equation}

Clearly, the time average of the added-mass component (5.15) vanishes exactly because it is the time derivative of products of periodic functions of $t$ with zero mean. Therefore, all the time-averaged thrust comes from the circulatory part (5.16), and can be written as

(5.18) \begin{align} \overline {C}_T = \overline {C}_{Tc} &= t_h h_0^2+t_a a_0^2 + t_{d_1} d_{1m}^2 + t_{d_2} d_{2m}^2 + t_{ha} h_0 a_0 +t_{hd1} h_0 d_{1m} \nonumber \\ &\quad +t_{hd2} h_0 d_{2m} + t_{ad1} a_0 d_{1m} +t_{ad2} a_0 d_{2m} + t_{d1d2} d_{1m} d_{2m} , \end{align}

where the coefficients $t_h, t_a, t_{d_1}, t_{d_2}, t_{ha}, t_{hd1}, t_{hd2}, t_{ad1}, t_{ad2}$ and $t_{d1d2}$ multiplying the heave, pitch and flexural deformation amplitudes, which are functions of the reduced frequency $k$ , and eventually of the phases $\phi$ , $\psi _1$ and $\psi _2$ , are given in Appendix E. For $d_{2m}=0$ this expression coincides with that in Fernandez-Feria & Alaminos-Quesada (Reference Fernandez-Feria and Alaminos-Quesada2021) for a pivot axis at the leading edge, though written in a slightly different notation.

An offset drag $C_{Dx}$ is usually subtracted to this theoretical expression coming from the linear potential theory to account for the streamwise viscous friction. Its value is obviously much smaller than the transversal drag coefficient $C_{Dz}$ used above to account for the nonlinear viscous damping of the flexural deflection amplitude near the resonant frequencies. Typically, $C_{Dx} \approx 0.05$ for small amplitude pitching and heaving motions (Mackowski & Williamson Reference Mackowski and Williamson2015; Senturk & Smits Reference Senturk and Smits2019). Its correction to $\overline {C}_T$ is usually quite irrelevant, except for very small frequencies, since the theoretical $\overline {C}_T \to 0$ as $k \to 0$ . However small, its use important to avoid an unphysical singularity in the efficiency as $k \to 0$ (e.g. Mackowski & Williamson Reference Mackowski and Williamson2015; Fernandez-Feria Reference Fernandez-Feria2017).

5.3. Power coefficient and efficiency

The power required to generate the pitching and heaving motion about the leading edge is given, in non-dimensional form, by

(5.19) \begin{equation} C_P(t)= \dot {h}(t) C_{L_p}(t) - 2 \dot {\alpha }(t) C_{M_p}(t) . \end{equation}

On using (2.10)–(2.11) and the results of §§ 3.1 and 3.2, it can be written as

(5.20) \begin{equation} C_P= C_{P_R} + C_{P_F} = C_{P_R}+ C_{P_0}+C_{P_c} , \end{equation}

where

(5.21) \begin{equation} C_{P_R} = R \left [ \dot {h} \left ( \ddot {h} - \ddot {\alpha } + \frac {96}{5} \ddot {d}_1 + \frac {416}{3} \ddot {d}_2 \right) - 2 \dot {\alpha } \left ( \frac {1}{2} \ddot {h} - \frac {2}{3} \ddot {\alpha } + \frac {208}{15} \ddot {d}_1 + \frac {704}{7} \ddot {d}_2 \right) \right] \end{equation}

is the contribution due to the inertia of the foil of mass ratio $R$ , whereas

(5.22) \begin{equation} C_{P_F} = - \dot {h} C_L + 2 \dot {\alpha } C_M \end{equation}

is the fluid contribution, which in turn can be decomposed into an added-mass part, associated with the fluid inertia, and circulatory part,

(5.23) \begin{equation} C_{P_0} = - \dot {h} C_{L_0} + 2 \dot {\alpha } C_{M_0} , \end{equation}

(5.24) \begin{equation} C_{P_c} = - \dot {h} C_{L_c} + 2 \dot {\alpha } C_{M_c} = \left ( - \dot {h} + \frac {1}{2} \dot {\alpha } \right) [\Gamma _0 \mathcal{C}(k)] , \end{equation}

where $C_{L_0}$ and $C_{M_0}$ are given by (3.12) and (3.19), respectively, and use has been made of (3.13) and (3.20) for $C_{L_c}$ and $C_{M_c}$ . The contribution to $C_P$ associated with the nonlinear damping term (4.12) in (2.10)–(2.11) has been neglected because it has already used to damp the resonant peaks in the flexural deformation amplitude through its inclusion in (2.12) and (2.13), so that the resulting cubic terms from this nonlinear damping in $\dot {h}C_{L_p}$ and $\dot {\alpha } C_{M_p}$ become negligible once $|d_1|$ and $|d_2|$ are small.

The time-averaged power coefficient,

(5.25) \begin{equation} \overline {C}_P = \frac {k}{2 \pi } \int _t^{t+ 2 \pi /k} C_P(t) {\textrm d} t , \end{equation}

can also be decomposed into the above three contributions. The contribution from the inertia of the foil vanishes for a rigid foil (Fernandez-Feria Reference Fernandez-Feria2023), i.e. when the flexural deflection amplitudes $d_{1m}$ and $d_{2m}$ are zero:

(5.26) \begin{align} \overline {C}_{P_R} &= - \frac {16}{105} R k^3 \left \{ d_{1m} \left [ 91 a_0 \sin (\phi -\psi _1) + 63 h_0 \sin \psi _1 \right] \right. \nonumber \\ &\quad\left. +\, d_{2m} \left [ 660 a_0 \sin (\phi -\psi _2) + 455 h_0 \sin \psi _2 \right] \right \} . \end{align}

This contribution is important at high frequencies when $R$ is not small, e.g. for a flapping foil in air. As it will seen in the results reported below, it may help to reduce the power consumption and thus to increase the propulsion efficiency for large $R$ .

In general, the total time-averaged power coefficient can be written as

(5.27) \begin{align} \overline {C}_{P}= p_h h_0^2+p_a a_0^2 + p_{ha} h_0 a_0 +p_{hd1} h_0 d_{1m} +p_{hd2} h_0 d_{2m} + p_{ad1} a_0 d_{1m} +p_{ad2} a_0 d_{2m} , \\[-15pt] \nonumber\end{align}

where the coefficients $p_h, p_a, p_{ha}, p_{hd1}, p_{hd2}, p_{ad1}$ and $p_{ad2}$ , which are functions of the reduced frequency $k$ , and eventually of the mass ratio $R$ and the phases $\phi$ , $\psi _1$ and $\psi _2$ , are given in Appendix F. These coefficients are written there separating the foil inertia contribution (5.26) (terms containing the parameter $R$ ), the circulatory contribution (terms containing the real, $\mathcal{F}$ , or imaginary, $\mathcal{G}$ , parts of Theodorsen’s function $\mathcal{C}$ ), and the added-mass contribution (remaining terms). It is worth noticing from the expressions given there that $\overline {C}_{P_0}$ also contains cubic terms in $k$ , like $\mathcal{C}_{P_R}$ , but also quadratic and linear terms, whereas the terms in $\overline {C}_{P_c}$ are, at most, quadratic in $k$ . This means that what could most affect power, and therefore, efficiency near the second natural frequency are the terms coming from the inertia, both of the fluid and of the solid.

Once the thrust and power are computed, the (Froude) propulsive efficiency can be obtained as $\eta = \overline {C}_T/\overline {C}_P$ . However, this expression is valid provided that $\overline {C}_P$ is positive, and this is not always so because some contributions to $C_P$ may become negative, especially those coming from the inertia (e.g. Yin & Luo Reference Yin and Luo2010; Fernandez-Feria Reference Fernandez-Feria2023), so that there is energy transfer from the foil back to the driving system. To avoid this situation, a conservative assumption is usually adopted (Yin & Luo Reference Yin and Luo2010), which the negative power in either one of the two terms in (5.19) is not reusable. Consequently, we define an alternative time-averaged power coefficient $\overline {C}_{P^+}$ where only the positive parts of $\dot {h}(t) C_{L_p}(t)$ and $- 2 \dot {\alpha }(t) C_{M_p}(t)$ in (5.19) are averaged over the entire flapping cycle. The propulsion efficiency is thus defined as

(5.28) \begin{equation} \eta = \frac {\overline {C}_T}{\overline {C}_{P+}} . \end{equation}

In what follows, $\overline {C}_{P}$ will be used to simplify the notation, but understanding that it means $\overline {C}_{P^+}$ , except otherwise specified.

5.4. Comparison with experimental results

Figure 6 compares the present results for the trailing-edge deflection amplitude, the thrust coefficient, the power coefficient and the efficiency with experimental results by Quinn et al. (Reference Quinn, Lauder and Smits2014) for a heaving foil with $h_0 \simeq 0.1$ in a water current with $U \simeq 0.1$ m s^–1. The experimental data, extracted from figure 9 of Paraz et al. (Reference Paraz, Schouveiler and Eloy2016) and plotted against $f/f_{r1}=k/k_{r1}$ , are for the panel A in Quinn et al. (Reference Quinn, Lauder and Smits2014), with bending stiffness per unit span $E \varepsilon ^3/12 \simeq 1.96$ N m and thickness-to-chord ratio $\varepsilon /c =0.0034$ , which correspond to $S \simeq 1269$ and $R\simeq 0.1$ in the present notation. The flexural deflection plotted in figure 6(a) is obtained with $C_{Dz}=12$ , a value also used by Paraz et al. (Reference Paraz, Schouveiler and Eloy2016) in their model, showing a very good agreement with the experimental data of Quinn et al. (Reference Quinn, Lauder and Smits2014). An offset drag coefficient $C_{Dx}=0.05$ is used in the theoretical $\overline {C}_T$ plotted in figure 6(b), in accordance with the experimental data reported in figure 13 of Quinn et al. (Reference Quinn, Lauder and Smits2014) as the frequency tends to zero. The corresponding time-averaged power coefficient $\overline {C}_P$ and efficiency $\eta$ are plotted in figures 6(c) and 6(d), respectively, showing also a good agreement with the experimental data. For reference sake, the theoretical results for the rigid foil counterpart are also plotted with dashed lines, showing the enhancement of the propulsive performance due to flexibility, especially around the two natural frequencies ( $k_{r1} \simeq 11.6$ and $k_{r2}\simeq 89.6$ in this case).

Figure 6. Comparison between the present theoretical results (continuous lines) for $A_T/A_0$ (a), $\overline {C}_T/(k^2 h_0^2)$ (b), $\overline {C}_P/(k^3 h_0^2)$ (c) and $\eta$ (d) with experimental data (symbols) by Quinn, Lauder & Smits (Reference Quinn, Lauder and Smits2014) for a flexible heaving foil with $h_0 \simeq 0.1$ . The experimental data are extracted from figure 9 of Paraz et al. (Reference Paraz, Schouveiler and Eloy2016). Dashed lines correspond to the results for the rigid foil counterpart.

Figure 7. Comparison between the present theoretical results for $\overline {C}_T/(k^2 h_0^2)$ at two high values of $S$ (continuous lines) with experimental data (symbols) by Paraz et al. (Reference Paraz, Schouveiler and Eloy2016) for a flexible heaving foil with $h_0 = 0.07$ .

The experimental results plotted in figure 6 do not reach frequencies around the second natural mode. Paraz et al. (Reference Paraz, Schouveiler and Eloy2016) report some experimental measurements of the thrust force for higher frequencies in their figure 8(a), but they are made ‘in water at rest’, which represents a difficulty for comparison with current theoretical results because $U$ is used here to make dimensionless $\overline {C}_T$ , $k$ and $S$ . Formally, their experimental data corresponds to $S \to \infty$ and, since these authors report results for $\overline {C}_T/( k^2 h_0^2)$ in the present notation, in whose non-dimensionalization does not intervene $U$ , figure 6 compares the experimental data for the case $h_0=0.07$ with the present theoretical results for two large values of the stiffness, $S=5 \times 10^3$ and $S=10^4$ . The mass ratio used is $R=0.16$ , that coincides with that in figure 4, and the same values of the drag coefficients $C_{Dz}$ and $C_{Dx}$ of the previous comparison are used. The present results underestimate the thrust around the first natural frequency and overestimate it around the second. But taking into account that the limit $U=0$ cannot be well captured by the present theoretical approximation, the agreement is not that bad, reproducing correctly the two thrust peaks near the first two natural frequencies. It should also be noted that the thrust is scaled with $k^2$ in figure 7, so that the greatest thrust peak really happens around the second mode.

All the other experimental data found by the author in the literature for frequencies near the second natural mode are for large oscillation amplitudes, or for a foil with small aspect ratio, or both, so that the present theory could not be applied.

Figure 8. Normalized flexural deflection amplitude at the trailing edge (a), thrust (b), power (c) and efficiency (d) for pure heave as $k$ and $S$ are varied when $R=0.1$ . Here $C_{Dz}=12$ and $C_{Dx}=0.05$ . The vertical dashed line marks the case depicted in figure 6. For the other lines related to the natural frequencies see caption of figure 2.