1. Introduction
The growing industrial interest in small-scale unmanned aerial vehicles (SUAVs) for sensing, reconnaissance and parcel delivery continues to spur scientific interest into novel aerodynamic design solutions for low-speed flows, inspired by biological fliers (Hassanalian & Abdelkefi Reference Hassanalian and Abdelkefi2017). A special focus has been given to the membrane wings of bats, who possess impressive maneuvering and gliding abilities without relying upon high-frequency wing flapping for lift (Hedenström & Johansson Reference Hedenström and Johansson2015). These membrane wings are lightweight and are therefore appealing for SUAV applications. However, the compliance of membrane wings couples their geometrical shape and dynamics to the surrounding fluid mechanics and complicates the prediction of their aerodynamic performance. Several computational and experimental studies have examined the aerodynamics of these flexible membrane wings in steady flow conditions over the last two decades (e.g. Song et al. Reference Song, Tian, Israeli, Galvao, Bishop, Swartz and Breuer2008; Gordnier Reference Gordnier2009; Rojratsirikul, Wang & Gursul Reference Rojratsirikul, Wang and Gursul2009; Arbos-Torrent, Ganapathisubramani & Palacios Reference Arbos-Torrent, Ganapathisubramani and Palacios2013; Serrano-Galiano, Sandham & Sandberg Reference Serrano-Galiano, Sandham and Sandberg2018). Gordnier (Reference Gordnier2009) and Rojratsirikul et al. (Reference Rojratsirikul, Wang and Gursul2009) showed that membrane wings in low-Reynolds-number flows delay stall and enhance the mean lift coefficient due to the onset of flow–membrane oscillations. The oscillations of the compliant membrane are essential to these aerodynamic benefits, as Gordnier (Reference Gordnier2009) found no significant advantage for a static membrane wing when compared with an equivalent (cambered) rigid aerofoil. Thus, the unsteady behaviour of membrane wings is a principal source of interest to discover mechanisms for aerodynamic performance improvement. Although many studies investigated the membrane wing response to steady flow conditions (see Tiomkin & Raveh Reference Tiomkin and Raveh2021 for a recent literature survey), few have studied its response to unsteady flow conditions or prescribed motions, where aeroelastic membrane deformation may yield further aerodynamic benefits.
The pursuit of novel mechanisms to improve the unsteady aerodynamic performance of membrane wings and understand their associated fluid mechanics has led to a recent research focus on flapping membrane wings. The combination of aerofoil flexibility with flapping motions can potentially eliminate flow separation along the aerofoil altogether and improve the aerodynamic maneuverability of the vehicle, as demonstrated in the context of bat flight (Muijres et al. Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström2008; Chin & Lentink Reference Chin and Lentink2016). Several computational (Gopalakrishnan & Tafti Reference Gopalakrishnan and Tafti2010; Jaworski & Gordnier Reference Jaworski and Gordnier2015) and experimental (Tregidgo, Wang & Gursul Reference Tregidgo, Wang and Gursul2013) studies investigated the membrane wing response to prescribed flapping motions, which are generally described as pitch or heave oscillations, or as a combination of the two. Gopalakrishnan & Tafti (Reference Gopalakrishnan and Tafti2010) used coupled large eddy simulations of a rectangular membrane to show that induced camber enhances both the lift and the thrust during a flapping pitching motions. These simulations identified the movement of the leading-edge vortex along the membrane aerofoil surface to be the main source of the increased lift and thrust relative to a flapping rigid wing, for which the leading-edge vortex detaches and moves away of the wing which causes a drop in the lift coefficient (Eldredge & Jones Reference Eldredge and Jones2019). These results are supported by the computations of Jaworski & Gordnier (Reference Jaworski and Gordnier2015) that focused on the role of prestress and elastic modulus in the propulsion of a flapping membrane aerofoil, which is enhanced by the interaction of the leading-edge vortex with the local elastic deformation. Experiments of Tregidgo et al. (Reference Tregidgo, Wang and Gursul2013) focused on the membrane dynamic response to a transient sinusoidal pitch maneuver of reduced frequency of
$k=0.022$ and amplitude of
$10^\circ$. Different vibrational modes were identified that depended on the stationary angle of attack about which the unsteady maneuver was carried out. For small stationary angles of attack
$(0^\circ \le \alpha \le 4^\circ )$, first mode oscillations were observed with a small lag in the membrane dynamic response relative to the prescribed motion. This delay was more pronounced for a larger stationary angle of attack of
$\alpha =10^\circ$, for which hysteresis was identified between the pitch-up and the pitch-down sections of the motion, which were accompanied by different vibrational modes.
The above studies collectively emphasise the complexity of the flapping membrane wing problem. However, due to their computational and experimental nature, their scope is limited to a few specific points in the parameter regime of flapping membrane wings, i.e. specific values of reduced frequency, mass ratio and membrane elasticity. An analytical solution of a simplified model problem is therefore desired to shed light on the role of each dimensionless group in the wide parameter space of flapping membrane wings. Physical insights from such an analytical solution are expected to inform future computational and experimental studies en route to obtaining a more complete understanding of the physics of flapping membrane wings.
Several analytical studies have been carried out that focus on propulsive thrust and efficiency predictions for flapping flexible wings. Alben (Reference Alben2008) presented an analytical solution for a flapping inextensible elastic sheet (with a free trailing edge), utilising unsteady thin aerofoil theory coupled to a beam structural model. Their work identified an optimal thrust condition at the resonance peaks for small pitching amplitudes. More recently, Alon Tzezana & Breuer (Reference Alon Tzezana and Breuer2019) showed via analytical solution that a heaving membrane wing transitions between thrust and drag near the membrane resonance frequency, as the reverse von Kármán wake transitions to a traditional von Kármán wake. In their reviews on flapping wing aerodynamics of biological and bio-inspired flyers, Shyy et al. (Reference Shyy, Aono, Kang and Liu2013, Reference Shyy, Kang, Chirarattananon, Ravi and Liu2016) highlighted the importance of using a time-domain approach to predict accurately the aerodynamic performance of flapping wings at the scale of bats and birds due to the inseparable flapping and body timescales of these flyers, which is not the case for smaller insect-scale flyers. Thus, although a quasi-steady model can make accurate predictions for insect-scale vehicles, this model assumption is not recommended for SUAV applications, where a time-dependent approach is essential to address vehicle stability and control. Furthermore, whilst the studies of Alben (Reference Alben2008) and Alon Tzezana & Breuer (Reference Alon Tzezana and Breuer2019) elucidate the propulsive potential of flapping flexible wings, a theoretical basis to understand the membrane wing aerodynamic performance in prescribed flapping motions remains underdeveloped, specifically in terms of the ability to predict its unsteady lift and structural dynamic response.
In addition to the unsteady lift and thrust enhancement mechanisms engendered by membrane wings under prescribed flapping motions, an understanding of the response of these wings to flow disturbances such as gusts is important to the design of membrane wing SUAVs. Due to their small size and slow flight speed, SUAVs are especially susceptible to flight disruption from small gusts typical of urban environments (Watkins et al. Reference Watkins, Milbank, Loxton and Melbourne2006; Elbanhawi et al. Reference Elbanhawi, Mohamed, Clothier, Palmer, Simic and Watkins2017; Jones, Cetiner & Smith Reference Jones, Cetiner and Smith2022). Classical linear unsteady aerodynamic theory (see von Kármán & Sears Reference von Kármán and Sears1938; Sears Reference Sears1940) predicts the transient lift response of a rigid aerofoil to transverse gusts of small gust ratios, where the gust ratio is the transverse gust amplitude divided by the freestream flow speed. This theory has long been utilised to predict the unsteady load on rigid wings in terms of lift amplitude and phase lag. However, when compliant membrane wings are considered, the lift response is composed of both the local change in angle of attack and the resulting deformation of the aerofoil. The membrane deformation couples aeroelastically to the aerodynamic load, which may amplify or attenuate the unsteady lift response. Initial results by Berci et al. (Reference Berci, Gaskell, Hewson and Toropov2013) from a semi-analytical state-space model indicate the appearance of structural oscillations in the massless membrane response to a sharp-edged transverse gust. However, these oscillations were described only in terms of the mid-chord membrane deformation and without consideration of the structural mode of oscillation and the lift response of the aerofoil. A complete analysis of the membrane response to unsteady flow is currently lacking in the literature.
The current study aims to fill this knowledge gap by presenting an unsteady analytical model and its solution for a membrane wing in inviscid incompressible flow, under the unsteady conditions of prescribed motions or transverse gust profiles. The transient membrane response is determined in the Laplace domain; steady-state harmonic oscillations of the membrane deformation and the unsteady lift response are investigated using a simplified solution in the frequency (Fourier) domain, which is convenient to compare against established rigid aerofoil theory. These solutions yield novel extensions to the classical unsteady aerodynamic functions for flexible membrane wings.
The remainder of this paper is organised as follows. Section 2 presents the mathematical problem for the generalised case of a membrane wing in arbitrary motion or gust, and for specific canonical unsteady flow scenarios. In § 3, the results of the theoretical model are presented in terms of membrane wing deformation and aerodynamic performance, as represented by extensions to the classical unsteady aerodynamic functions by Theodorsen, Wagner, Sears and Küssner. Section 4 closes with concluding remarks.
2. Formulation
2.1. Membrane wing
Consider an extensible membrane aerofoil of thickness
$h$ and density
$\rho _m$, which is held by simple supports at a distance
$2b$ from one edge to the other. The membrane is initially still and taut, and is immersed in a uniform and inviscid incompressible freestream of density
$\rho$ and speed
$U$, aligned parallel to the membrane chord (see figure 1a). Assuming small deformations of the membrane,
$|\tilde {y}_{\tilde {x}}|\ll 1$, the membrane dynamic equation is
where
$\tilde {y}$ denotes the membrane profile,
$\tilde {t}$ represents time,
$\tilde {x}$ is a coordinate along the chord and
$T$ and
$\Delta p$ are the tension and pressure difference along the membrane, respectively. While the membrane is extensible, we note that Tiomkin & Raveh (Reference Tiomkin and Raveh2017) showed that the tension can be considered constant to leading order for the small angles of attack and deformations assumed in the current study.
Sketch of the membrane aerofoil gust problem: (a) initial time
$t=0$; (b) later time
$t>0$. Note that the membrane is initially taut, at zero angle of attack, and deforms under transient (gust) fluid loads.

The non-dimensional form of the dynamic equation is
in which
$b$,
$b/U$,
$\rho$,
$U$,
$\frac {1}{2}\rho U^2$ and
$\rho U^2 b$ are used as the units of length, time, density, circulation (per unit length), pressure and force (per unit span). Note that
$b$ is used as the unit of length throughout the dynamic equation, but the mass ratio is normalised with
$c$ as the unit of length, namely
$\mu =\rho _m h/\rho c$, following the convention in previous membrane wing studies (Alon Tzezana & Breuer Reference Alon Tzezana and Breuer2019; Jaworski & Gordnier Reference Jaworski and Gordnier2012). The mass ratio
$\mu$ and tension coefficient
$C_T$ are fixed parameters in the present analysis, and the unsteady membrane deformation and pressure coefficient profiles,
$y$ and
$\Delta C_p$, respectively, are part of the solution. A schematic drawing of the membrane geometry in the non-dimensional form is presented in figure 1 for the gust response case; this coordinate system is used to describe the membrane deformation in all of the considered cases.
2.2. Incompressible potential flow
The extensible membrane aerofoil may encounter or produce an unsteady flow field that superposes on the uniform background flow. Inviscid, incompressible potential flow is considered with an initial angle of attack of
$\alpha =0^\circ$, which isolates the effects of unsteady angle of attack variations or transient gusts on the membrane dynamic response; this approach is similar to the traditional formulation available for the arbitrary motion of rigid aerofoils (e.g. Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman1996). For completeness, this section outlines the formulation of Tiomkin & Raveh (Reference Tiomkin and Raveh2017) for a membrane aerofoil in steady flow and extends it to include the dynamic membrane response to an unsteady flow.
The standard coordinate transformation
places the profile leading edge at
$x=-1$ (
$\theta =0$) and the trailing edge at
$x=1$ (
${\theta ={\rm \pi}}$). This coordinate transformation permits the membrane slope,
$y_x$, to be expressed as a Fourier cosine series expansion per Nielsen (Reference Nielsen1963), which is augmented here by allowing the Fourier coefficients to be time-dependent:
\begin{equation} y_x(t,\theta)=\frac{1}{2}F_0(t)+\sum_{n=1}^\infty F_n(t)\cos n\theta . \end{equation}
We proceed with expressing the membrane dynamic equation (2.2) in terms of the new coordinate
$\theta$, which yields a system of differential equations for the Fourier coefficients.
Integration of (2.4) along the horizontal coordinate, from the leading edge to a point
$x$ along the chord, yields the membrane profile:
\begin{align} y(t,\theta)&= \frac{1}{2}F_0(t)(1-\cos\theta)+\frac{1}{2}F_1(t)\sin^2\theta - \frac{1}{2}\sum_{n=2}^\infty F_n(t) \left(\frac{1}{n^2-1}\right)\nonumber\\ &\quad \times [2+(n-1)\cos (n+1)\theta - (n+1)\cos (n-1)\theta] , \end{align}which has to sustain the fixed boundary conditions of the membrane edges. The leading-edge boundary condition is automatically satisfied by (2.5). However, the fixed trailing-edge boundary condition imposes the constraint
\begin{equation} F_0(t)=2\sum_{\substack{m=1}}^\infty \frac{F_{2m}(t)}{(2m)^2-1} . \end{equation}In addition, the assumption of an initially still and taut membrane yields zero-valued initial conditions for the Fourier coefficients and their first time derivative.
The pressure difference across the membrane in (2.2) is obtained by using the unsteady vortex sheet method (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 274),
where
$\gamma$ is the normalised vortex sheet strength per unit length along the profile. The vorticity distribution along the aerofoil is determined by the fundamental equation of thin aerofoil theory,
where the dashed integral denotes the Cauchy principal value. Here
$w_a(t,x)$ is the normal velocity on the membrane surface (normalised by
$U$), and
$\gamma _w (t,\eta )$ describes the normalised vorticity per unit length at location
$\eta$ along the wake,
$\eta \in (1,\infty )$, at time
$t$. Wake vortices are assumed to be continuously shed from the trailing edge at the freestream velocity into a flat wake and have a fixed strength, which asserts that the wake vorticity distribution,
$\gamma _w (t,\eta )$, is equivalent to the vorticity at the trailing edge at time
$t-\eta +1$:
Application of Söhngen's inversion formula to (2.8) and enforcement of Kelvin's theorem (cf. Söhngen Reference Söhngen1939; Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 289) leads to
\begin{equation} 2\int_{{-}1}^{1} \sqrt{\frac{1+\xi}{1-\xi}} w_a(t,\xi)\,\mathrm{d}\xi ={-}\int_{1}^{1+t}\sqrt{\frac{\eta+1}{\eta-1}}\gamma_{_{T\!E}} (t-\eta+1) \,\mathrm{d}\eta . \end{equation}Tiomkin & Raveh (Reference Tiomkin and Raveh2017) showed that the application of the Laplace transform to (2.10) yields a closed-form expression for the wake vorticity distribution in the Laplace plane.
Provided that a solution for
$\gamma _{_{T\!E}}$ is obtainable in the time domain, the method of Schwarz (Reference Schwarz1940) (see also Iosilevskii Reference Iosilevskii2007) produces a general expression for the pressure difference along the aerofoil:
\begin{align} \Delta C_p(t,x)&={-}\frac{4}{\rm \pi}\sqrt{\frac{1-x}{1+x}}{\int\hskip -1,05em -\,}_{{-}1}^{1} \sqrt{\frac{1+\xi}{1-\xi}}\frac{w_a(t,\xi)}{x-\xi}\,\mathrm{d}\xi + \frac{4}{\rm \pi}{\int\hskip -1,05em -\,}_{{-}1}^{1}\varLambda_1 (x,\xi) w_{a_t}(t,\xi)\,\mathrm{d}\xi\nonumber\\ &\quad + \frac{2}{\rm \pi}\sqrt{\frac{1-x}{1+x}} \int_{1}^{1+t} \frac{\gamma_{_{T\!E}} (t-\eta+1)}{\sqrt{\eta^2-1}}\,\mathrm{d}\eta , \end{align}
where
$\varLambda _1$ is an auxiliary function expressed in (A3) of Appendix A. The first integral term describes the quasi-steady pressure difference, the second term is the apparent mass contribution (non-circulatory term) and the third term describes the contribution of the wake.
The contributions to expression (2.11) for the aerodynamic load along the membrane may be further separated and analysed by describing the normal velocity on the membrane surface as a superposition:
where
$w_{a_{f}}(t,x)$ is the contribution of the unsteady flow (i.e. prescribed chord motion or a traveling gust) to the normal flow velocity on the membrane, and
$w_{a_{d}}(t,x)$ is the respective contribution of the membrane deformation,
Substitution of (2.12) into (2.11) permits a separation of the effect of the membrane deformation,
$w_{a_d}$, from the effect of the unsteady flow,
$w_{a_f}$, on the aerodynamic load, namely
where the subscripts
$d$ and
$f$ denote terms due to membrane deformation and unsteady flow, respectively. Details of the analytical expressions for
$\Delta C_{p_d}$ (in the Laplace plane) are available in Appendix A of Tiomkin & Raveh (Reference Tiomkin and Raveh2017). We develop in § 2.4 the closed-form expressions for
$\Delta C_{p_f}$ that are necessary to complete the description of the aerodynamic load on a membrane undergoing prescribed chord motion or encountering a gust. Note that the membrane Fourier coefficients appear only in the expression for
$\Delta C_{p_d}$, whilst
$\Delta C_{p_f}$ depends only on the prescribed motion or gust.
The next section combines the terms obtained for the membrane deformation and the resulting aerodynamic load to produce a set of equations for the coupled aeroelastic problem for any arbitrary prescribed chord motion or gust.
2.3. Aeroelastic coupling and methods of solution
The coupled aeroelastic equation that describes the membrane response to unsteady flow conditions is obtained by substituting (2.4), (2.5) and (2.14) into (2.2). This procedure yields a matrix equation in which the unknowns are the Fourier coefficients that describe the membrane deformation. The aeroelastic equation is described and solved in the Laplace domain for generalised time-dependent cases, or in the frequency domain for harmonic motions or gusts. The details of these two methods are given next in §§ 2.3.1 and 2.3.2, respectively. Note that the overbar and hat symbols are used throughout to denote variables in the Laplace and frequency domains, respectively.
2.3.1. Laplace domain
The membrane dynamic equation (2.2) is expressed in the Laplace domain by applying the Laplace transform to (2.4), (2.5) and (2.14) and substituting the resulting expressions into the Laplace transform of (2.2). We then multiply the resulting equation by
$\sin \theta$ and use the mathematical relations (A1) and (A2) to construct a matrix system of equations,
where
$\bar {\boldsymbol {F}}$ is the vector of Fourier coefficients
$\bar {F}_n (s)$,
$n=1\ldots N$ and
$N$ is the number of coefficients chosen to represent the membrane-profile slope in (2.4), taken here as
$N=24$ following the numerical convergence studies of Nielsen (Reference Nielsen1963) and Tiomkin & Raveh (Reference Tiomkin and Raveh2017). Here, the overbar denotes the Laplace transform of the variable,
$\bar {F}_n(s) = \mathcal {L} \{F_n(t);s\}$. The matrices
${{\boldsymbol{\mathsf{A}}}}, {{\boldsymbol{\mathsf{B}}}}, {{\boldsymbol{\mathsf{E}}}}$ and the vector
$\boldsymbol {H}$ are obtained by matching the coefficients of the harmonics of
$\sin \theta$ in the dynamic equation (2.2). It is noted here that matrices
${{\boldsymbol{\mathsf{B}}}}$ and
${{\boldsymbol{\mathsf{E}}}}$ and vector
$\boldsymbol {H}$ depend on the Laplace variable
$s$, while
${{\boldsymbol{\mathsf{A}}}}$ is constant. We further note that the matrices
${{\boldsymbol{\mathsf{A}}}}$,
${{\boldsymbol{\mathsf{B}}}}$ and
${{\boldsymbol{\mathsf{E}}}}$ are obtained from the steady flow solution under zero angle of attack, i.e. by applying
$\Delta {C_p} = \Delta C_{p_d}$ to the dynamic equation; these matrices are detailed in Tiomkin & Raveh (Reference Tiomkin and Raveh2017). The effect of the unsteady flow appears only in
$\boldsymbol {H}$, on the right-hand side of the resulting dynamic equation, and is determined by
$\Delta C_{p_f}$. Thus,
$\Delta C_{p_f}$ acts as an excitation force that is applied to the membrane.
The Fourier coefficients,
$\bar {\boldsymbol {F}}$, can now be computed from (2.15), and their substitution into the Laplace transform of (2.5) produces the membrane dynamic solution in the Laplace domain. This approach predicts the membrane dynamic response to any arbitrary motion or gust. However, a numerical Laplace inversion is required to obtain a solution in the time domain, as no analytical expression is available for the inverse Laplace transform of our problem. We apply the numerical scheme of Valsa & Branc̆ik (Reference Valsa and Branc̆ik1998) to carry out this inversion, which is robust and reliable for both oscillatory and non-oscillatory functions.
A solution can alternatively be determined in the frequency (Fourier) domain by setting
$s=\mathrm {i} k$, where
$k$ is the reduced frequency (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 292). This approach computes readily the steady-state response of the membrane wing to harmonic gusts or motions. However, this method cannot obtain the transient response of the membrane and will therefore only be used here for the harmonic cases and as a means of verification of the indicial lift responses obtained in the Laplace domain. Details of the application of the frequency-domain method are presented in the next section.
2.3.2. Frequency domain
The assumption of harmonic motion for all variables converts the membrane dynamic solution to the frequency domain, where, for example
$y(t,x)=\hat {y}(k,x)\, \mathrm {e}^{\mathrm {i} kt}$, and the hat denotes a complex-valued amplitude. Assignment of
$s=\mathrm {i} k$ into (2.15) yields this equation in the frequency domain,
where
$\hat {{{\boldsymbol{\mathsf{B}}}}}={{\boldsymbol{\mathsf{B}}}}(s=\mathrm {i} k)$,
$\hat {{{\boldsymbol{\mathsf{E}}}}}={{\boldsymbol{\mathsf{E}}}}(s=\mathrm {i} k)$,
$\hat {\boldsymbol {H}}=\boldsymbol {H}(s=\mathrm {i} k)$ and
$\hat {\boldsymbol {F}}$ is the vector of complex amplitudes of the Fourier coefficients
$\hat {F}_n (s)$,
$n=1\ldots N$. Once determined by (2.16), these Fourier coefficients produce the resulting membrane deformation through (2.5). Note that the constant matrix
${{\boldsymbol{\mathsf{A}}}}$ is unaffected by the shift from the Laplace (2.15) to the frequency (2.16) domain.
Frequency domain analysis is a natural approach to study the canonical unsteady aerodynamic scenarios of Theodorsen (harmonic oscillations) and Sears (sinusoidal gust). For these two scenarios, the frequency-domain solution yields the membrane dynamic response and the aerodynamic lift response; these results are compared against the Laplace domain results for verification. Integration of the Theodorsen and Sears harmonic functions over the entire frequency domain yields the indicial lift responses to a step change in angle of attack (Wagner's function) and to a sharp-edged gust (Küssner's function), respectively (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, pp. 284–287; Baddoo, Hajian & Jaworski Reference Baddoo, Hajian and Jaworski2021). Thus, the frequency-domain solution can generate all four canonical functions for the membrane aerofoil, whilst the transient dynamic response of the membrane can only be studied through the Laplace-domain solution. Case-specific technical details for this approach are further discussed in §§ 2.5.1 and 2.5.3 for harmonic heave motions and sinusoidal gusts, respectively.
2.4. Unsteady flow conditions
We next describe the unsteady flow conditions that the membrane aerofoil encounters for two cases of prescribed motion: generalised and harmonic heave motions, and a step change in angle of attack. We describe the generalised problem of a membrane aerofoil that encounters a small-amplitude transverse gust, and then focus on two canonical gust profiles of sinusoidal or sharp-edged geometry. For each of these cases, the aerodynamic load due to the unsteady flow,
$\Delta C_{p_f}$, is derived by applying
$w_a=w_{a_f}$ to (2.11). The normal velocity at the aerofoil surface,
$w_{a_f}$, is defined according to von Kármán & Sears (Reference von Kármán and Sears1938) and Wagner (Reference Wagner1925), where the amplitudes of the prescribed motions and gusts sustain the small-disturbance requirement
$|\alpha _0|\ll 1$. The expression obtained for
$\Delta C_{p_f}$ forms the
$\boldsymbol {H}$ vector in (2.15) for the Laplace-domain solution, or the
$\hat {\boldsymbol {H}}$ vector in (2.16) for the frequency-domain solution, which concludes the formulation of our problem.
2.4.1. Prescribed heave motion
Consider a membrane aerofoil that performs a prescribed translatory motion,
$h(t)$, that is normal to the flight direction, where
$h$ is normalised by
$b$ and is positive downwards. Under these conditions, the aerodynamic load on the membrane due to heave motion is obtained by substituting
$w_a(t,x)=\dot {h}(t)$ into (2.10), where the overdot denotes a time derivative. The solution of (2.10) in the Laplace domain yields a closed-form expression for the wake vorticity distribution, which is substituted into the Laplace transform of (2.11). Subsequent application of the coordinate transformation (2.3) yields the Laplace transform of the aerodynamic load due to harmonic heave oscillations:
where
$\bar {h}(s)$ is the Laplace transform of the prescribed heave motion, and
$\bar {\varPhi }(s)$ is the Laplace transform of Wagner's function (e.g. Sears Reference Sears1940),
Here
${C}(s)$ is the generalised Theodorsen function (e.g. Edwards Reference Edwards1979),
and
$K_0$ and
$K_1$ are modified Bessel functions of the second kind.
In the frequency domain, assuming
$h(t) = h_0\, \mathrm {e}^{\mathrm {i} kt}$, the amplitude of the effective angle of attack is
$\alpha _0=\mathrm {i} k h_0$, and the aerodynamic load due to harmonic heave oscillations is
where
$C(k)$ is the frequency-domain Theodorsen's function,
\begin{equation} C(k) = \frac{H_1^{(2)}(k)}{H_1^{(2)}(k)+\mathrm{i} H_0^{(2)}(k)}, \end{equation}
and
$H_0^{(2)}$ and
$H_1^{(2)}$ are Hankel functions of the second kind. The corresponding lift response is
which is, in fact, the rigid plate response to heave oscillations (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 272). The aerodynamic load expressions (2.17) and (2.20) form the term
$\Delta C_{p_f}$ in (2.14) for prescribed heave motions in the Laplace and frequency domains, respectively.
2.4.2. Step angle of attack
The canonical unsteady aerodynamics problem for the indicial lift response of an aerofoil to a step in angle of attack was originally solved by Wagner (Reference Wagner1925) for a rigid flat plate. The aerodynamic load on a membrane due to a step change in angle of attack may be computed using
$w_a(t,x)=\alpha (t)$, where
and
$\mathcal {H}(t)$ is the Heaviside function. The procedure detailed in § 2.4.1 is repeated to derive an expression for the pressure difference coefficient,
The corresponding lift coefficient due to a step change in angle of attack is
which recovers the indicial lift response found by Wagner (Reference Wagner1925).
2.4.3. Generalised transverse gust
We consider a membrane aerofoil that encounters a vertically oriented gust with an arbitrary profile. The leading edge of the aerofoil encounters the gust front at time
$t=0$ (figure 1), and the gust amplitude is assumed to be small with respect to the freestream velocity. The aerodynamic load due to the imposed unsteady flow,
$\Delta C_{p_f}$, is briefly detailed here to complete the formulation of the problem, as it is, in fact, the rigid aerofoil response to a transverse gust presented in Appendix B of Iosilevskii (Reference Iosilevskii2007). Note that
$\Delta C_{p_f}$ in (2.14) is denoted
$\Delta C_{p_g}$ in the present context of gusts.
The gust model assumptions permit the substitution of
$w_a(t,x)=\alpha _g(t-x-1)$ into (2.10) to produce the gust effect on the vorticity distribution along the wake, which leads to a closed-form expression for the aerodynamic load along the aerofoil using (2.11):
where
is the Laplace transform of Küssner's function,
$\varPsi (t)$ (e.g. Sears Reference Sears1940). Note that the aerodynamic load due to an arbitrary transverse gust is obtained by a convolution of Küssner's function and the time derivative of the gust angle of attack,
Equation (2.26) forms an expression for
$\Delta C_{p_f}$ in the Laplace domain, for an arbitrary, small-amplitude gust profile.
2.4.4. Sinusoidal gust
A sinusoidal gust encountered at the leading edge at time
$t=0$ can be expressed as
The substitution of the Laplace transform of (2.29) into (2.26) yields the aerodynamic load along the aerofoil due to a sinusoidal gust,
It is natural to solve the steady-state problem for harmonic gusts in the frequency domain by assuming
which yields
Here
\begin{equation} \hat{\varPsi}(k)=\frac{2}{\rm \pi}\frac{\mathrm{e}^{-\mathrm{i} k}}{k^2}\left(\frac{1}{H_1^{(2)}(k)+\mathrm{i} H_0^{(2)}(k)}\right) \end{equation}
is the Fourier transform of the Küssner function, which is obtained by assigning
$s=\mathrm {i} k$ in (2.27). The lift response to sinusoidal gusts in the frequency domain is
where the modified Sears function,
$S(k)$, is given by Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, p. 287)
Here
$\tilde {S}(k)$ is the classical Sears function whose gust front is at the mid-chord location at time
$t=0$, and
$J_0$ and
$J_1$ are Bessel functions of the first kind.
2.4.5. Sharp-edged gust
A sharp-edged gust is similarly described by
The substitution of (2.36) into (2.26) yields the aerodynamic load along the aerofoil due to a sharp-edged gust,
The resulting lift response is
which is the expected classical indicial lift of a rigid plate due to a sharp-edged gust. The distributed aerodynamic load (2.37) is, in fact, the external force applied on the membrane through the term
$\boldsymbol {H}$ in (2.15) in the case of a sharp-edged gust. This applied force initiates a membrane deformation which brings about a change in the aerodynamic load through aeroelastic coupling.
2.5. Unsteady lift response functions
The membrane unsteady lift coefficient is derived by integration of the aerodynamic load (2.14) along the membrane chord-line, which leads to
where the normalised lift due to membrane deformation is
Here,
$\varPhi (t)$ is the time-domain Wagner function, and
$f(t)$ and
$g(t)$ are functions of the Fourier coefficients given by
$$\begin{gather} f(t) = \frac{1}{2}\mathcal{F}_{1}(t) - \frac{1}{2}\mathcal{F}_{0}(t) -\frac{1}{4} \dot{\mathcal{F}}_0(t) -\frac{1}{4} \dot{\mathcal{F}}_1(t) + \frac{1}{4} \dot{\mathcal{F}}_2(t) + \sum_{m=2}^{N/2} \frac{\dot{\mathcal{F}}_{2m-1}(t)}{(2m-1)^2-1}, \end{gather}$$
$$\begin{gather}g(t) ={-}\frac{1}{4}\dot{\mathcal{F}}_0(t) + \frac{1}{4}\dot{\mathcal{F}}_2(t) - \frac{3}{16} \ddot{\mathcal{F}}_1(t) + \frac{1}{8} \ddot{\mathcal{F}}_3(t) +\frac{1}{2} \sum_{m=3}^{N/2} \frac{\ddot{\mathcal{F}}_{2m-1}(t)}{(2m-1)^2-1} , \end{gather}$$
where
$\mathcal {F}_n=F_n/\alpha _0$. The term
${\alpha }_0$ is the unsteady angle of attack amplitude in the harmonic cases, or the steady angle of attack in the indicial cases; the reader may consult § 2.4 for details on the definition of
${\alpha }_0$ and the lift due to the unsteady flow,
$C_{l_f}$, for each case of prescribed chord motion or gust encounter considered here.
The Laplace transform of
$C_{l_d}(t)$,
obtains the lift coefficient due to membrane deformation in the Laplace domain for indicial scenarios. Subsequent numerical Laplace inversion yields the membrane indicial lift response functions in the time domain. For cases of harmonic oscillations, in which the membrane solution is obtained in the frequency domain, the lift coefficient due to membrane deformation is expressed in the frequency domain
where
$\hat {f}(k)$ and
$\hat {g}(k)$ are found by substituting
$\mathcal {F}_n(t)=\hat {\mathcal {F}}_n(k)\, \mathrm {e}^{\mathrm {i} kt}$ into (2.41) and (2.42), respectively, and the auxiliary functions in the time domain become
$f(t)=\hat {f}(k)\, \mathrm {e}^{\mathrm {i} kt}$ and
$g(t)=\hat {g}(k)\, \mathrm {e}^{\mathrm {i} kt}$.
Substitution of
$C_{l_d}$ ((2.43) or (2.44) for the indicial or harmonic scenarios, respectively) and the case-specific
$C_{l_f}$ into (2.39) yields a closed-form expression for the total membrane lift coefficient,
$C_{l_m}$, from which extensions to the classical unsteady aerodynamic functions are derived after a solution for the Fourier coefficients is obtained.
2.5.1. Equivalent Theodorsen function
An equivalent Theodorsen function is constructed for a flexible membrane wing following the classical approach presented in Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, p. 279). An extension for Theodorsen's function is obtained by computing the membrane response to prescribed heave oscillations in the frequency domain:
where
$L_h^C$ is the circulatory lift due to heave oscillations, and
$L_{s_\alpha }$ is the static (aeroelastic) membrane lift-curve slope. The membrane circulatory lift is obtained by superposition between the circulatory lift due to membrane deformation (first term in (2.44)) and the circulatory lift due to the unsteady flow (first term in (2.22)). Substitution of the membrane circulatory lift expression into (2.45) yields a closed-form expression for the membrane Theodorsen function in the frequency domain,
The static membrane lift slope,
is established by direct integration of the static pressure load given by Nielsen (Reference Nielsen1963). Here,
$\mathcal {F}_{s_n}$ are the static membrane Fourier coefficients, normalised by the static angle of attack, which depend solely on the tension coefficient. Note that as
$k \rightarrow 0$ the Fourier coefficients in the unsteady membrane solution converge to the static solution,
${\hat {\mathcal {F}}_n \xrightarrow []{} \mathcal {F}_{s_n}}$, where
$\hat {\mathcal {F}}_n = F_n/\alpha _0$ and
$\alpha _0=\mathrm {i} k h_0$. Application of this limit to (2.46) after substituting the leading term in the asymptotic expansion of
$\hat {f}(k)$ in (B6) recovers
In other words, in the limit of small reduced frequencies the equivalent Theodorsen function converges to the standard Theodorsen function, as expected.
Finally, we note that the Fourier coefficients of the membrane solution converge to zero for a very stiff membrane as
$C_T\rightarrow \infty$ at any reduced frequency other than the fluid-loaded resonance frequencies. Under these conditions, the auxiliary function
$f(t)$ (2.41) goes to zero, and the static membrane lift-curve slope converges on
$2{\rm \pi}$. Therefore, the equivalent Theodorsen function recovers the rigid plate function for
$C_T\rightarrow \infty$, as expected. Appendix B reports further details on this limit.
2.5.2. Equivalent Wagner function
We next derive the equivalent Wagner function for a flexible membrane wing following Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, pp. 284–287). Note that the equivalent Theodorsen function derived in § 2.5.1 enables the computation of the equivalent Wagner function in the time domain through (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 285)
where
$L_{\alpha _0}^C$ is the circulatory lift of the membrane due to a step change in angle of attack, expressed in the time domain. The above equation allows for the computation of the equivalent Wagner function from both the Laplace-domain solution (first equality) or the frequency-domain solution (second equality). However, whilst the frequency-domain solution is more efficient when computing the Theodorsen function, the application of this solution to (2.49) requires a broad range of frequencies to obtain an accurate initial (high-frequency) response of the membrane lift. Therefore, the equivalent Wagner function is presented using the Laplace-domain solution.
Following the procedure presented in § 2.5.1 for the equivalent Theodorsen function, a closed-form expression is derived for the equivalent Wagner function in the Laplace domain,
where
$\bar {f}(s)$ is the Laplace transform of
$f(t)$ from (2.41). In the time domain, the equivalent Wagner function becomes
We note that
$\varPhi _m(0) = {{\rm \pi} }/{C_{l_{s\alpha }}}$ at time
$t=0$, and the static membrane lift slope,
$C_{l_{s\alpha }}$, is higher than the rigid plate lift slope due to the static membrane camber. Therefore, the initial value of the equivalent Wagner function is smaller than the classical value of
$1/2$ predicted by the standard Wagner function for a rigid flat plate.
At long times
$(t\rightarrow \infty )$, we can use the fact that
$\dot {f}(\tau )$ rapidly converges to zero as the membrane profile converges to the appropriate static profile. Thus,
\begin{align} \lim_{t\rightarrow\infty} \left\{\int_0^t \varPhi(t-\tau) \dot{f}(\tau) \,\mathrm{d}\tau\right\} &\cong \lim_{t\rightarrow\infty} \left\{\varPhi(t)\int_0^t \dot{f}(\tau) \,\mathrm{d}\tau\right\}\nonumber\\ &\cong \lim_{t\rightarrow\infty} \{\varPhi(t) f(t)\} \cong \left(\frac{C_{l_{s\alpha}}}{2{\rm \pi}}-1\right) \varPhi(t) \end{align}
and the equivalent Wagner function converges to the standard Wagner function,
${\varPhi _m(t\rightarrow \infty )\cong \varPhi (t\rightarrow \infty )=1}$, as expected. For a very stiff membrane
$(C_T\rightarrow \infty )$, the static membrane lift slope becomes
$2{\rm \pi}$, and the equivalent Wagner function converges to the rigid-plate Wagner function during the entire time response.
2.5.3. Equivalent Sears function
Extension of the canonical modified Sears function is next derived for a flexible membrane wing in the frequency domain, following the classical approach presented in Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, pp. 286–287). The equivalent Sears function is obtained by calculating the membrane response to oscillating gusts (§ 2.4.4) and normalising the expression by the membrane static lift:
where
$\hat {L}_{gs}(k)$ is the (complex) amplitude of the lift due to harmonic gust with reduced frequency
$k$.
Substitution of the membrane lift expression (2.39) into (2.53), superposing the lift due to membrane deformation (2.44) with the lift due to the sinusoidal gust (2.34) in the frequency domain, leads to a closed-form expression for the membrane equivalent Sears function,
The first term in (2.54) describes the rigid aerofoil lift response, the second term is the circulatory lift response due to membrane deformation, and the third term represents the non-circulatory lift response due to membrane deformation. Note that the unsteady membrane solution to an encounter with sinusoidal gusts converges to the static membrane solution as
$k\rightarrow 0$, namely
$\hat {\mathcal {F}}_n \xrightarrow [k \rightarrow 0]{} \mathcal {F}_{s_n}$, where
$\hat {\mathcal {F}}_n = \hat {F}_n/\alpha _0$. Thus, for very low reduced frequencies the equivalent Sears function converges to
\begin{align} S_m(k) &\cong \frac{2{\rm \pi}}{C_{l_{s\alpha}}} \left\{ S(k) + \left(\frac{C_{l_{s\alpha}}}{2{\rm \pi}} -1\right) \left[ 1 + k\left(\mathrm{i}\ln{\frac{k}{2}} -\frac{\rm \pi}{2}\right)\right] \right\}\nonumber\\ &\quad + \textit{O}( k\hat{\mathcal{F}}_3 ; k\hat{\mathcal{F}}_4 ; k^2\ln{k} ; k^2 ) \quad \text{as }\ k\rightarrow 0. \end{align}
We note that for the limiting case of
$C_T\rightarrow \infty$ the static membrane lift slope approaches
$2{\rm \pi}$ and the lift due to the membrane deformation converges to zero for
$k\neq \omega _{r_n}$. Therefore, the equivalent Sears function converges to the standard modified Sears function for
${C_T\rightarrow \infty }$, as expected. Appendix B reports further details on these low-frequency limits.
2.5.4. Equivalent Küssner function
The aeroelastic membrane Küssner function in the time domain,
follows from the equivalent Sears function determined in § 2.5.3 using the procedure outlined in § 2.5.2 for the Wagner function (e.g. Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 287). Here
$L_{{seg}}$ is the time-domain membrane lift response to a sharp-edged gust. The above equation enables the computation of the equivalent Küssner function from both the Laplace-domain solution (first equality) or the frequency-domain solution (second equality).
A closed-form expression for the equivalent Küssner function is determined in the Laplace domain by superposing the rigid aerofoil indicial lift (2.38) and the lift due to membrane deformation (2.43), with substitution of the resultant unsteady lift into the Laplace transform of (2.56):
The equivalent Küssner function in the time domain is clearly
We note that the initial value of the equivalent Küssner function is
${\varPsi _m(0) = 0}$. For
$t\rightarrow \infty$ the equivalent Küssner function asymptotically converges to unity, and for a very stiff membrane of
$C_T\rightarrow \infty$ the rigid-plate Küssner function is recovered.
3. Results and discussion
The membrane response to prescribed chord motion or an incoming gust is derived in both the time domain (via inverse Laplace transform) and the frequency domain, the latter of which is used to study the steady-state response to harmonic motions or sinusoidal gusts. Results are shown for four canonical cases: harmonic heave oscillations, step change in angle of attack, sinusoidal gust and sharp-edged gust, for which extensions of the four respective classical unsteady lift functions are presented for a flexible membrane wing. We begin by studying the lift and dynamic response of a nominal membrane of
$\mu =1$ and
$C_T=2.5$, followed by analysis of the role of each of the membrane parameters.
3.1. Prescribed motion
The membrane response to prescribed chord motion is derived for two canonical problems: harmonic heave oscillations, from which Theodorsen's function is derived, and a step change in angle of attack (also known as Wagner's problem). In both cases the membrane is free to deform around the chord-line, which adheres to the prescribed motion. Extensions of the classical Theodorsen and Wagner functions are presented for flexible membrane wings, along with a discussion on the membrane dynamic response to these unsteady flow conditions and the role of the membrane parameters
$(\mu, C_T)$ in its aerodynamic performance.
3.1.1. Harmonic heave oscillations
To assess the membrane wing response to prescribed oscillations in heave, we compute first the membrane amplitude at various reduced frequencies of oscillation,
$k$, for various tension coefficients and two mass ratios (figure 2). In addition, the resonance frequencies of the fluid-loaded membrane system (
$\omega _{r_1}, \omega _{r_2},$ etc.) are computed from the homogeneous system of (2.16), following the method of Kornecki, Dowell & O'Brien (Reference Kornecki, Dowell and O'Brien1976). The parametric dependence of the fluid-loaded resonance frequencies on
$C_T$ and
$k$ is illustrated with dashed red lines in figure 2. The left column of figure 2 presents maps of the resulting maximum membrane amplitude, obtained for
$\mu =1$ (upper row) and
$\mu =18$ (lower row) for varying tension coefficient and reduced frequency. As expected, significant amplitudes of oscillation occur for frequencies near the resonance frequencies of the fluid-loaded membrane. In the heavy membrane case,
$\mu =18$, where the mass ratio is encroaching upon the flutter instability threshold, predicted by Tiomkin & Raveh (Reference Tiomkin and Raveh2017) at
$\mu \ge 18.8$ for
$C_T=2$, the peaks in the maximal amplitude map are more concentrated, with a significantly increased amplitude along the second fluid-loaded resonance frequency. Note that these narrow peaks in figure 2(c) prevent the addition of the resonance frequencies to this plot, as these lines cover the peaks entirely; the relevant fluid-loaded resonance frequencies are plotted in figure 2(d) for reference. We further note that the presence of aerodynamic damping leads to finite amplitudes of the membrane at resonance in this linear analysis. The amplitude peaks along the fluid-loaded resonance frequencies reach large values that are beyond the validity range of the current study (especially in the heavy membrane case). However, away from these very narrow peaks the results across the rest of the frequency regime satisfy the ansatz of linear dynamics assumed by the present work.
Membrane amplitude in response to heave oscillations: (a,c) maximum amplitude maps obtained for various values of
$C_T$ and
$k$ for
$\mu =1$ and
$\mu =18$, respectively. (b,d) Present the membrane amplitude profiles obtained at the fluid-loaded resonance frequencies for
$\mu =1$ and
$\mu =18$, respectively. Background colour depicts maximum amplitude of membrane deformation, normalised by the heave motion amplitude,
$h_0$. Red dashed lines describe the first, second and third resonance frequencies of the fluid-loaded membrane, and black dotted lines present membrane in vacuo natural frequencies. A large value of the mass ratio,
$\mu =18$, is chosen to assess the membrane response near the onset of flutter, which Tiomkin & Raveh (Reference Tiomkin and Raveh2017) predict to occur for
$\mu \ge 18.8$ when
$C_T=2$. Aerodynamic damping leads to finite membrane amplitudes at resonance.

The right column of figure 2 shows the membrane amplitude profiles obtained along the fluid-loaded resonance frequencies, with background colour used to indicate the maximum amplitude value. Dotted black lines indicate the in vacuo natural frequencies,
$k_1, k_2,$ etc., where
$k_n=n{\rm \pi} \sqrt {C_T/8\mu }$, which are compared against the fluid-loaded resonance frequencies. Significant differences between the in vacuo and the fluid-loaded membrane resonance frequencies are obtained for the lowest mass ratio,
$\mu =1$, in figure 2(b). This difference in frequencies evokes oscillations with membrane amplitude profiles that are noticeably different from the membrane in vacuo modes. In addition, the membrane amplitude increases with reduced frequency, following the behaviour of the excitation term in (2.20) whose amplitude increases monotonically with
$k$. As the mass ratio increases in figure 2(d), the gap between the resonance and the structural frequencies diminishes, and second-mode oscillations become dominant; this is the first dynamically unstable membrane mode (Nielsen Reference Nielsen1963; Tiomkin & Raveh Reference Tiomkin and Raveh2017). In this case the membrane amplitude increases as the tension decreases (although the resonance frequency also decreases), which is typical of the membrane-wing behaviour on the verge of instability (Tiomkin & Raveh Reference Tiomkin and Raveh2017; Mavroyiakoumou & Alben Reference Mavroyiakoumou and Alben2020, Reference Mavroyiakoumou and Alben2021).
The differences between the in vacuo natural frequencies and the fluid-loaded resonance frequencies are mainly due to the added mass of the surrounding fluid, which must be taken into account when computing the total inertia of the coupled system. This effect may be quantified in non-dimensional terms as an added mass ratio,
$\mu _{{add}}$, which can be computed by assuming
$\omega_{r_1}={\rm \pi}\sqrt{C_T/(8\mu+8\mu_{add})}$. For a rigid plate the added mass is commonly taken as
$\mu _{{add}}={\rm \pi} /4$ (Katz & Plotkin Reference Katz and Plotkin2001, pp. 385–387). For a membrane wing, Alon Tzezana & Breuer (Reference Alon Tzezana and Breuer2019) found a constant added mass value of
$\mu _{{add}}=0.5$, and Minami (Reference Minami1998) determined that
$\mu _{{add}} = 0.68$ for a membrane oscillating in quiescent air. Minami (Reference Minami1998) used standing modes to describe the membrane deformation, without considering the tension along the membrane. However, Yadykin, Tenetov & Levin (Reference Yadykin, Tenetov and Levin2003) showed that the added mass of flexible plates is strongly affected by the mode of vibration. Figures 2(b) and 2(d) reveal in the current investigation that the membrane parameters and the reduced frequency of the harmonic heave motion control the amplitude profile of the oscillating membrane. Therefore, the added mass ratio in fact depends on both the membrane mass ratio and tension coefficient when considering the coupled problem of the membrane passive deformation in response to unsteady flow.
Figure 3(a) presents the membrane added mass ratio values for
$\mu =1$ and
$\mu =18$ as a function of the tension coefficient over
$1.73\le C_T\le 100$, which are compared against the known rigid plate added mass ratio of
${\rm \pi} /4$ and the results of Jaworski & Gordnier (Reference Jaworski and Gordnier2015) obtained for
$\mu =1.2065$. It is evident that the added mass ratio increases with the membrane mass ratio and decreases with increase in tension coefficient, where an asymptotic approach to the rigid plate solution as
$C_T\to \infty$ is noted. In addition, a good agreement is obtained with the results computed by Jaworski & Gordnier (Reference Jaworski and Gordnier2015, Table 1) assuming quasi-steady aerodynamics. Jaworski & Gordnier (Reference Jaworski and Gordnier2015) argued that circulatory effects are negligible in the computation of the resonance frequency of the system, and the close agreement in figure 3(a) substantiates this claim.
The dependence of non-dimensional added mass and aerodynamic damping on membrane tension coefficient for membrane mass ratios of
$\mu =1, 18$: (a) added mass ratio,
$\mu _{{add}}$, compared against the rigid plate added mass value of
${\rm \pi} /4$, and the added mass obtained by Jaworski & Gordnier (Reference Jaworski and Gordnier2015) for
$\mu =1.2065$; (b) damping coefficient,
$\zeta$, computed from the frequency ratio at the lift amplitude peak that sustains
$(k/\omega _{r_1})_{{peak}}=\sqrt {1-2\zeta ^2}$. The membrane added mass approaches the rigid plate value as
$C_T\to \infty$, and is in close agreement with Jaworski & Gordnier (Reference Jaworski and Gordnier2015). The aerodynamic damping is effectively constant for all of the examined values of
$C_T$, except at the lowest values near the divergence instability threshold,
$C_T\cong 1.73$, where the damping coefficient approaches
$\zeta =1/\sqrt {2}$. Resonance peaks do not occur for
$\zeta >1/\sqrt {2}$ (Rao Reference Rao2007, p. 274).

An aerodynamic damping coefficient,
$\zeta$, may also be computed for the fluid-loaded membrane using the frequency ratio at the peak of the unsteady lift amplitude value,
$(k/\omega _{r_1})_{{peak}}=\sqrt {1-2\zeta ^2}$ (Rao Reference Rao2007, pp. 271–274). Figure 3(b) plots the damping coefficient for two values of mass ratio as a function of the tension coefficient and shows that it is practically constant for
$C_T\ge 3$. The damping coefficient approaches the limit of
$\zeta =1/\sqrt {2}$ as the tension coefficient is further reduced, which is near the divergence instability threshold of
$C_T\cong 1.73$ (Tiomkin & Raveh Reference Tiomkin and Raveh2017). This limit describes the aerodynamic damping value beyond which no resonance peak is obtained, as would be expected for any harmonically forced linear system (Rao Reference Rao2007, p. 274). In general, all of the examined cases possess substantial aerodynamic damping, which explains the finite amplitudes obtained at the fluid-loaded resonance frequency conditions in the present linear analysis.
The membrane lift response to heave oscillations is next evaluated by comparing the membrane equivalent Theodorsen function (2.46) to the standard Theodorsen function of a rigid flat plate. Figure 4 illustrates this comparison for the nominal membrane of
$\mu =1$ and
$C_T=2.5$ as an Argand diagram (figure 4a) and in terms of modulus and phase (figure 4b). This representation of a complex-valued function as a two-dimensional plot is used to describe the unsteady lift amplitude and the phase lag relative to the heaving motion of the aerofoil. When
${\rm Im}{[C_m(k)]}<0$ in the Argand diagram, the lift response lags the rigid motion (negative phase), whereas the lift precedes the heaving motion (positive phase) when
${{\rm Im}{[C_m(k)]}>0}$. Whilst figures 4(a) and 4(b) are, in fact, two ways to describe the same results, each of these figures highlights different characteristics of the membrane lift response, as will be elaborated next.
Membrane lift response to heave oscillations, obtained for
$C_T=2.5$,
$\mu =1$ in terms of the membrane equivalent Theodorsen function: (a) Argand diagram; (b) modulus and phase. Frequencies of inflection points are denoted with black circles (
$k_{{inv}_1}$) and crosses (
$k_{{inv}_2}$), and resonance frequencies are denoted with red pluses (
$\omega _{{r}_1}$) and magenta crosses (
$\omega _{{r}_2}$).

For low reduced frequencies, the membrane lift response follows the general behaviour of the rigid plate response, with reduced amplitude and increased phase lag. As the reduced frequency of heave oscillations increases, at some point (typically for
$k$ smaller than the first resonance frequency) the membrane equivalent Theodorsen function changes its direction abruptly in the complex plane, where the unsteady lift response amplitude increases, rather than converging to zero as it would for rigid aerofoils. We mark this inflection point by reduced frequency
$k_{{inv}_1}$. With a further increase in
$k$ beyond
$k_{{inv}_1}$, a circular path is obtained in the complex plain plot until the next inflection point is reached at
$k=k_{inv_2}$, and so on. In figure 4(a), we present results for reduced frequencies up to the second in vacuo natural frequency,
$k_2\cong 3.5$, for the sake of clarity. Each of these circular arcs, which are identifiable in the Argand plot (figure 4a), contains one of the system's resonance frequencies for which a local maximum is observed in the lift response amplitude (figure 4b). A region of special interest is revealed around the first resonance frequency, for
$0.65\le k\le 0.96$, where the membrane aerofoil demonstrates a substantial increase in lift magnitude over a rigid aerofoil. Interestingly, oscillations with lower or higher reduced frequency (in the examined range of
$k\le 3.5$) result in substantial deficit in unsteady lift amplitude. Viewed in a practical context, these results for the lift amplitude may be used to extract the maximum unsteady wing load in the design process. Therefore, cases where the flexible membrane presents higher maximum loads than a rigid aerofoil could be hazardous when using predictions of the standard Theodorsen function, for example.
To further examine the origin of the circular paths in the membrane Theodorsen function, we recall that the equivalent Theodorsen function (2.46) is, in fact, a product of the standard Theodorsen function and a function of the membrane Fourier coefficients. Figure 5 presents the evolution of the first two Fourier coefficients (normalised by
$\alpha _0=\mathrm {i} k h_0$) with varying reduced frequency, obtained for the nominal membrane case. These normalised Fourier coefficients are the most dominant coefficients in the Fourier series used to describe the membrane slope (2.4), and are plotted as an Argand diagram (figure 5a) and in terms of its modulus and phase (figure 5b). Note that the unsteady solution recovers the static aeroelastic membrane results for
$k\rightarrow 0$, as expected, where the static results are marked by pentagrams in the Argand diagram. As the reduced frequency is increased from the static limit, the amplitudes of both of the normalised Fourier coefficients decrease at first, yielding a smaller amplitude of the oscillating membrane shape, and then increase as the reduced frequency approaches the first fluid-loaded resonance frequency. The first inflection point in the equivalent Theodorsen function corresponds to the first local minimum of
$|\hat {\mathcal {F}}_1|$, which is followed by a circular path in the complex plane plots of all coefficients (figure 5a). This entire frequency regime, in which the first circle appears in the Fourier coefficients, is dominated by the membrane's first mode, as is evident by the dominance of the first Fourier coefficient in figure 5(b). This dominance is maximum near the first resonance frequency and diminishes as
$k$ increases, which continues until the second resonance frequency is approached and the second mode coefficient becomes dominant. The frequency at which
$|\hat {\mathcal {F}}_2|$ crosses
$|\hat {\mathcal {F}}_1|$ is in fact
$k_{{inv}_2}$, which marks the transition from the first circle to the second circle in the equivalent Theodorsen function in figure 4(a). This behaviour, in which the inflection points are identified by a switch of dominance between the membrane modes, also continues to higher modes as the frequency is further increased, but is not shown here for the sake of clarity and brevity. Thus, it can be concluded that the circular arcs in the equivalent Theodorsen function are due to the membrane dynamic response, where each circle is related to a different dominant mode, and the inflection points between circles occur at the intersection between the modulus functions of two consecutive normalised Fourier coefficients.
The first two (most dominant) complex-valued Fourier coefficients, normalised by
$\alpha _0=\mathrm {i} k h_0$, for a nominal membrane of
$C_T=2.5$,
$\mu =1$ undergoing harmonic heave oscillations of reduced frequency
$k$:(a) Argand diagram; (b) modulus and phase. Static solutions are denoted by pentagram markers and recovered by the unsteady results as
$k\rightarrow 0$. First inflection point is marked with black circles and identified by the first local minimum of
$|\hat {\mathcal {F}}_1|$. Second inflection point is denoted by black crosses, identified by an intersection between
$|\hat {\mathcal {F}}_1|$ and
$|\hat {\mathcal {F}}_2|$. The insets in (a) introduce an enlarged view of the first inflection point region for illuminating the trend of both functions as
$k$ passes
$k_{{inv}_1}$.

Figure 6 presents the membrane amplitude profiles computed for a nominal membrane undergoing heave oscillations of varying reduced frequency, where the profile represents the amplitude of oscillation at every point along the membrane chord. A contour plot of the amplitude profiles is presented in figure 6(a) for varying reduced frequency,
$k$, with black dashed and dotted lines denoting the resonance and inflection point frequencies, respectively. Small deformations relative to the heave amplitude are obtained along the entire frequency range, except in the vicinity of the system's (fluid-loaded) resonance frequencies, for which large-amplitude profiles are observed with a shape similar to the membrane structural modes. This observation suggests a close coupling between the unsteady lift amplitude and the membrane amplitude in response to harmonic heave oscillations. Furthermore, a favourable lift is clearly achievable only for odd resonance frequencies (
$\omega _{r_1}, \omega _{r_3},$ etc.), which correspond to shapes that are symmetric around the mid-chord point, rather than anti-symmetric in the even-mode cases, as illustrated by comparing the lift response in figure 4(b) with the membrane amplitude profiles in figure 6(a).
Membrane dynamic response to harmonic heave oscillations, obtained for a nominal membrane of
$C_T=2.5$,
$\mu =1$. (a) Amplitude profiles computed for varying values of reduced frequency. Resonance frequencies are denoted with dashed black lines and inflection point frequencies are marked with dotted black lines. (b,c) Membrane amplitude profiles obtained for reduced frequencies around the first and second inflection points, respectively.

Figures 6(b) and 6(c) provide a more detailed view of the membrane amplitude profiles obtained for reduced frequencies near
$k_{{inv}_1}$ and
$k_{{inv}_2}$, respectively. In accordance with the normalised Fourier coefficients (figure 5), we note that small-amplitude profiles are obtained for small values of reduced frequency (figure 6b). These membrane shapes are convex, with a maximum camber point at the fore section of the aerofoil, in accordance with the static membrane solution (Nielsen Reference Nielsen1963). As the reduced frequency increases, the amplitude of the profile decreases until for
$k_{{inv}_1}$ an inflection point appears in the membrane profile and drastic changes in the profile shape are evoked with further increase in reduced frequency. These deformations indicate the excitation of the membrane structural modes as the reduced frequency approaches the system's first resonance frequency. For larger reduced frequencies near the second inflection point,
$k_{{inv}_2}$, figure 6(c) shows that significant membrane oscillations are evoked, in which the inflection point (
$k_{{inv}_2}$) represents the shift in the membrane amplitude profile, from the fully convex shape obtained for the first resonance frequency to the second mode shape obtained for the second resonance frequency. This shift is identified by an inflection point that appears near the leading edge of the membrane profile for
$k\ge k_{{inv}_2}$, after which the second mode of the membrane becomes most dominant. This result is in accordance with the behaviour of the Fourier coefficients, presented in figure 5, affirming the conclusion that the inflection points in the complex plane plot of the equivalent Theodorsen function mark the shift in dominance between consecutive membrane modes.
The circular paths of the membrane Theodorsen function provide direct insight into the dynamic interplay of the membrane and surrounding flow dynamics through the lens of the unsteady lift. Specifically, the aeroelastic oscillations of the membrane are connected to changes in the surrounding fluid flow and the shed wake, which are reflected in the unsteady lift response. For reduced frequencies near the fluid-loaded resonances, larger membrane deformations alter the shed flow vorticity to bring about longer response times, which lead to an increased phase lag between the unsteady lift and the prescribed motion. Further increases of the reduced frequency beyond the fluid-loaded resonance frequency yield decreases in the lift amplitude and phase lag due to a smaller camber of the deformed membrane amplitude. In addition, we note that the even membrane modes produce smaller circles about a center closer to the origin of the complex plot of the membrane Theodorsen function, which indicate a smaller lift amplitude than the odd membrane modes, as expected aerodynamically in light of the shape of these modes (figure 2b).
Figures 7 and 8 illustrate the separate effects of the tension coefficient and the membrane mass ratio on the membrane Theodorsen function. Tension coefficients between
$2$ and
$4$, and mass ratios between
$0.5$ and
$2.5$ are chosen to represent realistic membrane wings (e.g. Rojratsirikul, Wang & Gursul Reference Rojratsirikul, Wang and Gursul2010; Tiomkin & Raveh Reference Tiomkin and Raveh2021), while still remaining in the membrane stable regime, as the membrane loses stability via divergence for
$C_T<1.73$ and loses stability via flutter only in the case of heavy membranes of
$\mu >18.8$ (see Tiomkin & Raveh Reference Tiomkin and Raveh2017 for details). Results are presented for reduced frequencies up to the second in vacuo frequency
$(k=k_2)$ for the sake of clarity. The membrane stiffness is strongly influenced by the tension coefficient, and figure 7(b) shows that the membrane lift response to low-frequency oscillations approaches the rigid plate response with increasing tension coefficient, as expected. This result is further validated by examining an extreme case of
$C_T=50$, presented in figure 7(b) with a dashed red line. For this high tension coefficient the membrane is practically rigid, and indeed the resulting equivalent Theodorsen function follows closely the standard Theodorsen function for a wide range of frequencies up to about
$k=2$, where differences in the amplitudes arise due to parametric proximity to the system's first resonance frequency. Because the system's resonance frequencies increase with
$C_T$, the range of frequencies for which the equivalent Theodorsen function follows the standard Theodorsen function increases with
$C_T$ as well, and the first resonance circle occurs at a larger reduced frequency (i.e.
$k_{{inv}_1}$ increases). In addition, the circle diameter increases with
$C_T$ for
$\mu =1$, indicating that an increase in tension coefficient leads locally to an increased amplitude of the unsteady lift response near the resonance frequency, as the membrane oscillation amplitude also increases (figure 2b). Thus, in a practical sense, the hazardous region where a substantial increase in unsteady lift amplitude is obtained is clearly controlled by the tension coefficient, suggesting the possibility of optimising the flapping wing performance by controlling the tension along the membrane.
Effect of membrane tension coefficient on the membrane equivalent Theodorsen function for
$\mu =1$: (a) Argand diagram; (b) modulus and phase. The rigid-plate Theodorsen function is presented with a black line for comparison with the membrane equivalent Theodorsen function. The solution for a very large tension coefficient of
$C_T=50$ is presented with a red dashed line in (b) and indicates asymptotic convergence of the unsteady solution to the rigid plate solution as
$C_T \rightarrow \infty$.

Effect of membrane mass ratio on the membrane equivalent Theodorsen function for
$C_T=2.5$: (a) Argand diagram; (b) modulus and phase. The rigid-plate Theodorsen function is presented with a black line for comparison with the membrane equivalent Theodorsen function.

Analysis of the mass ratio effect on the membrane response to prescribed heave oscillations (figure 8) shows that at small reduced frequencies the mass ratio has practically no effect on the lift response, as expected. For higher reduced frequencies the first resonance circle appears earlier (at lower
$k$) as the mass ratio is increased, in accordance with the decrease in resonance frequency, while the amplitude of the lift response is practically unaffected. The hazardous region is also controlled by the mass ratio, as the peak in the lift amplitude follows the movement of the resonance frequency. However, in practical applications, this variable is often harder to control than the tension coefficient.
Some insight into the scaling of the membrane lift response to harmonic heave oscillations is gleaned from figures 9(a) and 9(b), which present the effects of the tension coefficient and the membrane mass ratio on the membrane Theodorsen function, respectively, as a function of the normalised reduced frequency,
$k/\omega _{r_1}$. For varying mass ratio, all of the examined cases collapse to a single curve for reduced frequency ratios up to
$k/\omega _{r_1}\cong 2.6$, beyond which the second fluid-loaded resonance peak is approached. Variation in the tension coefficient shows that the membrane Theodorsen function modulus peak at the first fluid-loaded resonance frequency is linearly proportional to the tension coefficient for all membranes of
$C_T\ge 2.5$; note that the modulus is normalised by
$C_T$ in figure 9(a). This dependence of the amplitude peak on
$C_T$ suggests that the aerodynamic damping that controls the peak amplitude is effectively constant for
$C_T\ge 2.5$, as was shown in figure 3(b) for
$\mu =1$. For lower values of the tension coefficient, a larger aerodynamic damping is obtained (
$\zeta \to 1/\sqrt {2}$), and the amplitude peak reduces significantly.
Effect of (a) membrane tension coefficient and (b) mass ratio on the membrane-equivalent Theodorsen function, obtained for
$\mu =1$ and
$C_T=2.5$, respectively. The modulus and phase of the equivalent Theodorsen functions are plotted against the reduced frequency of the prescribed motion, normalised by the first fluid-loaded resonance frequency. Aerodynamic damping leads to a finite peak in the modulus of the equivalent Theodorsen function, which scales on
$C_T$, and to a peak frequency that is smaller than the fluid-loaded resonance frequency.

3.1.2. Step angle of attack
We next examine the response of the nominal membrane to a step in angle of attack, in terms of its dynamic and aerodynamic response (figure 10). The membrane dynamic response begins with an initially taut profile, followed by membrane oscillations as a result of the abrupt change in angle of attack. These oscillations decrease in amplitude with time until a steady-state profile is obtained that is identical to the respective static solution (figure 10a). The resulting lift coefficient history (figure 10b) presents a similar trend to the membrane deformation history and suggests a close coupling between the two. In addition, when comparing the nominal membrane lift history in response to a step in angle of attack,
$C_{l_m}$, with the rigid plate response,
$C_{l_f}$, we see that the membrane wing achieves a larger lift across almost the entire response, converging to a value more than double that of the rigid-plate lift, due to membrane camber. Interestingly, during the initial transient stage of the response (for
$t<1.4$) the membrane lift is lower than the rigid plate lift. Namely, the membrane deformation due to the abrupt change in angle of attack produces negative lift, as evident by the plot of
$C_{l_d}^C$ and
$C_{l_d}^{{NC}}$ in figure 10(b), which describe the circulatory and non-circulatory terms in
$C_{l_d}$ (2.40). As the membrane is initially still and taut, the membrane surface accelerates in response to the sudden change in flow conditions, leading to a negative apparent mass lift at
$t=0$ and zero circulatory lift. The membrane inflates with time and its acceleration reduces, which yields a decrease in the non-circulatory lift magnitude. The initial acceleration subsequently yields a negative circulatory lift, but it also causes an increase in the membrane velocity, which increases the circulatory lift. This trend continues until at time
$t=1.4$ the circulatory lift is able to compensate the lift deficit due to
$C_{l_d}^{{NC}}$ and the membrane lift exceeds the rigid plate lift. For longer times the membrane deformation yields a higher circulatory lift that contributes to the further increase in the total membrane lift, as it proceeds to converge to the static membrane lift coefficient,
$C_{l_s}$.
Membrane (a) dynamic response and (b) lift response to a step in angle of attack, obtained for a nominal membrane of
$C_T=2.5$ and
$\mu =1$. Black line in (a) and dashed black line in (b) denote the static solution (Nielsen Reference Nielsen1963). The membrane unsteady lift coefficient,
$C_{l_m}$, is computed by superposition between the rigid plate indicial lift,
$C_{l_f}$, and the lift due to membrane deformation,
$C_{l_d}$, which is composed of the circulatory and non-circulatory terms,
$C_{l_d}^C$ and
$C_{l_d}^{{NC}}$, respectively.

Figure 11 illustrates the separate effects of the membrane tension coefficient and mass ratio on the membrane lift response to a step in angle of attack in terms of the equivalent Wagner function, and compares it to the standard Wagner function for rigid aerofoils. The equivalent Wagner functions were computed in the Laplace domain (2.50) and transformed to the time domain via numerical Laplace inversion. In general, for all of the examined cases, the membrane lift response is slower than the rigid plate response. However, we recall that the static membrane lift-curve slope is determined by the tension coefficient (2.47), and can be significantly larger than the rigid-plate lift slope due to aeroelastic camber. Thus, for example, while for
$C_T=2$,
$\mu =1$, at time
$t=100$ the equivalent Wagner function reaches only
$93.8\,\%$ of its steady-state solution, compared with
$99\,\%$ for a rigid flat plate, the steady-state lift in this case is substantially higher for the membrane wing (
$C_{l_{s\alpha }}\cong 28$), yielding a lift that is more than
$4$ times larger than the rigid plate lift. As the tension along the membrane is increased, the membrane Wagner function approaches the classical solution for a rigid plate, as expected (figure 11a). In addition, the initial value of the equivalent Wagner function is significantly lower than its rigid value, as it recovers the result predicted in § 2.5.2,
$\varPhi _m(0)={\rm \pi} /C_{l_{s\alpha }}$, marked by pentagrams in figure 11.
Effect of the (a) tension coefficient and (b) mass ratio on the membrane equivalent Wagner function, obtained for
$\mu =1$ and
$C_T=2.5$, respectively, as compared with the standard Wagner function (black line). Results are obtained via Laplace-domain solution and are verified against the expected initial values marked with pentagrams.

The effect of the mass ratio on the equivalent Wagner function is much less pronounced than the tension coefficient effect (figure 11b) and is discernible only for short time periods (i.e. high frequencies) when inertial effects are important. For long time durations, the lift responses are equivalent for any practical use.
3.2. Gust response
The membrane response to encounters with transverse gusts is now studied for two canonical cases: (i) a sinusoidal gust and (ii) a sharp-edged gust. The sinusoidal and sharp-edged gusts produce unsteady lift responses described by the equivalent Sears and Küssner functions, respectively. These extensions of the classical Sears and Küssner functions are presented for flexible membrane wings, along with discussion on the membrane dynamic response to these unsteady flow conditions and the role of the membrane parameters
$(\mu, C_T)$ in its aerodynamic performance. We note that while the sharp-edged gust has no physical meaning by itself, it is a very useful tool when predicting the aerofoil's response to an arbitrary (small-amplitude) transverse gust by appeal to convolution theory (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 288).
3.2.1. Sinusoidal gust
The membrane lift and dynamic response to an encounter with a sinusoidal gust is controlled by the membrane tension coefficient, mass ratio and the gust reduced frequency. The response of a nominal membrane wing to sinusoidal gusts of varying reduced frequency is analysed first, followed by a separate analysis of the effect of each of the membrane parameters on the resulting unsteady lift and dynamic response of the membrane.
Figure 12 presents the lift response of the nominal membrane to sinusoidal gusts in terms of the equivalent Sears function (2.54). This equivalent Sears function is compared against the classical modified Sears function for a rigid flat plate in figure 12(a) using an Argand diagram and in figure 12(b) in terms of the squared magnitude and the phase. The choice of a squared amplitude plot rather than a modulus plot follows the convention originated by Drischler (Reference Drischler1956) for harmonic gusts. For low reduced frequencies of
$k< k_{{inv}_1}\cong 0.41$, the membrane equivalent Sears function closely follows the classical modified Sears function, with a slightly decreased amplitude and increased phase lag. However, at the point of inflection (
$k=k_{{inv}_1}$) the lift amplitude begins to increase significantly with reduced frequency, creating a circular path in the complex plane, in a manner similar to the equivalent Theodorsen function (cf. figure 4). The lift amplitude increases beyond the rigid plate response in the vicinity of the first fluid-loaded resonance frequency (
$k=\omega _{r_1}$), which is followed by a sharp decrease in amplitude at higher reduced frequencies. This behaviour suggests the existence of a reduced frequency region for which gust mitigation is achievable using flexible membrane wings, while also revealing a range of frequencies (near the first resonance frequency) for which membrane flexibility could have adverse results.
Membrane lift response to sinusoidal gusts of various frequencies, obtained for a nominal membrane of
$C_T=2.5$,
$\mu =1$ in terms of the membrane equivalent Sears function: (a) Argand diagram and (b) squared modulus and phase. Frequencies of the inflection points are denoted with black circles (
$k_{{inv}_1}$) and crosses (
$k_{{inv}_2}$), and resonance frequencies are denoted with red pluses (
$\omega _{{r}_1}$) and magenta cross signs (
$\omega _{{r}_2}$).

To further study the origin of the inflection points in the complex plane plot of the equivalent Sears function, we recall that the equivalent Sears function (2.54) depends on the standard modified Sears function, the standard Theodorsen function and the Fourier coefficients used to describe the membrane deformation. Figure 13 presents the behaviour of the first two normalised Fourier coefficients obtained for the nominal membrane in response to sinusoidal gusts of varying reduced frequency. These Fourier coefficients are the most dominant coefficients in the membrane dynamic response for the range of reduced frequencies examined. For
$k\rightarrow 0$, the Fourier coefficients converge to the appropriate static solution, marked by pentagrams in figure 13(a). As the reduced frequency increases, the changes in both Fourier coefficients resemble the behaviours of the normalised coefficients in the heaving membrane case (cf. figure 5), while the amplification at the resonance frequency is less pronounced for the sinusoidal gust response. The first and second inflection points in the equivalent Sears function are identified here, similarly to the heaving membrane case, by locating the first local minimum in
$|\hat {\mathcal {F}}_1|$, and the crossing between
$|\hat {\mathcal {F}}_1|$ and
$|\hat {\mathcal {F}}_2|$ as the dominance transfers from the first membrane mode to the second mode around
$k_{{inv}_2}$. Namely, the circles in the complex plane plot of the equivalent Sears function are due to the membrane dynamic response, just like in the equivalent Theodorsen function, where each circle corresponds to a different dominant mode in the membrane oscillations. However, both the location of these circles in the frequency domain and the magnitude of the lift amplification differ significantly from the heaving membrane case, as the gust encounter introduces different flow mechanisms due to shedding of the unsteady flow conditions along the aerofoil.
The first two (most dominant) complex-valued normalised Fourier coefficients obtained for a nominal membrane of
$C_T=2.5, \mu =1$ that encounters a sinusoidal gust of reduced frequency
$k$: (a) Argand diagram; (b) modulus and phase. Static solutions are denoted by pentagram markers and recovered by the unsteady results as
$k\rightarrow 0$. First inflection point is marked with black circles and identified by the first local minimum of
$|\hat {\mathcal {F}}_1|$. Second inflection point is denoted by black crosses, identified by an intersection between
$|\hat {\mathcal {F}}_1|$ and
$|\hat {\mathcal {F}}_2|$.

Figure 14(a) presents the membrane amplitude profiles computed during steady-state oscillations of a nominal membrane that encounters sinusoidal gusts of various reduced frequencies. For
$k\to 0$ a convex amplitude profile is obtained, in accordance with the static membrane solution. Then, as the reduced frequency is increased, the membrane amplitude profile is flattened until
$k$ approaches the first resonance frequency, for which a large maximum amplitude is obtained. Figure 14(b) illustrates the membrane amplitude profiles computed for sinusoidal gusts at reduced frequencies near the first inflection point frequency,
$k_{{inv}_1}$. For reduced frequencies smaller than
$k_{{inv}_1}$ a convex shape with a maximum camber point at the fore part of the aerofoil is obtained. As the reduced frequency increases to
$k_{{inv}_1}$ the maximum amplitude of the membrane decreases, and the maximum camber point slowly approaches the mid-chord location. A further increase in the reduced frequency beyond
$k_{{inv}_1}$ yields a sudden shift of the maximum camber point downstream, as the membrane profile bears a close resemblance to the first unstable eigenshape of the membrane in cases of divergence instability, as reported by Sygulski (Reference Sygulski2007) and Tiomkin & Raveh (Reference Tiomkin and Raveh2017). This shift in the membrane amplitude profile signals the excitation of the membrane structural modes as the first fluid-loaded resonance frequency is approached. As we further increase the reduced frequency to the vicinity of the second inflection point frequency,
$k_{{inv}_2}$, in figure 14(c) we see a clear change in the membrane amplitude profile from a shape that is dominated by the first structural mode to a shape in which the second structural mode is most dominant, as supported by the Fourier coefficients in figure 13. Thus, any inflection point in the complex plane plot of the equivalent Sears function is related to a shift in dominance between two consecutive membrane mode shapes. As the gust frequency is increased, higher membrane modes become dominant. We further note that for the frequency regime in which membrane oscillations amplify the lift response (around
$k=\omega _{r_1}$), large amplitude deformations are obtained with a convex amplitude profile. For higher reduced frequencies, the membrane amplitude profile is no longer convex due to the appearance of additional nodal points along the profile, and smaller maximum camber is obtained. Membrane oscillations with these amplitude profiles attenuate the aerofoil's lift response, presenting a reduced lift amplitude relative to the rigid plate lift (figure 12b).
Membrane dynamic response to an encounter with a sinusoidal gust, obtained for a nominal membrane of
${C_T=2.5}$,
${\mu =1}$. (a) Amplitude profile computed for varying values of reduced frequency. Resonance frequencies,
$\omega _r$, are denoted with dashed black lines and frequencies of inflection points,
$k_{inv}$, are denoted by dotted black lines. Membrane amplitude profiles near the (b) first and (c) second inflection points.

Figure 15 illustrates the effect of the tension coefficient on the equivalent Sears function, which is presented as a curve in the complex plane and in terms of its squared amplitude and phase, as compared with the classical modified Sears function. These results are obtained using the frequency-domain solution, which are verified by comparison with results of the Laplace domain solution (see dashed red line in figure 15a). The equivalent Sears function approaches the classical rigid aerofoil function uniformly at low reduced frequencies as the tension coefficient,
$C_T$, increases. An extreme case of
$C_T=50$, in which the membrane is practically rigid, is presented in figure 15(b) to confirm the approach of the analytical solution to that of the rigid aerofoil at large
$C_T$. The equivalent Sears function of this extreme case follows closely the rigid-plate modified Sears function for a wide range of reduced frequencies up to about
$k=2$ and validates our solution for large tension coefficients. As the tension coefficient is increased from two to four in figure 15, the first inflection point is delayed to a higher reduced frequency while also producing a larger unsteady lift amplification near the first resonance frequency. The beneficial region, where significant gust mitigation is obtained, is also shifted to higher values of
$k$ with an increase in the tension coefficient. This unique characteristic of the membrane wing presents a special opportunity to calibrate the pre-tension on the membrane to successfully alleviate gusts of specific frequency regime, which is known to be one of the crucial challenges in SUAV design (Mohamed et al. Reference Mohamed, Massey, Watkins and Clothier2014).
Effect of tension coefficient on the membrane equivalent Sears function, for
$\mu =1$: (a) Argand diagram; (b) squared amplitude and phase, as compared with the rigid-plate modified Sears function (black line). An additional solution for a very large tension coefficient of
$C_T=50$ is presented with a magenta dashed line in (b), to validate the unsteady solution by convergence to the rigid plate solution. All results are obtained with the frequency-domain solution, presented for frequencies up to
$k=k_2$ for clarity. Direct comparison with the Laplace-domain solution is made for
$C_T=2$ (red dashed line in (a)), where the strong agreement indicates the equivalence of the methods.

Looking into the effect of the mass ratio on the membrane lift response to a sinusoidal gust, figure 16 presents the equivalent Sears functions which are computed for varying values of the mass ratio,
$\mu$, and are compared with the classical modified Sears function. For very low frequencies, up to the point of inflection
$k_{{inv}_1}$ (whose value decreases with increase in
$\mu$), the lift response is practically unaffected by changes in the mass ratio, as predicted by the asymptotic analysis (2.55). Thus, the response can be considered quasi-steady for this frequency regime. As the gust frequency is increased beyond the inflection point, a circular path appears at the complex plane plot with a radius (and amplitude) that slightly increases with
$\mu$. The membrane oscillates with the gust frequency so that an increase in
$\mu$ leads to increase in the membrane inertia, which is proportional to
$\mu k^2$, resulting with a larger amplitude of oscillation for a given gust frequency. However, as the system's resonance frequencies decrease with increase in
$\mu$ (cf. figure 2), the lift amplification at the first resonance frequency increased only slightly with an increase in
$\mu$. We note that the in vacuo natural frequencies of the membrane are proportional to
$\mu ^{-\frac {1}{2}}$, indicating that for these frequencies the mass ratio has no effect on the membrane amplitude. Thus, the small effect of the mass ratio on the lift amplitude, illustrated in figure 15, is attributed to the added mass of the fluid-loaded membrane. For frequencies beyond the first resonance frequency, the gust mitigation region is controlled by the membrane mass ratio in a manner that is similar to the effect of the tension coefficient (cf. figure 15b), where a decrease in
$\mu$ delays the region to higher values of
$k$. However, in contrast to the tension coefficient, this variable is typically fixed in membrane wing applications, and thus is expected to be less useful for aerodynamic design purposes.
Effect of membrane mass ratio on the membrane equivalent Sears function plotted,
$C_T=2.5$: (a) Argand diagram; (b) squared amplitude and phase, as compared with the rigid-plate modified Sears function (black line). All results are obtained with the frequency-domain solution, presented for frequencies up to
$k=k_2$ for clarity. Results from the Laplace-domain solution are shown in (a) with a red dashed line, presenting good agreement with the frequency-domain solution.

3.2.2. Sharp-edged gust
The dynamic response of a nominal membrane wing to a sharp-edged gust is presented in figure 17 in terms of its deformation in time and its lift response. Here we choose to present the actual lift response first (instead of the normalised response represented by the Küssner function) to obtain quantitative conclusions on the difference in the lift produced by a flexible membrane wing (blue line) and a rigid flat plate (black line) during sharp-edged gust encounter. While the transient response up to
$t\cong 30$ introduces oscillations in the membrane deformation, resulting in oscillations in the lift response, at a later time a fully convex shape is obtained converging to the appropriate static solution (Nielsen Reference Nielsen1963). The lift response in figure 17(b) shows that, as expected, the membrane achieves a much higher lift coefficient than a rigid flat plate, converging to more than double the lift coefficient of the flat plate, as in the case of a step change in angle of attack (cf. figure 10). However, zooming in onto the transient lift response during gust penetration reveals that at initial stage, for
$t<1.7$, the membrane presents a lower lift coefficient than the rigid plate. The inset in figure 17(b) shows that the transient membrane deformation leads to a negative contribution to the lift coefficient for
$t<1.7$, which results in a reduced lift coefficient compared with the rigid plate indicial lift. As the gust front approaches the trailing edge, the circulatory lift due to membrane deformation,
$C_{l_d}^C$, increases first until it overcomes the non-circulatory lift,
$C_{l_d}^{{NC}}$, for
$t=1.7$. We observe that this time is longer than the time required to compensate for the initial lift deficit in the membrane response to a step change in angle of attack (cf. figure 10b). From this moment on, the membrane lift coefficient surpasses the rigid-plate lift and slowly converges to the static solution.
Membrane (a) dynamic response and (b) lift response to a sharp-edged gust, obtained for a nominal membrane of
$C_T=2.5$ and
$\mu =1$. Black line in (a) and dashed black line in (b) denote the static solution (Nielsen Reference Nielsen1963). The membrane unsteady lift coefficient,
$C_{l_m}$, is computed by superposition between the rigid plate indicial lift,
$C_{l_g}$, and the lift due to membrane deformation,
$C_{l_d}$, which is composed of the circulatory and non-circulatory terms,
$C_{l_d}^C$ and
$C_{l_d}^{{NC}}$, respectively. For
$t<1.7$ the membrane deformation yields negative lift that reduces the total membrane lift compared with the rigid plate indicial lift. After time
$t=1.7$ the membrane lift surpasses the rigid plate lift due to induced membrane camber and converges to the static solution.

Figure 18 presents the membrane profiles as obtained for
$t\le 2$, during gust penetration when the gust front travels along the chord. The membrane encounters the gust at the leading edge at time
$t=0$, in a taut initial position. As the gust front moves downstream, a small hump appears in the membrane profile near the leading edge, which increases in size and moves downstream with the advancement of the gust. At time
$t=1.7$, for which the membrane lift surpasses the rigid plate lift, the gust front has not yet reached the trailing edge, but a fully convex shape has developed with a maximum camber point at the aft part of the aerofoil. Note that the membrane profile obtains a positive camber during gust penetration, which would result with increased static lift compared to a rigid flat plate. However, the unsteady response of the membrane, and in particular its acceleration, produces a negative non-circulatory lift at the initial response of the membrane (figure 17b), highlighting the importance of a full unsteady aerodynamic model for predicting the unsteady lift response of the membrane.
Membrane deformation in response to an encounter with a sharp-edged gust, obtained during gust penetration (when the gust front travels along the chord during
$0\le t\le 2$) for a nominal membrane of
${C_T=2.5}$,
$\mu =1$. Red points denote maximum camber point at each time step, and a dashed magenta line is used to identify the membrane profile at time
$t=1.7$ (also denoted in figure 17b for the lift response).

The membrane equivalent Küssner function is presented in figure 19 to study separately the effect of the tension coefficient and mass ratio. These plots were computed using the Laplace-domain solution (2.57) and were verified by comparison to the frequency-domain solution (2.56). In general, the membrane response to a sharp-edged gust is slower than the rigid plate response, similarly to the Wagner function case. As the tension coefficient is reduced, a slower response is obtained since a larger camber profile is achieved at steady state, which takes a longer time to attain. Contrarily, the membrane mass ratio appears to affect only the initial oscillatory stage, while the rest of the response is practically independent of the mass ratio, as expected.
Effect of (a) tension coefficient and (b) mass ratio on the membrane equivalent Küssner function, obtained for
$\mu =1$ and
$C_T=2.5$, respectively, and compared with the standard Küssner function for a rigid aerofoil (black line). Results are obtained via Laplace domain solution and are validated against the frequency domain solution (red dashed line).

4. Concluding remarks
Unsteady aerodynamic theory is extended to include aeroelastic coupling for predicting the lift response of flexible membrane wings to small-amplitude chord motions or transverse gusts. The pressure loads and membrane deformations due to dynamic fluid–membrane coupling are determined in the Laplace domain for cases of a step change in angle of attack or gust profile, and in the frequency domain for harmonic heave motions and sinusoidal transverse gusts. In addition, the numerical inversion scheme of Valsa & Branc̆ik (Reference Valsa and Branc̆ik1998) is used to evaluate the indicial lift responses in the time domain. Lift responses computed for the canonical unsteady aerodynamic scenarios of harmonic aerofoil motions or gusts, as well as of a step change in the angle of attack or gust profile, constitute aeroelastic extensions to the classical Theodorsen, Sears, Wagner and Küssner functions, respectively, for a membrane aerofoil.
In the harmonic scenarios, the membrane-equivalent Theodorsen and Sears functions follow the parametric trends of their classical, rigid aerofoil counterparts at low reduced frequency,
$k$, albeit with a reduced unsteady lift amplitude and an increased phase lag. As the reduced frequency increases and approaches the first resonance of the fluid-loaded membrane, the membrane-equivalent Theodorsen and Sears functions introduce distinct circular paths in the complex plane, which initiate at different values of
$k$ for each function. Closed-form expressions for both functions reveal that these circular paths are related to the membrane dynamic response during oscillations, where each circle corresponds to a different dominant mode, and the inflection points that connect the circles represent the shift in dominance between two consecutive membrane modes. The membrane deformation shape differs for each of the unsteady scenarios, which is also reflected in the unsteady lift functions. For example, in the heaving membrane case, the unsteady lift lags the aerofoil motion for the considered membrane parameters, whereas the unsteady lift leads the sinusoidal gust for reduced frequencies near the second fluid-loaded resonance, albeit with a significantly smaller lift amplitude.
The model results for harmonic aerofoil motions or incoming flow disturbances suggest parametric regions where the aeroelastic response of the membrane could enable performance improvements for flapping flight or gust resilience. The unsteady lift amplitude of the membrane in each of these scenarios is higher than that of the rigid flat plate for a range of reduced frequencies in the neighbourhood of the first fluid-loaded membrane resonance. Thus, in this frequency range the standard Theodorsen and Sears functions underestimate the load on the aerofoil. This parametric region is controlled by the tension coefficient, whereby the aerodynamic load on the aerofoil may be enhanced (or reduced) through either passive or modest active control of the membrane pretension. For reduced frequencies in a finite range above the first resonance frequency, the present model predicts that membrane oscillations attenuate the unsteady lift response to sinusoidal gusts or harmonic heave motions. In addition, by increasing the membrane pretension, this beneficial region shifts to higher reduced frequencies, which could be exploited as a gust mitigation strategy in practice.
Closed-form expressions derived in the Laplace domain for the membrane-equivalent Wagner and Küssner functions reveal the direct dependence of these functions on the original rigid-plate functions and on the Fourier coefficients that describe the membrane profile. Each indicial lift response of the membrane aerofoil is slower relative to the rigid plate response but results in a significantly higher lift in the steady state due to aeroelastic membrane deformation. The membrane initial response in short times to an abrupt change in angle of attack or to an encounter with a sharp-edged gust produces a negative non-circulatory lift due to the acceleration of the membrane profile from a still and taut position. Therefore, the overall initial lift response of the membrane is smaller than the rigid plate response due to the gradual elastic reaction of the membrane to the changes in the fluid flow. However, at later times the non-circulatory lift due to the abrupt change in the flow field weakens and the circulatory lift increases, such that the membrane lift quickly overtakes the rigid plate lift response and converges to the static solution. While the membrane Küssner function always begins from a zero value at the initial time when the gust front meets the membrane leading edge, the initial value of the membrane Wagner function varies with the membrane tension coefficient, as it depends on the membrane static lift slope. As the tension coefficient is increased, the initial value of the membrane-equivalent Wagner function also increases and approaches the classical rigid-aerofoil value of
$1/2$.
Results from the present theoretical effort invite computational and experimental companion efforts to elucidate the practical range of validity of the model, including the influence of nonlinear flow effects that occur at large gust ratios (e.g. Andreu-Angulo et al. Reference Andreu-Angulo, Babinsky, Biler, Sedky and Jones2020; Jones Reference Jones2020), which are expected to inform future improvements to the predictive aeroelastic framework.
Funding
This work was supported by the Zuckerman-CHE STEM Leadership Program, with partial support from the National Science Foundation under award 1846852.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Mathematical identities
The following identities are used in the mathematical derivations of § 2:
$$\begin{gather} 1=\frac{2}{\rm \pi}\sum_{n=1}^\infty \frac{1-({-}1)^n}{n}\sin n\theta , \quad 0<\theta<{\rm \pi} , \end{gather}$$
$$\begin{gather}\cos\theta=\frac{8}{\rm \pi}\sum_{m=1}^\infty \frac{m}{4~m^2-1}\sin 2m\theta , \quad 0<\theta<{\rm \pi} , \end{gather}$$
$$\begin{gather}\varLambda_1(x,\xi)=\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{gather}$$Appendix B. Asymptotic analysis of the lift due to membrane deformation for low-frequency regime
This appendix presents an analysis of the lift due to membrane deformation in the asymptotic regime of low reduced frequencies. The analysis begins with an identification of the leading terms in
$k$ for the standard Theodorsen function,
$C(k)$, and the two auxiliary functions,
$\hat {f}(k)$ and
$\hat {g}(k)$, as
$k\rightarrow 0$. These terms are then combined to obtain the leading terms in the lift coefficient,
$\hat {C}_{l_d}$, for low reduced frequency. A series expansion is applied to the standard Theodorsen function,
$C(k)$, leading to
where
$\gamma _e$ is the Euler constant. For the auxiliary functions,
$\hat {f}(k)$ and
$\hat {g}(k)$, we derive the leading terms up to
$\textit {O}(k\hat {\mathcal {F}}_3 , k\hat {\mathcal {F}}_4)$ and
$\textit {O}(k^2)$, respectively,
We have used the fact that the magnitude of the Fourier coefficients
$\hat {\mathcal {F}}_n(k)$ is strongly reduced for
$n>2$ in the low reduced frequency regime, as seen in figure 20(a), in which the Fourier coefficients obtained for a nominal membrane that encounters a sinusoidal gust of low reduced frequency are presented. For
$k\rightarrow 0$ the first Fourier coefficient is an order of magnitude larger than the second Fourier coefficient, which is at least one order of magnitude larger than the rest of the coefficients. Figure 20(b) illustrates the behaviour of the product
$k|\hat {\mathcal {F}}_n|$ for the first seven Fourier coefficients in the low reduced frequency regime, and compares it with
$k^2$ and
$k^2\ln {k}$. As
$k\rightarrow 0$,
$k|\hat {\mathcal {F}}_3|$ and
$k|\hat {\mathcal {F}}_4|$ are comparable and are an order of magnitude larger than
$k^2\ln {k}$. However, as
$k$ increases this difference diminishes and at
$k\cong 0.005$ the three functions are comparable.
Fourier coefficients obtained for a nominal membrane wing that encounters a sinusoidal gust of low reduced frequency: (a) modulus, (b) modulus multiplied by the reduced frequency,
$k$, and compared against
$k^2$ and
$k^2\ln {k}$ denoted with a black dashed and dotted lines, respectively. Static solution is denoted in (a) by pentagram markers and is recovered by the unsteady solution for
$k\rightarrow 0$. Only seven of the
$24$ coefficients are presented for the sake of clarity, as the magnitude of higher mode coefficients is negligible for low reduced frequencies.

Further simplification of the above expansion is obtained by recalling that
and
which when substituted into (B2) and (B3), respectively, produces
Substitution of (B1), (B6) and (B7) into (2.44) yields
which represents the behaviour of the lift due to membrane deformation in response to a low-frequency flapping motion or gust encounter. We note that the lift amplitude is controlled by the difference between the static membrane lift slope and the rigid plate lift slope,
$C_{l_{s\alpha }}-2{\rm \pi}$, as this term describes the static lift due to membrane camber, which is recovered by applying
$k=0$ to the unsteady problem. The tension coefficient on the membrane controls the static membrane lift-slope and as the tension increases the lift slope decreases due to decreased camber (increased stiffness) and approaches
$2{\rm \pi}$. In the rigid-membrane limit of
$C_T\rightarrow \infty$, the lift due to membrane deformation converges to zero as the reduced frequency approaches zero, as expected. In addition, it is obvious that the membrane inertia has no notable role in the low-frequency lift response of the membrane.












































































































