Unsteady Propulsion of a Two-Dimensional Flapping Thin Airfoil in a Pulsating Stream

Unsteady Propulsion of a Two-Dimensional Flapping Thin Airfoil in a Pulsating Stream Ernesto Sanchez-Laulhe∗and Ramon Fernandez-Feria† University of Málaga, 29071 Málaga, Spain and Anibal Ollero‡ University of Seville, 41092 Seville, Spain https://doi.org/10.2514/1.J062925 The cruising velocity of animals, or robotic vehicles, that use flapping wings or fins to propel themselves is not constant but oscillates around a mean value with an amplitude usually much smaller than the mean, and a frequency that typically doubles the flapping frequency. Quantifying the effect that these velocity fluctuations may have on the propulsion of a flapping and oscillating airfoil is of great relevance to properly modeling the self-propelled performance of these animals or robotic vehicles. This is the objective of the present work, where the force and moment that an oscillating stream exerts on a two-dimensional pitching and heaving airfoil are obtained analytically using the vortical impulse theory in the linear potential flow limit. The thrust force of the flapping airfoil in a pulsating stream in this limit is obtained here for the first time. The lift force and moment derived here contain new terms in relation to the pioneering work by Greenberg (1947), which are shown quantitatively unimportant. The theoretical results obtained here are compared with existing computational data for flapping foils immersed in a stream with velocity oscillating sinusoidally about a mean value. I. Introduction ARENEWED interest in more precise models of the aerodynamic forces on oscillating foils, considering all the unsteady processes intervening in the fluid–foil interaction, has been motivated by the growing interest in developing biologically inspired aerial and aquatic vehicles self-propelled by flapping wings or fins [1–6]. Most aerodynamic models come from studies that consider the oscillating foil in a flow with constant freestream velocity. However, the cruising velocity of animals, or robotic vehicles, that use flapping wings or fins to propel themselves is not constant [3,7–11]. It oscillates around a mean value with an amplitude usually much smaller than the mean, and a frequency that typically doubles the flapping frequency. For instance, for a small flapping amplitude (airfoil chord length ratio, say, ϵ≪1), the cruising velocity oscillations amplitude has been found to be of the order of ϵ2∕3times the mean velocity for a simple aquatic locomotion model [11], and of the order of ϵ2for a simple ornithopter flight model [10]. The frequency doubling is nicely illustrated by the ability of a sinusoidally plunging airfoil to produce thrust by the Knoller–Betz effect [12]. Although it has been argued that the oscillations in the flow speed have little impact on the time-averaged aerodynamic force [13], and therefore on the time-averaged cruising velocity generated by an oscillating foil, results obtained in the present work that are compared with previous numerical results show that these oscillations may modify substantially not only the instantaneous aerodynamic forces, but also their time-averaged values. Thus, the oscillations in the cruising velocity may be relevant for the unsteady dynamics of the robotic vehicle propelled by the flapping foil, and therefore for its guidance and control. The objective of the present work is to include and quantify this effect of a periodic streamwise velocity about a mean value on the time-dependent thrust force generated by a pitching and heaving foil in the limit of small flapping amplitudes and high Reynolds numbers, so that linear potential flow theory can be applied, a limit of interest for the cruising regime of vehicles self-propelled by flapping wings or fins. The effect of a pulsating stream on the airfoil aerodynamic forces in the linear potential and two-dimensional (2D) flow limit was considered by Isaacs [14] and Greenberg [15], partially inspired by unsteady problems related to helicopter aerodynamics. Isaacs [14] analyzed the general problem of an airfoil at constant angle of attack with a variable stream velocity, deriving an explicit formula for the lift in the particular case of a pulsating stream, consisting of a constant velocity plus a sinusoidal variation of any amplitude and frequency. The work by Greenberg [15] extended that of Theodorsen [16] for the unsteady lift force and moment on a pitching and heaving foil in a uniform current by considering a pulsating stream, adding new terms to Theodorsen’s lift and moment containing time derivatives of the streamwise velocity. However, Greenberg did not consider the thrust force, and, as we shall be see below, an additional frequency term was missing in Greenberg’s lift and moment. Here we shall use the linearized vortical impulse theory, already used in the pioneering work by von Kármán and Sears [17], who reproduced Theodorsen’s lift force and moment for a pitching and heaving foil in a uniform stream. The impulse theory was generalized for an incompressible flow by Wu [18], and recently applied to obtain the thrust force generated by a pitching and heaving foil in a uniform stream in the linear limit [19]. This vortical impulse theory is used here to obtain the two components of the force and the moment exerted by a 2D sinusoidal head-on stream on a pitching and heaving foil, recovering Greenberg’s expressions for the lift and moment, but with additional terms, and obtaining a new expression for the thrust force. We also include a nonvanishing mean angle of attack in the formulation, to cover the case of a steady foil in a streamwise pulsating freestream, and also the effect of that mean angle of attack in the expressions for the lift, moment, and thrust of a flapping foil, not included in previous theoretical formulations. The aerodynamic performance of a 2D airfoil in a streamwise sinusoidal flow has been analyzed numerically by Lian and Shyy [20,21] as a simple model to investigate the performance of stationary and flapping wings under this particular class of gusty environments. These authors found that flapping wings can alleviate the undesirable gust effect on the lift force better than an otherwise identical stationary foil for certain flapping kinematics. These numerical results are used here to comparewith the theoretical results derived in the present work. In summary, the novelty of the present work is to use the linearized vortical impulse theory to derive general analytical expressions for Received 20 February 2023; revision received 12 May 2023; accepted for publication 18 May 2023; published online 26 June 2023. Copyright © 2023 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved. All requests for copying and permission to reprint should be submitted to CCC at www.copyright.com; employ the eISSN 1533-385X to initiate your request. See also AIAA Rights and Permissions www.aiaa.org/ randp. *Graduate Student, Fluid Mechanics Group; [email protected]. †Professor, Fluid Mechanics Group; [email protected]. ‡Professor, GRVC Robotics Laboratory; [email protected]. Article in Advance / 1 AIAA JOURNAL Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 Preprint version of: Sanchez-Laulhe E., Fernandez-Feria R., and Ollero A. (2023), Unsteady Propulsion of a Two-Dimensional Flapping Thin Airfoil in a pulsating stream In AIAA Journal the lift, thrust, and moment of a pitching and heaving thin airfoil with a fluctuating cruising velocity. The lift and moment expressions recover previous ones, but with new terms that are shown to be quantitatively of small relevance in relation to all the other additional terms relating to the freestream fluctuations. The main result is a new general analytical expression for the thrust force, which is essential to model properly the propulsive performance of animals or robotic vehicles that use flapping wings or fins to propel themselves. Another novelty of the present work not included in previous theoretical formulations is the effect of a nonvanishing mean angle of attack in the derived general expressions for the lift, thrust, and moment, thus increasing their ranges of applicability. The problem is formulated in the next section. General expressions for the lift, thrust, and moment are derived in Sec. III, and their corresponding time averages are given in Sec. IV. Comparison with published numerical results is presented in Sec. V, concluding in Sec. VI with a brief summary and scope of the results. II. Formulation We consider a 2D and slender rigid foil of chord length cundergoing harmonic pitching and heaving motions with circular frequency ωand advancing in a fluid of density ρwith a horizontal velocity that also oscillates harmonically, but now with a frequency ω1, around a constant value Us. We use nondimensional quantities scaled with the half-chord length c∕2, the velocity Us, and the density ρ. The nondimensional heaving and pitching motions, and the nondimensional velocity in the −xdirection (see Fig. 1) are written, respectively, as hth0Reikt;αtαsα0Reiktϕ; Ut1σReik1tϕ1(1) where Rmeans real part; kωc 2Us and k1ω1c 2Us (2) are the nondimensional, or reduced, frequencies; tis the nondimensional time, scaled with c∕2Us;h0,α0, and σare the nondimensional amplitudes of the harmonic heaving, pitching, and horizontal motions, respectively; and ϕand ϕ1are the phase shifts of these two last motions in relation to the foil’s heave. Also added to the formulation is a constant mean angle of attack αs. We assume that the Reynolds number based on Usand cis large enough and that h0≪1,jαsj≪1, and α0≪1, so that the 2D linear, potential flow theory can be applied. With the last assumptions, i.e., for jαtj ≪1, the foil approximately lies in the interval of the xaxis −1≤x≤1in a reference frame moving with it along the xdirection. The pivot axis is located at an arbitrary point xa(see Fig. 1). The amplitude σ≡Uσ∕Usof the horizontal velocity oscillations will also be assumed small (see Sec. III.C below for the exact assumption on σ). In this reference frame, the nondimensional vertical displacement and velocity of the foil are zsx; tht−x−aαt;−1≤x≤1 (3) v0x; t−Utαt_ ht−x−a_ αt(4) respectively, where a dot denotes differentiation with respect to t. Note that hand v0are both positive upwards, while αis positive clockwise. III. Unsteady Lift, Thrust, and Moment The vortical impulse theory of an incompressible and unbounded flow [18,22,23] is used to obtain the forces and moment on the foil. In the linearized potential flow limit, the vorticity is concentrated at the foil interface and the trailing wake, both considered as vortex sheets along the plane z0. The expressions for the lift force and moment from this formulation were first derived by von Kármán and Sears [17], who reproduced previous results by Theodorsen [16] using a more standard linear potential flow theory, while the thrust (or minus the drag) force was derived using the vortical impulse theory in [19]. In dimensionless form they can be written as CLt Lt 1∕2ρU2 sc−d dt Z1 −1 xϖsdxZ∞ 1 xϖedx(5) CTt Tt 1∕2ρU2 sc−d dt Z1 −1 zsϖsdxZ∞ 1 zeϖedx(6) CMt Mt 1∕2ρU2 sc2 1 4 d dt Z1 −1x−a2ϖsdxZ∞ 1x−a2ϖedx(7) where Ltand Ttare the forces components (per unit span) in the directions zand −x, respectively, and Mtis the moment (per unit span and positive when clockwise) in relation to the pitching axis xa. In the last expression it must be taken into account that da∕dt Utin a stationary reference frame with the fluid at rest far from the foil. In these expressions, ϖsx; t,−1≤x≤1, is the nondimensional vorticity density distribution on the foil; ϖex; tis the nondimensional vorticity density distribution in the trailing wake, both scaled with Us; and zex; tis the vertical position of each point of the vortex wake. We consider the long-time behavior in which the vortex wake sheet extends many chord lengths downstream of the airfoil, so that, in first approximation, 1≤x≤∞for both ϖex; t and zex; t. As aforementioned, jzsj≪1and jzej≪1. Following von Kármán and Sears [17], the contribution from the vortex-sheet wake to ϖsx; t, denoted by ϖ1, is separated from the bound circulation ϖ0that would produce the motion of the foil as if the wake had no effect: ϖsx; tϖ0x; tϖ1x; t;−1≤x≤1 (8) with Γ0tZ1 −1 ϖ0x; tdx(9) being the nondimensional circulation (scaled with Usc∕2) that would be obtained from the quasi-steady foil theory, without unsteady wake. The details of the derivation of the vorticity distributions ϖ0,ϖ1, and ϖecan be found in Ref. [17] (see also Refs. [19,24]). Just a summary is given here to better appreciate the differences introduced by the unsteady stream velocity Ut. A. Vorticity Distribution The contribution ϖ1from the unsteady planar wake ϖe, after applying Kutta’s condition at the trailing edge x1, is (using the dummy variable ξfor the integration [17]) ϖ1x; t1 πZ∞ 1 ξ1 ξ−1 s 1−x 1x rϖeξ;t ξ−xdξ(10) Fig. 1 Schematic of the pitching and heaving foil (nondimensional quantities). 2Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 Kelvin’s total-circulation conservation theorem requires that Γ0tΓ1tZ∞ 1 ϖex;tdx0;with Γ1tZ1 −1 ϖ1x;tdxZ∞ 10 @ ξ1 ξ−1 s−11 Aϖeξ;tdξ(11) so that it constitutes a general relation between Γ0tand ϖex; t: Γ0tZ∞ 1 ξ1 ξ−1 sϖeξ;tdξ0 (12) The vorticity ϖ0is obtained writing the velocity field induced by the whole vorticity distribution, and imposing that the normal velocity on the foil (i.e., at z0,−1≤x≤1, in the present linear approximation) is the vertical velocity (4). After applying Kutta’s condition at the trailing edge to fix Γ0t, we obtain ϖ0x; t 1  1−x2 pΓ0 π2_ ha_ α−Uαx_ α1–2x2(13) Γ0t2πUα−_ h−a−1 2_ α(14) Once the vorticity distribution is known in terms of the foil motion and ϖex; t, which is related to Γ0t(and therefore to the foil motion) through Eq. (12), general expressions can be obtained for Eqs. (5–7), taking into account that the wake is convected downstream with velocity Ut[17]: ϖeξ;tϖeX;z eξ;tzeX;with Xξ−Zt ti Utdt(15) so that both ϖeand zeremain constant in a reference frame where the foil moves with velocity Ut; i.e., the fluid at infinity is at rest. Thus, for any function fξsatisfying f10, d dt Z∞ 1 fξϖeXdξUtZ∞ 1 dfξ dξϖeXdξ(16) This property is valid even for the present case where Udepends on time. B. General Expressions for CL,CM, and CT Starting with the lift coefficient, after substituting the above vorticity distributions into Eq. (5), it can be written as a sum of three terms: CLtCL0tCL1tCL2t(17) with CL0t−UtZ∞ 1 ξ1 p  ξ−1 pϖeξ;tdξUtΓ0t(18) CL1t−d dt Z1 −1 xϖ0x; tdxπ_ UαU_ α− h−a α(19) CL2tUtZ∞ 1 ϖeξ;t  ξ2−1 pdξ(20) It is observed that, in addition to the fact that now the velocity Utis not constant (is not unity in dimensionless form), an additional term related to its temporal derivative appears in the added mass lift (19) in relation to Theodorsen’s (or von Kármán and Sears’) lift coefficient. This new term was already obtained by Greenberg [15] using Theodorsen’s approach instead of von Kármán and Sears’one used here. Additional terms related to the temporal variation of Uwill appear once the wake vorticity ϖeis resolved below for the harmonic motion (1). But before that it is instructive to write down the general expressions of CTtand CMtin terms of ϖex; t, as just done for CLt. Following with CMt, substituting Eqs. (8), (10), and (12–14) into Eq. (7), it can be written as CMtCM0tCM1tCM2ta 2CLt(21) where CM0t1 2UtZ1 −1 xϖ0x; tdx1 4UΓ0−πU_ α π 2UUα−_ h−a_ α(22) CM1t−1 4 d dt Z1 −1x2−1 2ϖ0x; tdxπ 16  α(23) CM2t−1 4UtZ∞ 1 ϖeξ;t  ξ2−1 pdξ−1 4CL2t(24) This expression is formally the same as that derived by Theodorsen, except that now Udepends on time (and, of course, the additional term proportional to _ Uappearing in CLcommented on above). The derivation of the thrust CTtis a bit more complex because it involves the vertical position of the wake ze[19]. Taking into account the property (15), the position of the wake follows the successive locations of the foil’s trailing edge in a reference frame where the fluid at infinity is at rest. Thus, on using Eq. (16), d dt Z∞ 1 zex; tϖex; tdxUtht−1−aαtϖex1;t (25) where ϖex1;tis the wake vorticity evaluated at the trailing edge (x1) at any instant of time t. The other term in Eq. (6) can be written as (for simplicity we do not write explicitly the dependencies on tand/or xof the different variables) d dt Z1 −1 zsϖsdxd dt haαZ1 −1 ϖ0dx−αZ1 −1 xϖ0dx d dt haαZ1 −1 ϖ1dx−αZ1 −1 xϖ1dx _ ha_ αΓ0haαdΓ0 dt −_ αZ1 −1 xϖ0dx αCL1_ ha_ αΓ1haαdΓ1 dt −_ αZ1 −1 xϖ1dx−αd dt Z1 −1 xϖ1dx(26) Taking into account that dΓ1 dt −dΓ0 dt −Uϖex1;t; Z1 −1 xϖ1dxZ∞ 1 ξ2−1 p−ξϖedξ(27) d dtZ1 −1 xϖ1dxUZ∞ 1 ξ  ξ2−1 pϖedξ−UZ∞ 1 ϖedξ−Uϖex1;t (28) Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO 3 Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 and the above expressions for Γ0,Γ1, and CL, the general expression for CTcan be written as CTt−αtCLtCT1tCT2t(29) with CT1t_ αZ1 −1 xϖ0dxπ_ α_ ha_ α−Uα(30) CT2tZ∞ 1h_ ha_ α−αU_ α ξ2−1 p−ξiϖedξ(31) This expression formally coincides with that derived in Ref. [19], but now Udepends on t, with an additional term proportional to _ Uin CL mentioned above and with new terms in ϖeassociated to the time dependence of Ut, as discussed next for the harmonic motion (1). It is worth mentioning here that the expression for the added-mass term CT1in the right-hand side of Eq. (30), derived in Ref. [19], coincides with that recently obtained by Limacher [25] for a thin airfoil, being absent in the classical thrust expression by Garrick [26]. C. Expressions of CL,CM, and CTfor the Harmonic Motion (1) Substituting Eq. (1) into Eq. (14), Γ0tcan be written as Γ0t2πRαs−ikh0α0eiϕ−a−1 2ikα0eiϕeikt αsσeiϕ1eik1t1 2α0σeiϕϕ1eik2t1 2α0σeiϕ−ϕ1eik3t; k2kk1;k 3k−k1(32) where, in addition to the frequencies kand k1, two new frequencies k2 and k3appear associated to the product term Uα. With this circulation, Kelvin’s theorem (12) and the property (15) yield the wake vorticity ϖe. It depends on the variable Xξ−Zt ti Utdtξ−Zt ti1σcosk1tϕ1dt ξ−t−σ k1 sink1tϕ1constant (33) Since σis assumed small, following Greenberg [15] we also assume that σ k12Uσ ω1c≡k−1 σ≪1 (34) i.e., the reduced frequency associated to the frequency and amplitude of the pulsating stream velocity, kσ, is very large, so that the wake’s phase velocity in relation to the foil is basically that of the uniform flow, X≃ξ−tconstant. Therefore, the wake vorticity can be written as the sum of five terms, corresponding to the five different dependencies with tof Γ0tin Eq. (32): ϖeξ;t−2παsδξ−∞geikt−ξg1eik1t−ξ g2eik2t−ξg3eik3t−ξ(35) where the first term corresponds to the starting point vortex far from the foil, δbeing Dirac’s delta function, and the constants of the following terms are given, according to Eqs. (12) and (32), by g−2πf−ikh0α0eiϕ−a−1∕2ikα0eiϕg R∞ 1 ξ1∕ξ−1 pe−ikξdξ 4f−ikh0α0eiϕ−a−1∕2ikα0eiϕg iH2 0kH2 1k(36) g1−2παsσeiϕ R∞ 1 ξ1∕ξ−1 pe−ik1ξdξ4αsσeiϕ iH2 0k1H2 1k1 (37) g2−πα0σeiϕϕ1 R∞ 1 ξ1∕ξ−1 pe−ik2ξdξ2α0σeiϕϕ1 iH2 0k2H2 1k2 (38) g3−πα0σeiϕ−ϕ1 R∞ 1 ξ1∕ξ−1 pe−ik3ξdξ2α0σeiϕ−ϕ1 iH2 0k3H2 1k3 (39) One has to take the real part of these expressions, where H2 nz Jnz−iYnz,n0;1, are the Hankel functions of the second kind, related to the Bessel functions of the first and second kind Jn and Yn, respectively [27]. It must be noticed that in the Hankel functions in g3we actually use jk3jinstead of k3, with ϕ−ϕ1when k3>0, and ϕ1−ϕwhen k3<0, in accordance with Eq. (32). Substituting this wake vorticity into the general expressions for CLt,CTt, and CMtof the above Sec. III.B, the following expressions are obtained: CLtπ_ UαU_ α− h−a α2πUαs URG0Ckeikt G01Ck1eik1tG02Ck2eik2t G03Ck3eik3t(40) CMtπ 2a_ Uαa−1 2U_ α−1 8a2 α−a hπa1 2Uαs a1 2U 2RG0Ckeikt G01Ck1eik1tG02Ck2eik2t G03Ck3eik3t(41) CTt−αCLπ_ α_ ha_ α−Uα−_ ha_ α−Uα2παs −_ ha_ α−UαR2i πG0C1keikt 2i πG01C1k1eik1t 2i πG02C1k2eik2t2i πG03C1k3eik3t −_ αRG0−2 πk1ikC1k−i kCkeikt G01−2 πk11ik1C1k1−i k1 Ck1eik1t G02−2 πk21ik2C1k2−i k2 Ck2eik2t G03−2 πk31ik3C1k3−i k3 Ck3eik3t(42) where Cz≡FziGz H2 1z iH2 0zH2 1z; C1z≡F1ziG1z 1∕ke−iz iH2 0zH2 1z(43) with Czthe well-known Theodorsen function [16,26], now applied to different reduced frequencies (zk; k1;k 2or k3), and G02π−ikh0α0eiϕ−ikα0eiϕa−1 2 (44) 4Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 G01 2πσαseiϕ1;G 02 πσα0eiϕϕ1; G03 πσα0eiϕ−ϕ1(45) taking into account the above comment on the signs of k3and ϕ−ϕ1 in Ck3and G03. Clearly, expressions (40–42) coincide with the previously known ones for a constant stream velocity [16,17,19] by setting σ0(i.e., U1and G01 G02 G03 0), except for the fact that we have also included the contributions from a nonzero mean angle of attack αs, which appear in the above expressions just after the corresponding added-mass terms. For a steady airfoil (h0α00) in a uniform current (σ0) one obviously recovers the classical expressions CL2παs,CMa1∕2CL∕2, and CT0. The expressions for the lift and moment (40) and (41) basically coincide with those obtained by Greenberg [15] from a different approach within the linearized potential flow theory (following the work by Theodorsen [16]), except for the terms with k3in Eqs. (40) and (41), i.e., the terms multiplied by G03, which are missing in Greenberg’s expressions for the lift and moment. Probably, when considering the products of pitching and pulsating-stream terms, this author only considered the sum of frequencies, missing the frequency difference terms. However, as shown in Fig. 2 for a representative case of a pitching and heaving airfoil in a pulsating stream with σ 0.2 that will be used below in the validation section, the differences in the results when these k3terms are neglected are quite small. On the other hand, Greenberg did not obtain the thrust force, derived here together with the lift and moment, though the present work has single frequency for pitching and heaving, whereas Greenberg had separate frequencies for each. The results for a uniform freestream (σ0), i.e., Theodorsen’s lift and moment for the same pitching and heaving motion, are also plotted in Fig. 2 (as well as in many figures in Sec. V below) to emphasize the substantial effect that the oscillating freestream may have on the aerodynamic forces and moment, not only on their instantaneous values but also on their time averages. The contributions of the new k3terms in CLand CMare negligible in comparison in all the cases we have considered. In particular, in the example plotted in Fig. 2, the maximum relative error in the instantaneous values is about 9% for CL, with an average relative error of just 5% for CLand slightly smaller for CM. Another quantity of interest is the input power, or energy transferred from the fluid to the foil per unit time and unit span, which in dimensionless form can be written as CPit−CLt_ h−2CMt_ α(46) D. Stationary Airfoil in a Pulsating Stream A relevant particular case of the above expressions is that of a stationary airfoil with a (small) angle of attack αsimmersed in the pulsating stream Utdefined in Eq. (1) (but now with ϕ10since there is no other harmonic motion). Remember that the present linear potential flow theory is approximately valid for jαj≪1, implying that the results are approximately valid for jαsj≪1, in the same sense that 2παsis the approximate lift coefficient of a stationary airfoil in a uniform freestream for small angle of attack αs. Setting ααs, h00, and α00in Eqs. (40–42) and rearranging some terms, the coefficients for a stationary airfoil in a pulsating stream can be written as CLtπ_ Uαs2πUαs1σRfCk1eik1tg (47) CMtπa _ U 2αsπa1 2Uαs1σRfCk1ek1tg (48) CTt−π_ Uα2 s−2σUα2 sRfπCk1−2iC1k1eik1tg(49) This problem was also considered by Isaacs [14] from the linearized potential flow theory, but without assuming small variations of the stream velocity, i.e., without the assumption (34), and evaluating only the lift. In fact, Isaacs compared his general lift expression with the simplest, quasi-steady one where the lift is assumed proportional to the square of the unsteady stream velocity, showing with a particular example of interest in rotary wing aircraft that the results are not very different. This case, which corresponds to k10.0424 and σ0.4, is considered in Fig. 3, comparing the present result (47) for CL(divided by the mean 0 0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1 1.2 1.4 1.6 Present Greenberg ([15]) =0 a) 0 0.2 0.4 0.6 0.8 1 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Present Greenberg ([15]) =0 b) Fig. 2 CLt(a) and CMt(b) from Eqs. (40) and (41) compared with the same expressions without the k3terms (Greenberg): h00.5,α08.42°, ϕ90°, k0.25,αs8°, σ0.2,k1k∕10, and ϕ1−90°. Also plotted results for a uniform stream (σ0). 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 Present Isaacs ([14]; exact) Isaacs ([14]; approx.) Fig. 3 CLtfrom Eq. (47) normalized with 2παscompared with Isaac’s exact and approximate (U2) solutions [14] for σ0.4and k10.0424. Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO 5 Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 2παs) with both the exact and the approximate solutions given in Isaacs’s final numerical example in Ref. [14]. Note that 2π∕k1is the nondimensional period of the pulsating stream. It is observed that CLtfrom Eq. (47) is very close to Isaacs’s exact solution, despite the fact that σ∕k1is not small, significantly closer than Isaacs’s approximate solution, which in the present notation is just U2t. IV. Time-Averaged Quantities In many situations one is more interested in the above expressions time-averaged. Since there are different frequencies involved, one cannot average over just one of the periods. An alternative is to retain only the nonperiodic terms in the expressions, which will coincide with the time average over a cycle in pitching and heaving foils in a uniform stream. The quantities thus obtained are denoted in the expressions below with an overbar. Only the terms containing σ and/or αsare written explicitly, using a subscript 0 to denote the corresponding values in a uniform stream with αs0(notice that  CL0 CM00) [16,19,28].  CL2παs1σ2Fk1 2(50)  CM1 2a1 22παs1σ2Fk1 2(51)  CT CT0−σαs2πFk12G1k1 −σα02 4πFk22G1k2πFk32G1k3 (52)  CPi  CPi0−πkσ2h0α0 4fsinϕFk2Fk3cosϕGk2 Gk3g−a1 2πkσ2α2 0 4Gk2Gk3 (53) Another quantity of interest is the (Froude) efficiency, η  CPo  CPi (54) where the nondimensional, time-averaged power output coefficient is (notice that  CPo0 CT0)  CPo  UCT  CT−πσ2α2 sFk1−πσ2α2 0Fk 2Fk2 4Fk3 4 kGk 2a−1 2−2σ2α2 sG1k1−σ2α2 0G1k −πkσ2h0α0 2Gkcos ϕ−Fksin ϕ−kσ2h0α0 2F1k2 F1k3cos ϕ−G1k2G1k3sin ϕ −akσ2α2 0 2F1k2F1k3−σ2α2 0 2G1k2G1k3 kσ2h0α0F1kcos ϕG1ksin ϕ kσ2α2 0F1ka−1 2kσ2α2 0 42F1k22F1k3 2G1k2 k22G1k3 k3πFk2 k2πFk3 k3(55) It is remarkable that the pulsating flow has no effect on the timeaveraged coefficients for pure heave (α00) with αs0. Some of the above time-averaged expressions are not valid when kk1,k12k, and k2k1. These special cases are summarized in the Appendix. V. Validation As a validation of the above theoretical results we compare them with the numerical results by Lian and Shyy [20,21]. These authors considered the effect of a sinusoidal head-on gust Utlike that in Eq. (1) (with fluctuation amplitudes σ0.2 and 0.5 and several values of k1and ϕ1) on the aerodynamic performance of a NACA0012 airfoil at Reynolds number 4×104, both for a stationary airfoil with a given angle of attack αs[20] and for plunging and pitching airfoils with several values of h0,α0,αs,k, and ϕ, and the pitch axis located at quarter chord (a−1∕2) [21]. The numerical code and results were validated with experimental data by Anderson et al. [29] for the case of a uniform flow. The main aim of the numerical study by Lian [21] was to investigate the gust effect on the aerodynamic performance of the flapping airfoil, with the objective of finding out the best kinematics conditions to alleviate this effect and ensure a stable flight. Here, we use some of these numerical results, both for a stationary and for a flapping airfoil with the smallest pitch and heave amplitudes, to compare with the present theoretical results for CLand CT. No comparison to Navier–Stokes calculations of Lian and Shy corresponding to the highest amplitudes reported by these authors are shown because in these cases the leading and/or trailing edge separation, not considered in the present potential flow method, becomes very important and, consequently, the theory ceases to be valid and the agreement becomes poorer. Figure 4 shows the comparison of the present theoretical expressions for CLand CTwith Lian and Shyy’s numerical results for a stationary NACA0012 airfoil [20], with and angle of attack αs4° in a pulsating stream with σ0.2, and k10.085.U2tis included with a dashed red line in this figure and in the following ones because it is a relevant reference for two reasons: it relates the plotted results to the velocity magnitude of the pulsating oncoming flow, and it corresponds to the quasi-static lift coefficient normalized with its timeaveraged value for a given αs. In the case of CT, we subtract a quasistatic drag CD00.05 from the theoretical expression, as previously done in some related comparisons of potential flow theory results with numerical and experimental data [30,31]. This value is selected from the numerical results by Senturk and Smits [32] for a Reynolds number about 4×104. The same drag offset is used in the theoretical CTin all the subsequent comparisons with the numerical results by Lian [21] for a flapping foil with the same Reynolds number. It is observed in Fig. 4 that the agreement between theoretical and numerical results is quite good for both CLtand CTt(the largest relative errors are about 10% for CLtand 15% for CTt,and the average relative errors are approximately 5% for CLand 7% for CT). Figure 5 compares the present results for CLtand CTtwith numerical ones by Lian [21] for a pitching and heaving NACA0012 airfoil with h00.5,α08.42°, ϕ90°, k0.25, and αs8°, in a pulsating stream with σ0.2,k1k∕10, and ϕ1−90° [notice that Lian uses sink1tin the head-on gust instead of the cosk1tϕ1in Eq. (1)]. One cycle of the pulsating stream is shown, comprising 10 flapping cycles since kis 10 times k1. There is a good agreement between the patterns of both sets of results, showing the same undulatory pattern following U2t(also shown in the figure). This pattern is quite different from the purely sinusoidal force histories in a uniform stream, also shown in Fig. 5 by setting σ0in Eqs. (40) and (42) with the same pitching and heaving parameters. The instantaneous largest relative errors are not small [about 58% for CLtand 98% for CTt], but the average relative errors over the gust cycle are much smaller (about 17% for CLand 35% for CT). To better appreciate this qualitative difference of the force histories for pulsating and uniform streams, Fig. 6 shows the corresponding timeaveraged histories of the temporal signals plotted in Fig. 5 by computing the time averages every two flapping cycles, which are plotted in the figure as five connected symbols for each case. For a 6Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 uniform stream (not shown), all the symbols would lie on a horizontal straight line at y1(note that the results are normalized with the time-averaged over the complete gust cycle). The agreement is remarkable between the numerical and theoretical time-averaged lift histories [largest relative errors of about 6% for CLtand 22% for CTt, and average relative errors of approximately 4% for CLand 14% for CT]. The agreement is apparently much better than in the instantaneous results plotted in Fig. 5, probably an effect of dividing by the overall time-averaged lift and thrust in each case. But this normalization highlights the effect of the fluctuating velocity on the aerodynamic forces, as done in the results plotted in Lian’s Fig. 8 with which we compare. Also noticeable is the strong effect of a pulsating stream with an amplitude of just σ0.2 on both thrust and lift over a gust cycle, both following the freestream velocity pattern in this case. 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Present Lian ([21]) =0 U2 a) 0 0.2 0.4 0.6 0.8 1 -0.5 0 0.5 1 1.5 Present Lian ([21]) =0 U2 b) Fig. 5 CLt(a) and CTt(b) from Eqs. (40) and (42) compared with numerical results from Fig. 9 of Lian [21] for a NACA0012 airfoil with h00.5, α08.42°, ϕ90°, k0.25, and αs8° in a stream with σ0.2,k1k∕10, and ϕ1−90°. Also plotted results for a uniform stream (σ0) and U2 for reference sake. 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 1.2 1.4 1.6 Present Lian and Shyy ([20]) U2 a) 0 0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1 1.2 1.4 1.6 Present with CD0=0.05 Lian and Shyy ([20]) U2 b) Fig. 4 Normalized CLt(a) and −CTt(b) from Eqs. (47) and (49) compared with numerical results from Fig. 20 of Lian and Shyy [20] for a stationary NACA0012 airfoil with αs4°, σ0.2, and k10.085. 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 Present Lian ([21]) U2 a) 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 Present Lian ([21]) U2 b) Fig. 6 Normalized time-averaged force histories corresponding to the case plotted in Fig. 5. The averages are computed every two flapping cycles, with circles corresponding to the present theoretical results and squares to Lian’s numerical results (from Fig. 8 in Ref. [21]). Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO 7 Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 However, this is not always so, as shown in Fig. 7 for a different flapping configuration [h00.75,α08.42°, ϕ75°, k0.63, and αs0] in the same pulsating stream. In this case, the pattern of the time-averaged thrust history shown in Fig. 7b is opposite to that of U2, with a smaller amplitude of the time-averaged thrust during the gust cycle than in the previous case, as it is more evident in the thrust histories shown in Fig. 7a, where the case with σ0is also included. The time-averaged thrust shows less variation along the gust cycle compared to the case considered in Fig. 6, so that the gust effect of the pulsating stream seems to be alleviated with this particular flapping kinematics. To show that the results remain accurate even for higher amplitudes and frequencies of the pulsating stream, Fig. 8 compares timeaveraged lift and thrust histories for a given flapping kinematics in a pulsating stream with two amplitudes, σ0.2 and 0.5, and a lower frequency, k1k∕20. The time averages are now computed every four cycles, so that five points per gust cycle are shown in Fig. 8 in each case, as in previous figures with twice k1∕k. It is remarkable that there is excellent agreement between theoretical and numerical lift, which practically coincide for both values of σ(for both σ0.2 and σ0.5, the largest relative errors are about 2% and the mean relative errors about 1%), undoubtedly due to the smaller heaving amplitude in this case (h00.25), thus more according with the present linear theory. Obviously, the variation along the gust cycle of the timeaveraged force histories increases with the gust amplitude σ, but there are no significant differences in the comparison between theoretical and numerical results, with a much better agreement for the lift than for the thrust (the largest relative error for CTis about 24% for σ 0.2 and 53% for σ0.5, while the mean relative errors are approximately 16 and 39%, respectively). This is probably due to viscous effects that have a larger influence on the forces in the streamwise direction (due to the viscous drag) compared to the vertical (lift) direction. But in both cases, theoretical and numerical results follow the same pattern during the gust cycle. For σ0.5, the peaks of the time-averaged lift and thrust histories during the gust cycle are about twice their corresponding meanvalues (notice that the results in Fig. 8 are normalized with the corresponding mean values along the complete gust cycle). VI. Conclusions General expressions for the lift, thrust, and moment of a pitching and heaving airfoil immersed in a pulsating freestream are derived from linear potential flow theory using a vortical impulse formulation. The results coincide with previous ones by Greenberg [15] for the lift and moment, derived from a more standard linearized potential flow theory, except for a new circulatory term missing in Greenberg’s results, which we show is quantitatively unimportant. What is more relevant is that we add here general expressions for the instantaneous and time-averaged thrust force and the corresponding propulsive efficiency. The theoretical results are also compared with previous ones for a stationary airfoil with a given angle of attack, and validated with available numerical results by Lian and Shyy [20,21] for Reynolds number 4×104, several pitching and heaving kinematics and different amplitudes and frequencies of the pulsating stream, with pulsating flow frequencies 10 or more times larger than the flapping frequency. A reasonably good agreement is found with these numerical results, especially for the smallest values of the pitch and heave amplitudes considered by these authors, in accordance with the 0 0.2 0.4 0.6 0.8 1 -2 0 2 4 6Present Lian ([21]) = 0 U2 a) 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 Present Lian ([21]) U2 b) Fig. 7 Comparison between the normalized CTt(a) and time-averaged CT(b) with numerical results from Fig. 11 of Lian [21] for a NACA0012 airfoil with h00.75,α08.42°, ϕ75°, k0.63, and αs0in a stream with σ0.2,k1k∕10, and ϕ1−90°. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 Present =0.2 Lian ([21]) =0.2 Present =0.5 Lian ([21]) =0.5 a) 0 0.2 0.4 0.6 0.8 1 -0.5 0 0.5 1 1.5 2 2.5 Present =0.2 Lian ([21]) =0.2 Present =0.5 Lian ([21]) =0.5 b) Fig. 8 Comparison between normalized time-averaged histories of CL(a) and CT(b) (circles) with numerical results from Fig. 8 of Lian [21] (squares) for a NACA0012 airfoil with h00.25,α08.42°, ϕ90°, k0.25,andαs8°inastreamwithk1k∕20,ϕ1−90°, and two values of σ, as indicated. 8Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925 linearized character of the present theory. The agreement is much better for the lift than for the thrust force, particularly for the timeaveraged history. This is a common feature with previous linearized potential flow theories, for the viscous effects, not considered in these theories, affects the thrust more than the lift. To palliate this, a constant offset representing the viscous static drag is usually subtracted from the theoretical thrust to present a more meaningful comparison with measured or numerically computed data, as done here, but which never substitutes accurately the actual, time-varying viscous drag [9]. It is remarkable that the agreement with the numerical results remains good even for amplitudes of the pulsating stream as large as 50% of the mean freestream velocity, where peaks of the time-averaged lift and thrust over a flapping cycle may become about twice the time-averaged lift and thrust over the complete pulsating stream cycle. These strong unsteady effects in the aerodynamic forces associated to the nonstationary freestream velocity, or, equivalently, to the nonstationary cruising velocity of a self-propelled foil, may have a large impact on the propulsive performance of a flapping foil. Hence the convenience of using the expressions derived here, instead of the more conventional ones obtained for a constant stream velocity, to model the self-propulsion by a small-amplitude flapping foil. Appendix: Time-Averaged Coefficients for Special Cases The time-averaged coefficients for some special (integer) combinations between kand k1for which the general expressions in Sec. IV are not valid are summarized here. Only the coefficients that differ from those in Sec. IV are written below. A. kk1  CL2παsσ2Fk 2πσα0cosϕ−ϕ1πkσh0Gkcosϕ1 −Fksinϕ1πσα0Fkcosϕ−ϕ1−Gksinϕ−ϕ1 πkσα0Gkcosϕ−ϕ1Fksinϕ−ϕ1a−1 2 (A1)  CT CT0−σαs2πFk2G1k4kσαsh0G1ksinϕ1 −σα02 2π 2F2kG12k−πkσαsα0 2sinϕ1−ϕ −πσαsα0Fkcosϕ1−ϕGksinϕ1−ϕ −πkσαsα0Gkcosϕ1−ϕ−Fksinϕ1−ϕa−1 2 2σαsα0F1ksinϕ1−ϕ−G1kcosϕ1−ϕ kσαsα0F1kcosϕ1−ϕ4a−3G1ksinϕ1−ϕ (A2) B. k2k1  CL2παsσ2Fk 2πσ2α0 2Fk1cosϕ−2ϕ1 −Gk1sinϕ−2ϕ1 (A3)  CT CT0−σαs22G1k1πFk1 −σα02 2π 2F3k1π 2Fk1G13k1G1k1 −2σ2αsα0π 2Fk1Gk1cosϕ−2ϕ1(A4) C. 2kk1  CT CT0−σα02 2π 2F3k1π 2Fk1G13k1G1k1 −σαs22G1k1πFk1 kσh0α02F1k−π 2Gkcosϕ1−ϕ −π 2Fk2G1ksinϕ1−ϕ σα2 0πk 43π 2Gk−3F1k−π 2kFka−1 2 −kG1k2a−3 2sinϕ1−2ϕ −3π 2Fk−3G1k−π 2kGka−1 2 kF1k2a−3 2cosϕ1−2ϕ(A5) Acknowledgments The authors acknowledge support from the Advanced Grant of the European Research Council GRIFFIN, Action 788247, and from the Junta de Andalucía, Spain, Grant UMA18-FEDER-JA-047. Ernesto Sanchez-Laulhe also acknowledges his predoctoral contract at the University of Malaga. References [1] Lauder, G. V., Anderson, E. J., Tangorra, J., and Madden, P. G. A., “Fish Biorobotics: Kinematics and Hydrodynamics of Self-Propul- sion,”Journal of Experimental Biology, Vol. 210, No. 16, 2007, pp. 2767–2780. https://doi.org/10.1242/jeb.000265 [2] Jones, K. D., and Platzer, M. F., “Design and Development Considerations for Biologically Inspired Flapping-Wing Micro Air Vehicles,” Experiments in Fluids, Vol. 46, May 2009, pp. 799–810. https://doi.org/10.1007/s00348-009-0654-1 [3] Wen, L., Wang, T. M., Wu, G. H., and Liang, J. H., “Hydrodynamic Investigation of a Self-Propelled Robotic Fish Based on a Force- Feedback Control Method,”Bioinspiration & Biomimetics, Vol. 7, No. 3, 2012, Paper 032012. https://doi.org/10.1088/1748-3182/7/3/036012 [4] Gerdes, J., Holness, A., Perez-Rosado, A., Roberts, L., Greisinger, A., Barnett, E., Kempny, J., Lingam, D., Yeh, C. H., Bruck, H. A., and Gupta, S. K., “Robo Raven: A Flapping-Wing Air Vehicle with Highly Compliant and Independently Controlled Wings,”Soft Robotics, Vol. 1, No. 4, 2014, pp. 275–288. https://doi.org/10.1089/soro.2014.0019 [5] Folkertsma, G. A., Straatman, W., Nijenhuis, N., Venner, C. H., and Stramigioli, S., “Robird: A Robotic Bird of Prey,”IEEE Robotics & Automation Magazine, Vol. 24, No. 3, 2017, pp. 22–29. https://doi.org/10.1109/MRA.2016.2636368 [6] Zhu, J., White, C., Wainwright, D. K., Di Santo, V., Lauder, G. V., and Bart-Smith, H., “Tuna Robotics: A High-Frequency Experimental Platform Exploring the Performance Space of Swimming Fishes,”Science Robotics, Vol. 4, No. 34, 2019, Paper eaax4615. https://doi.org/10.1126/scirobotics.aax4615 [7] Alben, S., and Shelley, M., “Coherent Locomotion as an Attracting State for a Free Flapping Body,”Proceedings of the National Academy of Sciences, Vol. 102, No. 32, 2005, pp. 11,163–11,166. https://doi.org/10.1073/pnas.0505064102 [8] Sánchez-Rodríguez, J., Raufaste, C., and Argentina, M., “A Minimal Model of Self Propelled Locomotion,”Journal of Fluids and Structures, Vol. 97, Aug. 2020, Paper 103071. https://doi.org/10.1016/j.jfluidstructs.2020.103071 [9] Fernandez-Feria, R., Sanmiguel-Rojas, E., and Lopez-Tello, P. E., “Numerical Validation of Simple Non-Stationary Models for Self- Propelled Pitching Foils,”Ocean Engineering, Vol. 260, Sept. 2022, Article in Advance / SANCHEZ-LAULHE, FERNANDEZ-FERIA, AND OLLERO 9 Downloaded by 79.147.122.31 on June 26, 2023 | http://arc.aiaa.org | DOI: 10.2514/1.J062925