2.1. Theoretical model

We consider a thin membrane aerofoil that is held by supports at a distance $c=2b$ from the leading to the trailing edge. The membrane is immersed in a uniform inviscid incompressible flow of fluid density $\rho$ and variable speed $U(t)$ at an angle of attack $\alpha$ (see figure 2b), as encountered by an aerofoil section along the span of a hovering membrane wing (figures 1b and 2a). The unsteady lift acting on the membrane aerofoil is derived here assuming small membrane deformation $(|y_x|\ll 1)$ and attached flow at high angles of attack. These assumptions are justified by the experimental results of Gehrke & Mulleners (Reference Gehrke and Mulleners2025), in which unsteady membrane deformation was shown to suppress the formation of leading-edge vortices and flow separation, even at angles of attack as high as $55^\circ$ for large parts of the flapping cycle (see figures 1d and 1f).

Figure 2. Representation of (a) the experimental membrane wing section and surrounding flow field and (b) the theoretical model used to derive the unsteady lift coefficient acting on a flexible membrane aerofoil in variable free-stream velocity. The PIV results shown on the experimental model support the assumptions of attached flow and wake extension along the chord line, utilised in the theoretical model.

Under the above assumptions, the normal velocity along the aerofoil due to membrane deformation (normalised by $U(t)$ ) becomes

(2.1) \begin{align} w_{a_d}(\eta ,x) = -\cos {\alpha }\, y_x (\eta ,x) - y_\eta (\eta ,x) , \qquad x\in \left (-1,1\right )\!, \end{align}

where $x$ is a coordinate along the chord, $y$ denotes the membrane profile, $\eta =\int _0^t V(t')\,\mathrm{d} t'$ is the distance travelled by the aerofoil in terms of semi-chord lengths, $V(t)=U(t)/\overline {U}$ is the non-dimensional free-stream velocity and $\overline {U}$ is the mean free-stream velocity. The above change of variables from the standard use of time $t$ to the travelling distance $\eta$ , first introduced by Wu (Reference Wu1971), allows the derivation of closed-form expressions for the unsteady aerodynamic load and the resulting lift coefficient in cases of variable free-stream velocity, as detailed below. Note that for a constant free-stream velocity we obtain $V(t)=1$ and $\eta (t) = t$ , and the generalised solution presented here will converge to the solution of Tiomkin & Jaworski (Reference Tiomkin and Jaworski2022) in cases of small angles of attack.

The fundamental equation of thin aerofoil theory, expressed in terms of $\eta$ , is given by

(2.2) \begin{align} \frac {1}{2\pi }\mathop{\int \kern-10pt -} \nolimits_{-1}^1\frac {\gamma (\eta ,\xi )}{x-\xi }\mathrm{d}\xi = w_{a_d}(\eta ,x) - \frac {1}{2\pi }\int _{1}^{1+\eta }\frac {\gamma _w\left (\eta ,\xi \right )}{x-\xi }\mathrm{d}\xi , \qquad x\in \left (-1,1\right )\!, \end{align}

in which $b,\ b/\overline {U}$ and $U(t)$ are used as the units of length, time and circulation per unit length, respectively. Here $\gamma$ is the normalised bound vorticity per unit length along the profile, the dashed integral denotes the Cauchy principal value and $\gamma _w (\eta ,\xi )$ describes the normalised vorticity per unit length at location $\xi$ along the wake, $\xi \in [1,1+\eta ]$ , at time $t$ when the aerofoil has travelled a total distance of $\eta (t)$ (see figure 2b).

Wake vortices are assumed to be continuously shed from the trailing edge into a flat wake at the instantaneous free-stream velocity with a fixed strength. Namely, the wake vorticity distribution, $\gamma _w (\eta ,\xi )$ , is equivalent to the vorticity at the trailing edge at an aerofoil displacement of $\eta -\xi +1$ :

(2.3) \begin{align} \gamma _w (\eta ,\xi )=\gamma _w (\eta -\xi +1,1)\triangleq \gamma _{_{T\!E}}(\eta -\xi +1). \end{align}

The above assumptions are supported by experimental observations of the wake (e.g. figure 2a), in which the wake extends along the aerofoil chord line in the vicinity of the membrane trailing edge. By following the standard application of Söhngen’s inversion formula to (2.2) and enforcing Kelvin’s theorem (cf. Söhngen Reference Söhngen1939; Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman1996, p. 289) we obtain

(2.4) \begin{align} 2\int _{-1}^{1} \sqrt {\frac {1+\xi }{1-\xi }}\,w_{a_d}(\eta ,\xi )\,\mathrm{d}\xi = -\int _{0}^{\eta }\sqrt {\frac {\zeta +2}{\zeta }}\gamma _{_{T\!E}} (\eta -\zeta )\, \mathrm{d}\zeta . \end{align}

Tiomkin & Raveh (Reference Tiomkin and Raveh2017) showed that for a constant free-stream velocity, for which $V(t)\equiv 1$ and $\eta \equiv t$ , application of the standard Laplace transform to (2.4) yields a closed-form expression for the wake vorticity distribution in the Laplace plane. In the current case, we utilise the Laplace transform suggested by Wu (Reference Wu1971) to derive a similar closed-form expression for the wake vorticity in the variable velocity case, where the Laplace transform is defined as

(2.5) \begin{align} \bar {F}(s) = \mathcal{L} \left \{F(\eta );s\right \} = \int _0^\infty F(\eta ) \mathrm{e}^{-s\eta }\mathrm{d}\eta . \end{align}

The solution is obtained by utilising a Fourier cosine series expansion to describe the membrane slope in terms of the aerofoil displacement, $\eta$ :

(2.6) \begin{align} y_x(\eta ,\theta )=\frac {1}{2}\mathcal{F}_0(\eta )+\sum _{n=1}^\infty \mathcal{F}_n(\eta )\cos n\theta , \end{align}

where we apply the standard coordinate transformation, $x=-\cos \theta$ , which places the aerofoil leading edge at $x=-1$ ( $\theta =0$ ) and the trailing edge at $x=1$ ( $\theta =\pi$ ). Note that the odd-index Fourier coefficients denote here the amplitude of the corresponding membrane mode-shape in the decomposition of the membrane deformation, while a combination of the even-index Fourier coefficients provides the even mode-shapes (see Tiomkin & Raveh Reference Tiomkin and Raveh2017). By substituting (2.6) into the normal velocity expression (2.1) and applying the Laplace transform to (2.4), we obtain an expression for the wake vorticity distribution in the Laplace domain:

(2.7) \begin{align} \bar {\gamma }_{T\!E}(s)= -\frac {2\pi \mathrm{e}^{-s}}{K_0(s)+K_1(s)} \bar {f}_m(s) = -2\pi \, \bar {\varPsi }(s) s^2 \bar {f}_m(s) , \end{align}

where $\bar {\varPsi }(s)$ is the Laplace transform of Küssner’s function (e.g. Sears Reference Sears1940), $K_0$ and $K_1$ are modified Bessel functions of the second kind and $\bar {f}_m(s)$ is an auxiliary function in the Laplace domain given by

(2.8)

The aerodynamic load along the aerofoil is then expressed by adapting the method of Schwarz (Reference Schwarz1940) to the variable free-stream velocity formulation, yielding

(2.9) \begin{eqnarray} \frac {\Delta C_p(\eta ,x)}{V(t)^2}&=& -\frac {4}{\pi }\sqrt {\frac {1-x}{1+x}} \mathop{\int \kern-10pt -}\nolimits_{-1}^{1} \sqrt {\frac {1+\xi }{1-\xi }}\,\frac {w_{a_d}(\eta ,\xi )}{x-\xi }\mathrm{d}\xi + \frac {4}{\pi } \mathop{\int \kern-10pt -}\nolimits_{-1}^{1}\varLambda _1 (x,\xi ) \frac {\partial w_{a_d}}{\partial \eta }\mathrm{d}\xi \nonumber \\&&\mbox{} + \frac {2}{\pi }\sqrt {\frac {1-x}{1+x}} \int _{0}^{\eta } \frac {\gamma _{_{T\!E}} (\eta -\zeta )}{\sqrt {\zeta (\zeta +2)}}\mathrm{d}\zeta , \end{eqnarray}

where $\varLambda _1$ is an auxiliary function given by

(2.10) \begin{align} \varLambda _1\left (x,\xi \right )=\ln \left |\frac {\sqrt {(1-x)(1+\xi )}+\sqrt {(1+x)(1-\xi )}} {\sqrt {(1-x)(1+\xi )}-\sqrt {(1+x)(1-\xi )}}\right | \!. \end{align}

It is convenient to denote the integral terms on the right-hand side of (2.9) by $\widetilde {\Delta C}_{\!p_0}$ , $\widetilde {\Delta C}_{\!p_1}$ and $\widetilde {\Delta C}_{\!p_2}$ , which, respectively, represent the quasi-steady, apparent mass and wake contributions to the aerodynamic load along the aerofoil. Closed-form expressions are readily obtained for these terms by substituting (2.1) into the $\widetilde {\Delta C}_{\!p_0}$ and $\widetilde {\Delta C}_{\!p_1}$ terms, while the wake term, $\widetilde {\Delta C}_{\!p_2}$ , is most conveniently expressed in the Laplace plane by applying the convolution theorem and utilising the wake vorticity expression given in the Laplace domain (2.7). Inverse Laplace transform is then required to obtain the total aerodynamic load in the time (or $\eta$ ) domain.

For non-small angles of attack, as considered here, integration on the aerodynamic load (2.9) along the membrane chord line leads to the normal aerodynamic force coefficient due to membrane deformation:

(2.11) \begin{align} C_{n_d}(t) = 2\pi V(t)^2 \left \{ \int _0^{\eta (t)} \varPhi (\eta (t)-\tau )\dot {f}_m(\tau )\,\mathrm{d}\tau +g_m(\eta (t)) \right \} \!, \end{align}

where $\varPhi (t)$ is the time-domain Wagner function, an overdot represents a derivative with respect to $\eta$ and $f_m(\eta )$ and $g_m(\eta )$ are auxiliary functions given by

(2.12)

(2.13)

where $f_m(\eta )$ is the $\eta$ -domain representation of (2.8). Note that $g_m(\eta )$ represents the normalised non-circulatory lift coefficient, while $f_m(\eta )$ describes the circulatory term, as it is directly related to the wake vorticity through (2.7). For a constant free stream at small angles of attack $(\cos {\alpha }\cong 1, \eta \equiv t)$ the solution described above in (2.11)–(2.13) coincides with the unsteady lift coefficient given by Tiomkin & Jaworski (Reference Tiomkin and Jaworski2022).

In cases of harmonic aerofoil motions of reduced frequency $k$ , the normal force coefficient in (2.11) becomes

(2.14) \begin{align} C_{n_d}(t) = 2\pi V(t)^2\left \{ C(k){f}_m\left (\eta \left (t\right )\right ) +g_m\left (\eta \left (t\right )\right ) \right \}\! , \end{align}

and the unsteady lift coefficient due to membrane deformation is

(2.15) \begin{align} C_{l_d}(t) = C_{n_d}\cos {\alpha } = 2\pi \cos {\alpha }\, V^2(t)\left \{ C(k){f}_m\left (\eta \left (t\right )\right ) +g_m\left (\eta \left (t\right )\right ) \right \}\! , \end{align}

where $C(k)$ is the frequency-domain Theodorsen’s function:

(2.16) \begin{align} C(k) = \frac {H_1^{(2)}(k)}{H_1^{(2)}(k)+\mathrm{i} H_0^{(2)}(k)}, \end{align}

in which $H_0^{(2)}$ and $H_1^{(2)}$ are Hankel functions of the second kind.

We further note that in cases where the first mode is significantly more dominant than all other modes $(\mathcal{F}_1\gg \mathcal{F}_n\; \forall \, n\neq 1)$ , the membrane deformation can be represented by $y(\eta ,\theta )=({1}/{2})\mathcal{F}_1(\eta ) \sin ^2{\theta }$ , where $({1}/{2})\mathcal{F}_1(\eta )$ is the maximum camber of the aerofoil. Under these conditions a simplified expression is readily obtained for the lift coefficient that depends only on the wing kinematics and the first Fourier coefficient (or maximum membrane camber):

(2.17)

The above expressions (2.11)–(2.15) were derived assuming small membrane camber, $|y_x|\ll 1$ , constant non-small angle of attack and attached flow. For scenarios in which these assumptions apply, the unsteady lift coefficient produced by the membrane deformation is given by (2.15) and is determined by the flow parameters ( $\alpha$ , $V(t)$ and $k$ ) and the Fourier coefficients that decompose the membrane deformation in time. While the flow parameters (or wing kinematics) are generally given for every scenario considered, the Fourier coefficients can be obtained either through a numerical solution of the coupled membrane–fluid problem or by using membrane deformation measurements to predict the lift.

Figure 3. (a) Sketch of the flexible membrane wing profile. (b) Angle of attack ( $\alpha$ ) and local wing velocity ( $U$ ) profiles as a function of normalised time $t/T$ . (c) Membrane camberline, $y(x)$ , obtained at different times ( $t/T = 0$ – $0.5$ ). (d) Membrane Fourier coefficients, $\mathcal{F}_n$ , as a function of time $t/T$ . (e) Measured lift coefficients $C_L$ as a function of time $t/T$ for a membrane and a rigid wing. (f) Finite-wing lift coefficient due to membrane deformation $C_{L_d}$ as a function of time $t/T$ ; comparing the experimental measurements with the theory prediction, along with simplified model predictions assuming single-mode oscillations ( $\mathcal{F}_1$ ), constant free-stream velocity ( $\overline {U}$ ) and a quasi-steady aerodynamic model (QS). Very good agreement is observed between the theoretical and measured lift coefficients, computed for $\hat {\alpha }=35^\circ , Ae=2.5$ , which is also captured by the single-mode approximation. However, the simplified constant- $\overline {U}$ model and the quasi-steady model do not agree well with measurements.

In the current work we choose the latter approach and apply the Fourier series expansion (2.6) to experimental data of the membrane deformation in time, from which time-dependent Fourier coefficients are computed (figure 3). These Fourier coefficients, $\mathcal{F}_n(\eta )$ , when substituted in (2.12), (2.13) and (2.15) provide a theoretical prediction of the sectional lift coefficient due to membrane deformation in the experimental conditions. Therefore, this method loosens the requirements of a constant tension in the membrane and small angles of attack, as commonly assumed in theoretical studies, while still limiting the validity of the results to scenarios of small membrane deformations and attached flow.

2.2. Experimental configuration

The experimental set-up of Gehrke et al. (Reference Gehrke, Richeux, Uksul and Mulleners2022) is briefly described here in the context of the current work, as this work serves for validation of the novel theoretical model presented in § 2.1. The experimental configuration features a membrane wing that is attached to a rigid frame that performs a prescribed flapping motion in a water tank. The membrane wing includes rigid leading and trailing edges that are free to rotate, allowing the wing edges to align favourably with the flow during each flapping cycle. An additional degree of freedom is maintained at the trailing edge by permitting free translation of the edge in the chordwise direction, varying the membrane chord length in time.

This configuration enables passive deformation of the thin aerofoil shape under aerodynamic loading, creating a natural synergy between structural flexibility and aerodynamic forces, characterised by the aeroelastic number $ \textit{Ae} = \textit{Eh}/(({{1}/{2})\rho \overline {U}^2 c})$ , where $E$ is the Young modulus of the membrane. The model wing performs a sinusoidal stroke along the horizontal plane, and a trapezoidal angle-of-attack ( $\alpha$ ) profile, similar to the hovering kinematics of small bats and insects (see figures 1b and 3b showing the experimental set-up and wing kinematics). The experimental study focused on capturing the membrane dynamic response and the unsteady aerodynamic forces produced by the flapping wing for a wide range of aeroelastic numbers $ \textit{Ae} = 0.25{-}12$ and angle-of-attack amplitudes of $\hat {\alpha }=15^\circ {-}75^\circ$ . Aerodynamic forces and moments were recorded at the wing root with a six-axis force/torque transducer, while the wing’s deformation was captured throughout the flapping cycle with two machine vision cameras in stereo configuration and reconstructed using photogrammetry. Three-dimensional deformation results showed no significant spanwise variation in membrane profile, allowing its representation by a single two-dimensional profile (Gehrke & Mulleners Reference Gehrke and Mulleners2025). Finally, the unsteady flow field around the membrane wing was recorded using planar particle image velocimetry (PIV) at a fixed spanwise position, located at the radius of the second moment of area ( $r = 0.55R$ , where $R$ is the total wingspan). More information about the experimental set-up, data acquisition and processing and in-depth aerodynamic analysis can be found in Gehrke et al. (Reference Gehrke, Richeux, Uksul and Mulleners2022) and Gehrke & Mulleners (Reference Gehrke and Mulleners2025).

2.3. Connecting theory to measurements

Seeking to validate the unsteady aerodynamic theory derived in § 2.1, membrane deformation measurements obtained by Gehrke et al. (Reference Gehrke, Richeux, Uksul and Mulleners2022) along the mid-wingspan section are converted first to body-fixed coordinates, normalised by the instantaneous semi-chord length, and decomposed into a Fourier series expansion (2.6) (see figures 3c and 3d, respectively). The resultant Fourier coefficients are then used to compute the theoretical sectional lift coefficient, $C_{l_d}(t)$ , via (2.15). This result yields the finite-wing lift contribution from membrane deformation:

(2.18) \begin{align} C_{L_d}(t) = \frac {C_{L_\alpha }}{2\pi } \frac {c(t)}{\overline {c}} C_{l_d}(t) , \end{align}

where $c(t)$ and $\overline {c}$ are the instantaneous and reference chord length, respectively. The finite-wing lift slope $C_{L_\alpha }$ is obtained using the two-station approximation (e.g. Anderson Reference Anderson2024, pp. 451–454) for a rectangular wing of aspect ratio , matching the wing geometry used in the experiments of Gehrke et al. (Reference Gehrke, Richeux, Uksul and Mulleners2022).

The corresponding experimental value for the lift coefficient due to membrane deformation is obtained by calculating the difference between the lift coefficient measured for the membrane wing and the lift coefficient obtained on a rigid flat plate that undergoes an equivalent prescribed flapping motion (see figure 3e). The resultant experimental lift coefficient due to membrane deformation is then compared with the corresponding theoretical value (2.18) for every experimental scenario considered (e.g. figure 3f). Note that for cases of dominant first-mode oscillations the sectional lift coefficient in (2.18) can be obtained from the approximate expression (2.17). In the scenario presented in figure 3, the first Fourier coefficient is clearly the most dominant coefficient and therefore the lift coefficient predicted with the single-mode approximation is practically indistinguishable from the full-model prediction. However, for more general cases, in which higher-mode oscillations are observable in the membrane dynamic response, the full model (2.15) is necessary to predict the sectional lift coefficient. Such higher-mode oscillations are expected, for example, in cases of higher reduced frequency, beyond the first natural frequency of the membrane. An alternative use of the standard unsteady aerodynamic model with constant free-stream velocity significantly under-predicts the lift coefficient, while a quasi-steady model strongly over-predicts the response (figure 3f), highlighting the importance of including the effect of the variable free-stream velocity and the shedding wake in the predictive model. Note that the grey areas in figure 3(b,d–f) mark the time intervals during which the angle of attack varies in time in the experiments, violating the theory assumptions. Therefore, the theoretical prediction, in which only effects of the instantaneous angle of attack are considered, is expected to agree with the experimental results for $0.08\lt t/T\lt 0.42$ , with a less favourable agreement in the grey-shaded areas.