Analytical results for the propulsion performance of a flexible foil with prescribed pitching and heaving motions and passive small deflection

Abstract The propulsive performance of a flexible foil with prescribed pitching and heaving motions about any pivot point location and passive chordwise flexural deflection is analysed within the framework of the linear potential flow theory and the Euler–Bernoulli beam equation using a quartic approximation for the deflection. The amplitude of the flexural component of the deflection and its phase, the thrust force, input power and propulsive efficiency are computed analytically in terms of the stiffness and mass ratio of the plate, frequency, pivot point location and remaining kinematic parameters. It is found that the maximum flexural deflection amplitude, thrust and input power are related to the first fluid–structure natural frequency of the system, corresponding to the deflection approximation considered. The same relation is observed for the propulsive efficiency when an offset drag is included in the analytical expressions. These results, which are valid for small amplitude and sufficiently large stiffness of the foil, are compared favourably with previous related results when the foil pivots about the leading edge. The configurations generating maximum thrust and efficiency enhancement by flexibility are analysed in relation to those of an otherwise identical rigid foil.


Introduction
Flexibility is a relevant feature of the flapping appendages of natural fliers and swimmers, which can lead to greatly improved propulsive performance. This enhancement of the propulsive capabilities of flexible flapping foils has been analysed and explained in a † Email address for correspondence: ramon.fernandez@uma.es large number of works, especially for the simplest configuration of a two-dimensional flapping foil with chordwise flexibility, as in the pioneering experimental works of Heathcote, Martin & Gursul (2004) and Heathcote & Gursul (2007) (see, for example, Wu (2011), Wang & Zhang (2016) and Smits (2019) for comprehensive reviews). To the large number of dimensionless parameters characterizing the kinematics of flapping rigid foils, flexibility, and in general the fluid-structure interaction, adds some new ones such as the dimensionless stiffness S, or ratio of elastic to fluid forces, and the mass ratio R, or ratio of solid to fluid inertia, in addition to a number of geometric parameters characterizing the deformation of the foil. Therefore, the analysis of the propulsive performance of flexible foils becomes cumbersome, even for two-dimensional foils and once the numerical or experimental difficulties associated with solving or measuring the fluid-structure interaction have been circumvented, since many simulations or measurements have to be made to cover the influence of all the parameters on the propulsive performance (see, for example, Quinn, Lauder & Smits (2015) and Olivier & Dumas (2016) for representative experimental and numerical studies, respectively). Thus, to dispose of analytical solutions characterizing the fluid-structure interaction of flexible flapping foils, even for very simplified configurations, is very helpful, both to understand the propulsive performance of flying and swimming animals and to design and improve bio-inspired flying and swimming robotic models. These analytical approximations have to be searched within the framework of the two-dimensional linearized inviscid flow theory for small deformations of the foil, pioneered by Wu (1961), who considered a flapping plate that incorporated flexibility, allowing it to deform according to the fluid and elastic forces it experiences. Passive flexibility changes the thrust that a flapping plate generates, and consequently its propulsive efficiency. Particularly, compared with rigid foils, it has generally been found from these inviscid flow theories that flexibility produces greater thrust when actuated near a fluid-structure natural frequency, but less otherwise, with propulsive efficiency greater than that of a rigid foil over a broad range of frequencies and stiffnesses (Katz & Weihs 1978;Alben 2008;Michelin & Llewellyn Smith 2009;Alben et al. 2012;Floryan & Rowley 2018).
An interesting analytical approach, within the linear inviscid limit considered in the present work, was recently presented by Moore (2014), who considered the particular case of a rigid foil with a prescribed heaving displacement and a passive pitching motion about its leading edge, relating the finite torsional stiffness of the leading edge with the kinematics of the rigid pitching motion and, hence, with the propulsive performance of the foil. More general approaches based on Wu's inviscid theory (Wu 1961) and the Euler-Bernoulli beam equation for small amplitudes have appeared recently, where the displacement of the foil and the fluid motion are decomposed into a Chebyshev series with a large number (infinite in theory) of parameters that have to be obtained numerically (Alben 2008;Moore 2017;Floryan & Rowley 2018;Tzezana & Breuer 2019;Floryan & Rowley 2020). All these works considered heaving and/or pitching motions about the leading edge of the foil and computed the thrust force from the projection of the pressure force into the flight direction plus the leading-edge suction (Garrick 1936;Wu 1961Wu , 1971. Here, however, we follow a simpler approach that provides analytical results by just considering a quartic deflection of small amplitude superimposed on the heaving and pitching motions of the foil, now about any pivot point location along the foil, not just the leading edge. Technically, the approach may be considered as a low-order expansion, or minimal model of a flexible foil. For this kinematics with low-amplitude flexural motion, which structurally is valid when the stiffness of the foil is sufficiently large, we use analytical results derived very recently for the propulsive performance of a pitching and heaving foil with a small flexural deflection (Alaminos-Quesada & Fernandez-Feria 2020), computing the thrust force from a general vortical impulse formulation (Fernandez-Feria 2016). These kinematic-based results are coupled with the dynamics of the foil through the Euler-Bernoulli beam equation to compute the amplitude d m and the phase shift ψ of the flexural deflection in terms of R, S and the kinematic parameters of the prescribed heaving and pitching motions about any given pivot point. The thrust, lift and moment exerted by the fluid on the foil are then obtained using the analytic general results for prescribed kinematics and deformation based on the impulse formulation (Alaminos-Quesada & Fernandez-Feria 2020). Thus, the present paper puts in a proper dynamical context the results in Alaminos-Quesada & Fernandez-Feria (2020), which considered just the kinematic problem, i.e. these former results provided the propulsion performance of a flapping foil when not only the heaving and pitching components of the motion were prescribed (forced), but also its flexural deformation as a quadratic deflection. Here this flexural deflection, now as a quartic polynomial approximation, is obtained from the dynamics of the foil in terms of its structural properties (stiffness, density, thickness, etc.) through the beam bending equation, for any given heaving and pitching motion applied at an arbitrary pivot point. Thus, the effect of flexibility on the propulsion performance of the foil is finally obtained in terms of the dimensionless stiffness and mass ratio of the foil, in addition to all the kinematic parameters.
In the approximation presented in the present work, the input power is obtained directly from the dynamic beam equation in terms of the input force and torque exerted externally on the pivot point, plus the effect of the inertia of the foil through the mass ratio R. The results are valid for a wide range of mass ratios, covering the typical values of both flapping swimming and flapping flight with flexible foils.
The paper is organized as follows. The formulation of the problem, including the general moments of the Euler-Bernoulli equation relating the input force and torque with the lift, moment and flexural coefficients, is given in § 2 together with the kinematics of the foil and the corresponding analytical expressions for the above-mentioned coefficients. Section 3 presents the analytical expressions for the first natural frequency of the system, as well as the amplitude and phase shift of the passive flexural deformation, and the corresponding thrust coefficient, input power coefficient and propulsive efficiency. These expressions are used in § 4 to characterize the propulsive performance of heaving and pitching foils about any pivot point in terms of the forcing frequency, stiffness and mass ratio of the foil, with more detail for a purely heaving motion forced on any pivot point, not just the leading edge as in most previous works. Finally, the conclusions are presented in § 5.

Formulation of the problem
We consider a two-dimensional foil of chord length c immersed in an incompressible and nearly inviscid flow with constant free-stream speed U along the x axis. We use non-dimensional variables scaled with the half-chord length c/2 and the velocity U. A harmonic heaving motion along the z axis, h(t), is imposed on a given point x = a of the foil, superimposed to a pitching motion with angle α(t) around this pivot point (see figure 1). These motions are generated by a vertical force and a torque actuating locally at the pivot point x = a, which in non-dimensional form are C Li and C Mi , respectively, defined below and in principle are unknowns.
The amplitudes of the heaving and pitching motions are assumed small compared with the chord length c, i.e. |h| 1 and |α| 1, so that the foil, and every point of the trail of vortices that it leaves behind, may be considered to be on the plane approximation, with the plate extending from x = −1 to x = 1. In addition, we consider the flexibility of the foil assuming a large stiffness, so that its deflection in relation to a rigid foil is also very small.
where z s (x, t) is the displacement in the z direction of the centreline of the foil (see figure 1) and C p (x, t) is the pressure difference between the lower and upper sides of the foil, p = p − − p + , scaled with ρU 2 , where ρ is the fluid density. To this distribution of fluid pressure force, we have added on the right-hand side of (2.1) the aforementioned external punctual forcing C Li (t) and C Mi (t), which are the non-dimensional force and torque per unit span, scaled with ρU 2 c/2 and ρU 2 c 2 /2, respectively, actuating locally at x = a to generate the prescribed heaving and pitching motions, respectively, where δ(x − a) is Dirac's delta function centred at x = a and δ its derivative. (See appendix A for more details.) Note that the moment is positive when counterclockwise. Finally, the non-dimensional quantities are the mass ratio (ratio of solid to fluid inertia) and dimensionless stiffness (ratio of elastic to fluid forces) of the foil (e.g. Moore 2017), respectively, where ρ s is the foil's density and ε and E its thickness and elastic modulus, respectively. These quantities are allowed to vary chordwise, though we consider here only the case in which R and S are constants. Equation (2.1) is valid when ε/c 1, |z s | 1 and |∂z s /∂x| 1. It is convenient to derive equations for the lift and the moment coefficients from the fluid-foil interaction, defined as where the stiffness terms have been integrated by parts (note that no simplifying assumptions about S(x) or R(x) have been made so far). The integral of (2.1) multiplied by (x − a) 2 will also be of interest, which will be related to the flexural deflection of the foil: (2.6) Here C F may be termed as the flexural coefficient, since it is to the flexural motion what C M is to the pitching motion, or C L to the heaving motion of the foil. It must be noted that (2.4)-(2.6) are valid for −1 < a < 1 since the point about which the delta function is centred must be inside of the domain of integration. Thus, results given below where a = −1 or a = 1 are in fact for a = −1 + or a = 1 − , respectively, with 1. Equation (2.1), but without the terms with C Li and C Mi , has been recently solved numerically in the present limit of linear potential flow for a general foil's kinematics z s (x, t) and constant stiffness S by Moore (2017) and by Floryan & Rowley (2018). Here we follow a different approach by using the three moment equations (2.4)-(2.6) to obtain the passive flexural motion of the foil at its lowest order for given heaving and pitching motions about x = a. In a previous work (Alaminos-Quesada & Fernandez-Feria 2020) we used a quadratic approximation for z s such as to compute the coefficients C L , C M and C F in terms of the heaving, h(t), pitching, α(t), and flexural quadratic, d(t), motions. In this foil's kinematics, |d|, like |h| and |α|, are assumed small, which would be a valid approximation for sufficiently large stiffness S of the foil. However, for constant S, the stiffness term in (2.1) vanishes for such a quadratic deflection, and, consequently, all the corresponding terms in the moment equations (2.4)-(2.6), so that these equations cannot relate the quadratic deflection d with the stiffness S, if S is constant.
To circumvent this difficulty we assume a quartic approximation for the deflection to account for the effect of S on d(t) at the lowest order. As we shall see, this approximation constitutes a minimal model of the flexible foil that yields quite accurately the first natural frequency of the system, with the propulsion results agreeing quite well with previous ones for particular values of R and a (Floryan & Rowley 2018) (namely R = 0.01 and a = −1).
In this approximation of z s as a quartic polynomial, the cubic and the quartic terms are related to the quadratic one by the boundary conditions of a free trailing edge, namely ∂ 2 z s /∂x 2 = ∂ 3 z s /∂x 3 = 0 at x = 1, in addition to the conditions z s = h and ∂z s /∂x = −α at x = a already satisfied by (2.7). This yields (2.8) With this approximation for z s (x, t) and constant S and R, (2.4)-(2.6) become where a dot denotes derivative with respect to t and the three functions of a multiplyingd are (2.14) In what follows we assume a harmonic motion of the foil: is the reduced frequency and Re means real part. We assume that h 0 is real and with φ the phase shift between the heaving and pitching motions of the foil, ψ the phase shift between the heaving and flexural motions, a 0 the maximum pitch amplitude of the plate and d m the amplitude of the flexural component of the deflection. We work with the complex expressions knowing that we have to take the real part of the results. As already mentioned, for the harmonic motion of the foil with quadratic deflection, i.e. when (2.7) with (2.15a-c) are used, the coefficients C L , C M and C F have been obtained analytically in terms of h, α and d ((4.1), (4.6) and (4.8), respectively, in Alaminos-Quesada & Fernandez-Feria (2020), with p = a and δ = d in the present notation, and now C M is assumed positive when counterclockwise). The addition of cubic and quartic terms in (2.8) obviously modifies these expressions, which are obtained here following the same procedure developed in Alaminos-Quesada & Fernandez-Feria (2020). By writing the foil's deflection (2.8) as and the corresponding foil's velocity in the z direction as these coefficients can be written as is the quasi-steady circulation and is Theodorsen's function (Theodorsen 1935), with H (2) n (z) = J n (z) − iY n (z), n = 0, 1, being Hankel's function of the second kind and order n, related to the Bessel functions of the first and second kind J n (z) and Y n (z) (Olver et al. 2010). Remember that one has to take the real part of these expressions.

Analytic expressions for the deflection and the propulsive performance
To analyse the propulsion performance of a heaving and pitching plate, with given functions h(t) and α(t), when a passive flexural motion with small amplitude d(t) is allowed, one has just to solve (2.11) for d with C F given by (2.26), since (2.9) and (2.10) would provide the input force and torque, C Li and C Mi , needed to generate such heaving and pitching motions, and therefore to compute the input power coefficient, as described below.
Equation (2.11) with (2.26) can formally be written as where the complex functions F 1 and F 2 are (3.4) The flexural deflection amplitude d m obviously vanishes when the stiffness S → ∞, since F 2 is an affine function of S. It is interesting to note that the algebraic function F 2 provides the first natural frequency of the system, which is the only one recovered by the present approximation. In fact, when one neglects the fluid-structure interaction, i.e. when C F = 0, equation F 2 = 0 provides the first resonant frequency of a plate in vacuo pivoting at x = a in the present approximation: For a plate pivoting at the leading edge (a = −1), this expression yields k r0 0.497 √ S/R, in remarkable agreement with the exact result for the first resonant frequency of a beam clamped at x = −1, ω 1 = λ 2 1 EI/(ρ s εbc 4 ), with λ 1 = 1.87510, where I = εb/12 is the area moment of inertia and b the span (e.g. Timoshenko, Young & Weaver 1974), which in the present (dimensionless) notation can be written as k 1 = 0.50750 √ S/R. When the fluid-structure interaction is taken into account (i.e. considering the effect of C F ), F 2 is complex and does not vanish, but the frequency that minimizes |F 2 | provides the present approximation for the first natural frequency of the system, k r , for a given set of parameters. Figure 2 compares this frequency k r thus obtained with the first natural frequency computed by Floryan & Rowley (2018) from a more general theory as the mass ratio R is varied. The frequency k r0 given by (3.5) is also included, which obviously is a good approximation for large R, when the fluid-structure interaction becomes negligible. This agreement between k r obtained here by minimizing |F 2 | and that obtained numerically by Floryan & Rowley (2018) from a more general theory is a good test of the present algebraic approximation, especially considering that peak thrust forces are produced by the flexible foil at, or near, natural frequencies of the system, as shown by Floryan & Rowley (2018) and corroborated by our results in § 4 below.
Once the flexural deflection is obtained, one may compute C L and C M using (2.24) and (2.25). These quantities are needed to obtain C Li and C Mi from (2.9) and (2.10), respectively, necessary to compute the input power coefficient, which in turn is needed to obtain the propulsive efficiency. Note from (2.9) and (2.10) that this expression formally coincides with the usual one, for R = 0. (Note also that here the sign of C M has been changed in relation to Alaminos-Quesada & Fernandez-Feria (2020) and that, according to (2.8) and contrary to the moment, the sign of α is defined positive when clockwise to follow the usual convention in aerodynamics.) Additionally, one may obtain the propulsion, or thrust, coefficient C T (t) from the resulting deflection d and the given heaving and pitching kinematics, as done above for C L (t), C M (t) and C F (t) using the quartic approximation for the deflection. In the case of a quadratic deflection like (2.7), C T (t) is given by (4.5) in Alaminos-Quesada & Fernandez-Feria (2020). For the quartic deflection (2.8) the thrust coefficient is derived in a similar way, but the resulting expression becomes cumbersome. However, we have compared the time-averaged thrust values obtained with both approximations for some particular cases and found that differences are quite small, provided that the deflection d is computed with the quartic approximation, as described above. For this reason, we use the expressions for C T in Alaminos-Quesada & Fernandez-Feria (2020) in all the reported results.
The propulsive (Froude) efficiency is defined as where the bar denotes time-average over a period of the foil's harmonic motion, and similarly forC P . The mean value of the thrust coefficient C T (t) can be written as given in Alaminos-Quesada & Fernandez-Feria (2020), but with p = a and d m /(1 − a) 2 substituted by d m , which are reproduced in appendix B for easy reference. In this expression use has been made of the non-dimensional parameters where θ is the well-known Lighthill's (1969) feathering parameter. For θ hd = 0, we have the mean thrust coefficientC 0 T of a rigid foil with the same heaving and pitching motion (Fernandez-Feria 2016. From (2.9)-(2.10) and (2.15a-c)-(2.28) substituted into (3.6), the time-averaged power coefficient can be written as C P (kh 0 ) 2 = p h + p hp θ + p p θ 2 + p hd + p pd θ θ hd + R p rhd + p rpd θ θ hd where the functions p h (k), p hp (k, a, φ), p p (k, a), p hd (k, a, ψ), p pd (k, a, φ, ψ), p rhd (k, a, ψ), p rpd (k, a, φ, ψ), p shd (k, a, ψ) and p spd (k, a, φ, ψ) are also given in appendix B.    figure 3 but as a function of reduced frequency k and mass ratio R for a plate with stiffness S = 10. The dashed lines correspond to the first natural frequency of the system k r that minimizes |F 2 |, while the dashed-and-dotted line is k r0 (∝ R −1/2 ) given by (3.5), valid for for R 1.
Note that the first three terms constitute Theodorsen's power coefficient for a pitching and heaving rigid foil (Theodorsen 1935 . The remaining terms are the contributions to the input power of the passive flexural deformation associated with the prescribed heaving and pitching motions (i.e. terms containing d in C L and C M in the contribution of −ḣC L + 2αC M to C P ), and those associated with the inertia (R) and with the flexibility (S) of the foil.

Purely heaving motion with chordwise deflection
In the simplest case when no pitching motion is imposed on the foil (a 0 = 0), the expressions (3.3)-(3.4), (3.9) and (3.11) become particularly simple. Figures 3-5 show d m /h 0 = kθ hd and ψ as functions of S and k, R and k, and a and k, respectively, for particular values of the remaining parameter in each case (a = −1, S = 10 and R = 0.01).
For sufficiently large S and k in figure 3, or sufficiently large R and k in figure 4, the amplitude of the flexural component of the deflection d m shows a marked maximum that approximately follows the first natural mode of the system, which is also plotted in these figures when obtained analytically in the present approximation by minimizing |F 2 (R, S, k, a)|. This relation between resonance and the maximum of the flexural deflection, and then with the maximum of the thrust (see below), is in agreement with previous results of Floryan & Rowley (2018). In the present analysis, only the first natural frequency is captured since the foil's flexibility has only one degree of freedom, with a quartic deflection that approximates the shape of the first mode of an Euler-Bernoulli beam. Higher resonances are obtained by Floryan & Rowley (2018), among other authors, following the work of Wu (1961), where a large number of degrees of freedom are considered by expanding the deflection of the foil, and the corresponding inviscid motion around it, in a Chebyshev series with N coefficients (see also Alben 2008; Moore 2015, 2017; Tzezana & Breuer 2019). The present is a lowest-order model to account for the flexibility of the foil. In spite of this limitation, the great advantage of the present analysis is that an analytical solution is obtained. This situation is somewhat analogous to that considered by Moore (2014) for the case of a plate with infinite stiffness and finite torsional stiffness at the leading edge, where the variable is the pitching angle a 0 instead of the flexural deflection amplitude d m , which is obtained, also analytically, in terms of the finite torsional stiffness of the leading edge instead of the finite foil's stiffness S considered here. The first resonant mode of the system is, in addition, the most relevant one physically because the following modes are activated at much higher frequencies. Moreover, as shown by Alben (2008), the maximum possible thrust coefficient is always achieved by the flexible foil operating at or near this first resonance. Figure 3(a) also includes the natural frequency of the second mode of vibration extracted from figure 15 of Floryan & Rowley (2018) (note that their f * is the present k/π), above which the present approximation is certainly not valid.
In figure 4, the first resonance mode for large mass ratio R is proportional to R −1/2 , in accordance with k r0 given by (3.5) for constant S. Figure 5 shows that the maxima of the flexural component of the deflection amplitude for the constant values of S and R considered are obtained for pivots located just downstream of the mid-chord point.
Except for a pivot point close to the trailing edge, the local maxima in the amplitude of the flexural component of the deflection roughly correspond to local maxima in thrust and in power coefficients, since, according to (3.9)-(3.11), the thrust and power input are quadratic and linear functions, respectively, of d m . Although this is not always so because the effect of the phase, figures 6 and 7, where we plot in the stiffness-frequency plane the ratios of the time-averaged thrust and power coefficients of a flexible foil to the mean thrust and power coefficients of an otherwise identical rigid foil, respectively, clearly show that, for the selected values of a, the maxima inC T andC P approximately follow the first natural frequency, as does d m . We select two different mass ratios, R = 0.01 and R = 1, and three different locations of the imposed heaving motion, the leading edge (a = −1), the quarter-chord point (a = −1/2) and the mid-chord point (a = 0). Though we use h 0 = 0.02 in all reported computations, except otherwise specified, the results plotted in these figures are independent of h 0 in the present linear theory.
In all cases plotted in figure 6, or in figure 7, there exists a region in the stiffness-frequency plane where the thrust, or the input power, of the flexible foil is substantially larger than that of the rigid foil counterpart. As aforementioned, these regions are always centred around the first resonance mode of the system (plotted in figures 6 and 7 ln(C P /C P 0 ), R = 0.01, a = 0 ln(C P /C P 0 ), R = 1, a = 0 ln(C P /C P 0 ), R = 1, a = -1/2 ln(C P /C P 0 ), R = 0.01, a = -1/2 ln(C P /C P 0 ), R = 1, a = -1 ln(C P /C P 0 ), R = 0.01, a = -1 The regions of significant thrust enhancement by flexibility in figure 6 always lie in the large stiffness part of the plane, where the present small deflection approximation is more likely to be valid. For the case plotted in figure 6(a) (R = 0.01 and a = −1), this region practically coincides with that in figure 9(a) of Floryan & Rowley (2018) for the first resonance mode of the system (remember that their reduced frequency is our k/π). Though Floryan & Rowley (2018) obtain results valid for higher resonance modes, the present ones are obtained analytically for given S, R and the kinematics parameters. In figures 6(a) and 7(a) we also include the natural frequency of the second mode of vibration for the case R = 0.01 with a = −1, extracted from figure 15 of Floryan & Rowley (2018), above which the present approximation is not valid.
For a given pivot x = a, the region of thrust enhancement displaces to lower frequencies and becomes somewhat narrower as the mass ratio R increases. On the other hand, as the pivot point displaces downstream along the foil for given R, the magnitude of the thrust enhancement by flexibility increases for the selected values of a in figure 6. For a = 0, a region of net drag is also found for frequencies higher than the first natural frequency, where the present approximation is not valid. Similar trends are found for the power coefficient shown in figure 7, which becomes much larger than that of the rigid foil counterpart around the first resonance mode. This is also in accordance with previous results of Floryan & Rowley (2018) for a = −1 and R = 0.01 (their figure 10(a) practically coincides with figure 7(a) around the first natural mode of the system).
To see these trends more clearly for the case ofC T , figure 8 shows the maxima in the thrust enhancement, (C T /C 0 T ) max , whereC 0 T is the thrust of an otherwise identical rigid foil, as the pivot point location is varied for a stiffness S = 10 and for three different mass ratios, together with the corresponding values of the frequency k opt , which practically coincides with k r . For R small and unity the trends are similar, but with the maxima for each a occurring at lower frequencies as R increases, and with a marked peak at a ≈ 0.1. For R = 10 the thrust enhancement is smaller and its maximum is reached at approximately a = 0.4, with almost no thrust enhancement for pivot point locations upstream of the mid-chord point.
The propulsive efficiency does not necessarily follow the same trends ofC T andC P because it is the ratio of two quantities with similar qualitative behaviours, with maxima at, or very close to, the first natural frequencies. In fact, it turns out that the efficiency is larger at low frequencies, almost independently of the stiffness ratio, which is also the behaviour of the rigid foil counterpart (now independently of S, of course). This trend of the efficiency is typical of the purely heaving motion in the linearized potential limit (e.g. Garrick's efficiency for a rigid heaving foil tends to its maximum value unity as k → 0). Obviously, this is not physically valid because viscous drag is neglected in the present theory, and it is especially relevant in relation to the thrust force at low frequencies, substantially reducing the efficiency at these small frequencies, that typically becomes negative due to the viscous drag.
However, the inviscid results for the efficiency may conveniently be corrected at low frequencies by just subtracting an offset drag C D0 to the thrust coefficientC T obtained with the linear potential theory (Mackowski & Williamson 2015;Fernandez-Feria 2017;Floryan et al. 2017). Offset drag C D0 basically depends on the heaving amplitude h 0 and  Figure 9. Efficiency η as a function of frequency and stiffness for a heaving foil with a = 0, R = 0.01 (a) and R = 1 (b), and using h 0 = 0.02 with C D0 = 0.05. The dashed lines correspond to the first natural frequency of the system k r that minimizes |F 2 |, while the dashed-and-dotted line is k r0 given by (3.5). Areas with η < −1 have been whited out. the Reynolds number, and to a lesser degree on the pivot point location, but for the small amplitudes and high Reynolds numbers for which the present theory is valid one may consider a constant value. Figure 9 shows the propulsive efficiency for a foil pivoting at the mid-chord point (a = 0) and two values of the mass ratio (R = 0.01 and 1) when an offset drag C D0 = 0.05 is used with h 0 = 0.02. It is observed that the efficiency remains positive at high frequencies, almost independently of the stiffness, a consequence of the displacement to higher frequencies of the maximum efficiency by the effect of drag. But the main effect of the offset drag on the efficiency is the appearance of a localized maximum at lower frequencies, close to the region with maximum thrust and input power around the first resonance mode of the system. The existence of these resonant peaks in efficiency as a consequence of drag is explained by Floryan & Rowley (2018), though, as forC T andC P , we only observe here the peak associated with the first resonant mode due to the present flexural approximation. It must be noted that choosing a different value of C D0 , within its physical range around the present selected value (Floryan et al. 2017), only modifies the efficiency at very low frequencies, leaving practically unchanged the optimal regions depicted in figure 9. If one computes the relative efficiency η − η 0 , where η 0 is the efficiency of an otherwise identical rigid foil, one finds that the maxima always lie in the region of negative, or very small, efficiencies, at lower frequencies than those corresponding to the first natural mode, so that flexibility practically does not increase the propulsive efficiency in the region of large thrust enhancement. Thus, this region of thrust enhancement by flexibility near the first natural frequency for high enough stiffness almost coincides with the region of maximum efficiency, but with no relevant improvement of the efficiency in relation to a rigid foil.
Based on this, in figure 10 we plot η at the first natural frequency k r , which is a good approximation to the maximum efficiency, as a function of the pivot point location a for S = 50 and several values of the mass ratio R, together with the corresponding natural frequencies. Observed is an unphysical singularity at a ≈ 0.1, also observed in figure 8 for the thrust, which is an artifact of the present linear potential model associated with a vanishing small power coefficient. Disregarding this singularity, for R = 0.01 the model predicts a maximum efficiency in this case (h 0 = 0.02 and C D0 = 0.05) of about 0.4 when the pivot point is close to the quarter-chord (a = −1/2). As R increases, the behaviour is similar, but with smaller efficiencies and a maximum attained for increasing a when a < 0. For a 0.1, the efficiency is negative or very small. In fact, the singularity roughly marks This behaviour is corroborated for all values of the stiffness S in figure 11, which shows η at k = k r in the S-a plane for R = 0.01 and R = 1, together with the corresponding values of k r . For −1 < a 0, the maximum efficiency for constant R is attained at values of the stiffness S that decrease almost linearly with a. Curiously, this linear band of maximum efficiency in the S-a plane practically corresponds to a constant value of k r . For R = 0.01, the maximum of η within this band is reached for a pivot point close to the the mid-chord with S ≈ 10 (disregarding the singularity at a ≈ 0.1 mentioned above, which has been whited out in figure 11, together with negative values of the efficiency), while for R = 1 the maximum is attained for a pivot point approaching the leading edge for a value of S larger than 100. In fact, a local maximum of η exists at a = −1 for both values of R at sufficiently high values of S (larger than 100). Given its physical relevance within the present approach, valid for large S, this local maximum is characterized in figure 12. This figure shows that the maximum efficiency for a foil pivoting at the leading edge in the present case is approximately 0.44, attained for R 10 and S ≈ 430 with an operating frequency k 3.13. Figure 11 also shows that for a 0.1 the efficiency is negative or very low, as already discussed in relation to figure 10 for S = 10, except for small S for which the present theory is not valid.

Purely pitching motion with chordwise deflection
For a foil with mass ratio R undergoing a purely pitching motion of small amplitude a 0 about a pivot point location a, the effect of flexibility on the propulsive performance is qualitatively similar in the stiffness-frequency plane to that described above for a heaving motion with the same a and R. The pitching motion just generates smaller maxima of the flexural deflection amplitude than that from the equivalent heaving motion, in relation to their rigid foil counterparts, and consequently less thrust enhancement by flexibility.
As an example, figure 13(a) shows, in the stiffness-frequency plane, the ratio of the time-averaged thrust of a flexible foil to the mean thrust coefficient of an otherwise identical rigid foil for R = 0.01, a = −1 and a 0 = 2 • , though the plotted results are independent of the pitching amplitude a 0 . The region of substantial thrust enhancement by flexibility is similar to that in figure 9(b) of Floryan & Rowley (2018) around the first resonance mode of the system (which is also plotted as a dotted line). Figure 13(b) shows the propulsive efficiency for the same case, but now using an offset drag C D0 = 0.05. As in the case of purely heaving motion, there is a region of maximum efficiency close to the thrust enhancement region around the first natural mode. Figure 14 shows the propulsive efficiency at k = k r as a is varied for S = 10 and several values of R, together with the corresponding values of k r . It is observed that the efficiency predicted by the model becomes singular in certain narrow intervals of a containing a = 0 that depends on R and S; within these ranges of a the present approximation with a small quartic deflection is not valid. For pivot locations upstream of these singular intervals, the efficiency reaches a maximum that decreases with increasing R at practically the same value of a ( 0.17 in the plotted cases). For pivot locations downstream of the singular narrow region, the efficiency is negative or very low. Therefore, the singularity marks again the pivot point locations downstream of which the propulsive efficiency of the flexible (now pitching) foil becomes poor. The difference from the heaving foil case discussed above is that now there exists a narrow region of pivot point locations that depends on S and R, instead of being an almost constant pivot point location (a ≈ 0.1) as S and R are varied. Indeed, this is observed in figure 15, which shows the effect of the stiffness S on the efficiency at k = k r for the case with R = 0.01 in figure 14. The singular regions and those with negative efficiency have been whited out in the a-S plane. As the pivot location moves from the leading edge to the singularity close to the mid-chord point, the maxima in the efficiency are attained at decreasing values of S, starting close to 100 for a = −1. For pivot point locations downstream of the singular region, the efficiency is negative or very low for all the range of values of S considered. Figure 16 shows the evolution with the mass ratio R of the maximum in the efficiency when the pivot point is at the leading edge (a = −1), together with the corresponding values of S and k r .

4.3.
Pitching and heaving motion about the leading edge Finally, as an example of combined pitching and heaving motions we only present results for the case in which they are forced at the leading edge (a = −1), and just to see how an appropriate addition of pitching to a heaving motion may significantly increase both the thrust and the propulsive efficiency of the foil, additionally to that caused by flexibility discussed above, as is well known for rigid flapping foils (Anderson et al. 1998).
To that end we select the case with the highest propulsive efficiency for a heaving foil with R = 0.01 when h 0 = 0.02, C D0 = 0.05 and a = −1, which according to figure 12 is approximately 40 % when k 6.9 and S about 130. Figure 17 shows the thrust and efficiency enhancement in relation to the corresponding heaving-only foil,C T /C Th and η − η h , respectively, when pitching motions of increasing amplitude a 0 and phase shift φ are added.
Any amount of pitching, even of infinitesimal amplitude, increases the thrust generated by the foil provided that its phase shift φ in relation to the heaving motion is appropriate. For small a 0 , this happens in this case when 90 • φ 270 • , with the thrust enhancement maximized at φ ≈ 180 • . In relation to the efficiency, figure 17(b) shows that there exists a threshold in the pitch amplitude above which the efficiency improves in relation to the heaving-only motion. For the present case with very small h 0 , this minimum amplitude is also very small, about 0.2 • , generating thrust enhancement for a phase shift φ of about 204 • . As the amplitude a 0 increases, the optimal phase shift for efficiency enhancement increases, approaching 360 • (i.e. 0 • ) for sufficiently high amplitude.

Concluding remarks
We derive analytical expressions for the propulsive performance of a two-dimensional foil with prescribed pitching and heaving motions about arbitrary pivot points when a quartic chordwise deflection is allowed to accommodate passively to the fluid-structure interactions. The formulation is limited to high Reynolds numbers (inviscid flow), small amplitudes of the foil motion (no flow separation) and large stiffness combined with low-to-moderate reduced frequencies and mass ratios, so that only the first resonant mode of the system is activated. The formulation incorporates the effect of the inertia of the solid foil on both the thrust and input power, and therefore on the propulsive efficiency, directly from the Euler-Bernoulli beam equation, an effect that becomes relevant when the mass ratio is not small. Thus, the results are valid for modelling both swimming and flying applications of flapping foils.
The present formulation uses a slightly modified form of the purely kinematic results in Alaminos-Quesada & Fernandez-Feria (2020) for a quadratic deflection to compute the flexural deflection amplitude and its phase in terms of the solid foil properties, namely its dimensionless stiffness and mass ratio, in addition to all the kinematic parameters associated with the forced heaving and pitching motion applied at an arbitrary pivot point. The maximum flexural deflection amplitude is found to be attained at, or near, the first natural frequency of the system, which is also obtained here algebraically, being the only natural frequency recovered in the present approximation. At this frequency, which depends on the stiffness, mass ratio and the location of the pivot point, the thrust and power input also reach their respective maximum values, approximately, though not necessarily the propulsive efficiency. These results are in agreement with previous numerical results for purely heaving and pitching motions about the leading edge (a = −1) for low mass ratios (R = 0.01) reported by Floryan & Rowley (2018), though the formulation by those authors accounts for the whole series of natural frequencies of the system, depending on the number of terms retained in the Chebyshev expansion. In this sense, the present approach is equivalent to the lowest order of this expansion formulated through the first three moments of the Euler-Bernoulli beam equation, thus constituting a minimal model of the flexible foil. Though the limitation to (small-amplitude) quartic deformations restricts the theory to include only the first resonant frequency of the system, as shown by Alben (2008) from a similar linear potential flow theory, the maximum possible thrust coefficient is always achieved by the flexible foil operating at or near this first resonance. The cases of purely heaving and purely pitching motions are analysed with detail in the present work, characterizing the parametric ranges generating maximum thrust and maximum propulsive efficiency enhancement as the stiffness of the foil is varied, including the effect of the pivot point location for a wide range of mass ratios.
Though results are presented for a foil with uniform thickness, chordwise rigidity and density, the analytical results are easily generalized to account for any chordwise distribution of the non-dimensional stiffness S(x) and the mass ratio R(x) using (2.4)-(2.6). More complicated motions of the foil, with additional kinematic parameters, could also be considered within the present approach by including additional moments of the Euler-Bernoulli equation, though then the resulting analytical expressions, if attainable, would become much more involved and probably would not have any advantage over the method based on Chebyshev expansions.