# Instability and the post-critical behaviour of two-dimensional inverted flags in axial flow

Mohammad Tavallaeinejad, Michael Païdoussis, Mathias Legrand, Mojtaba Kheiri

► To cite this version:

Mohammad Tavallaeinejad, Michael Païdoussis, Mathias Legrand, Mojtaba Kheiri. Instability and the post-critical behaviour of two-dimensional inverted flags in axial flow. Journal of Fluid Mechanics, 2020, 890, 10.1017/jfm.2020.111. hal-02520870

HAL Id: hal-02520870

https://hal.science/hal-02520870

Submitted on 27 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Instability and the post-critical behaviour of two-dimensional inverted flags in axial flow

Mohammad Tavallaeinejad<fnref n="1" />, Michael P.Païdoussis<fnref n="1" />†, Mathias Legrand<fnref n="1" />, and Mojtaba Kheiri<fnref n="2" />

Inverted flags – clamped-free elastic thin plates subjected to a fluid flowing axially and directed from the free end towards the clamped end – have been observed experimentally and computationally to exhibit large-amplitude flapping beyond a critical flow velocity. The motivation for further research on the dynamics of this system is partly due to its presence in some engineering and biological systems, and partly because of the very rich dynamics it displays. In the present paper, the goal is to develop a nonlinear analytical model for the dynamics and stability of high aspect ratio (i.e. height to length ratio) flags. The inviscid fluid flow is modelled via the quasi-steady version of Theodorsen’s unsteady aerodynamic theory, and the Polhamus leading-edge suction analogy is utilized to model flow separation effects from the free end (leading edge) at moderate angles of attack. Gear’s backward differentiation formula and a pseudo-arclength continuation technique are employed to solve the governing equations. Numerical results suggest that *fluidelastic instability* may be the underlying mechanism for the flapping motion of high aspect ratio *heavy* inverted flags. In other words, flapping may be viewed as a *self-excited vibration*. It was found from numerical results that the undeflected static equilibrium of the inverted flag is stable at low flow velocities, prior to the occurrence of a supercritical pitchfork bifurcation. The pitchfork bifurcation is associated with static divergence (buckling) of the flag. At higher flow velocities, past the pitchfork bifurcation, a supercritical Hopf bifurcation materialises, generating a flapping motion around the deflected static equilibrium. At even higher flow velocities, flapping motion becomes symmetric, around the undeflected static equilibrium. Interestingly, it was also found that heavy flags may exhibit large-amplitude flapping right after the initial static equilibrium, provided that they are subjected to a sufficiently large disturbance. Moreover, inverted flags with a non-zero initial angle of attack were found to be less stable than their perfectly flow-aligned counterparts.

**Key words:** fluid-structure interactions, inverted flags

## 1. Introduction

Exploring the dynamics of a thin, flexible plate-like structure subjected to axial flow has been of constant interest to researchers as it can enhance our understanding of natural processes, such as the reconfiguration of plants in wind (Gosselin et al. 2010), snoring in humans (Huang 1995a; 1995b) and flight of insects (Sane 2003). It can also help in the design of trouble-free engineering systems, such as aircraft and paper handling machines. The body of literature on the subject is enormous and cannot be fully covered in this paper. For a comprehensive reviews of the subject, the reader is referred to Dowell (1975), Mei et al. (1999), Shelley & Zhang (2011) and Païdoussis

<footnote n="1">¹Department of Mechanical Engineering, McGill University, Montréal, Québec, Canada H3A 0C3</footnote>

<footnote n="2">²Department of Mechanical, Industrial and Aerospace Engineering, Concordia University, Montréal, Québec, Canada H3G 1M8</footnote>

<footnote>† Email address for correspondence: michael.paidoussis@mcgill.ca</footnote>
(2016). Here, we focus only on a subset of fluid-structure interaction (FSI) problems involving a flexible thin plate in axial flow: a cantilevered thin plate (or a flag) of length *L* and height *H* subjected to a fluid flowing axially with velocity *U* and directed from the free end towards the clamped one, otherwise known as an ‘inverted flag’; see figure 1.

The first investigation on the dynamics and stability of inverted flags may be attributed to Guo & Païdoussis (2000). They studied theoretically the linear stability of rectangular plates with various boundary conditions, such as clamped-free, free-free and free-clamped, in inviscid, incompressible two-dimensional (2D) flow. They found that a free-clamped (or inverted) plate flutters via a Hopf bifurcation at a non-zero critical flow velocity which is inversely proportional to the fluid-to-plate mass ratio defined as $\mu = \rho_f L/\rho_p h$; $\rho_f$ and $\rho_p$ being the mass density of the fluid and plate, respectively, and $h$ is the thickness of the plate.† This differs from the divergence instability proposed by Kim et al. (2013), and reported by Sader et al. (2016b), Gurugubelli & Jaiman (2015) and Goza et al. (2018). Recently, interest in the dynamics of inverted flags has been revived mainly due to their potential application to small-scale energy harvesting; see, for example, the experimental study by Orrego et al. (2017) on inverted piezoelectric flags. Moreover, exploring the dynamics of inverted flags is expected to improve our understanding of the flapping of biological structures, such as plant leaves in wind (Sader et al. 2016a; Fan et al. 2019).

Kim et al. (2013) explored the dynamics of inverted flags in an open-loop wind tunnel at high Reynolds numbers $Re \equiv UL/\nu \approx 10^4 - 10^5$, $\nu$ being the kinematic viscosity of the fluid. They identified three consecutive regimes/modes of dynamical behaviour, namely the stretched-straight mode, the flapping mode and the fully deflected mode, as the flow velocity in the test-section was increased. It was found that, in contrast to a typical flag which undergoes flapping beyond a critical flow velocity, the inverted flag performs flapping only within a finite range of flow velocities, somewhat analogously to vortex-induced vibration of a spring-supported cylinder in cross-flow. It was also found that the inverted flag displays larger peak-to-peak flapping amplitude *A* (up to $A/L = 1.8$) than the conventional flag.

Following the seminal work by Kim et al. (2013), several researchers employed fully computational approaches to solve the fluid-structure interaction problem of inverted flags. For example, Ryu et al. (2015) used the immersed boundary method (IBM) to investigate the flapping dynamics of 2D inverted flags at relatively low Reynolds numbers (up to $Re = 250$). For $Re < 50$, the stretched-straight and flapping modes did not arise. For $100 \le Re \le 250$, however, the stretched-straight, flapping and deflected modes occurred in a regular sequence as the bending rigidity parameter of the flag, which is the reciprocal of the dimensionless flow velocity, was decreased. They found numerically that the deformation of the free end of the flag grows steadily, as the bending rigidity is decreased in the range between the stretched-straight and the flapping modes. They called this intermediate regime “the biased mode”.

vAlmost concurrently with Ryu et al. (2015), two other research groups performed computational investigations on the dynamics of inverted flags. Gilmanov et al. (2015) developed an FSI solver by coupling the curvilinear immersed boundary method and a free-rotation finite element model to solve systems undergoing arbitrarily large deformation. One of the case studies they considered was an inverted flag of aspect ratio $AR = H/L = 1$ and $\mu \approx 0.31$ flapping in a high Reynolds number axial flow ($Re \sim 10^5$). They obtained excellent agreement with the experimental results reported by Kim et al. (2013), particularly on the time history of the flag free-end displacement and the maximum and mean values of the drag coefficient. Tang et al. (2015) performed computational studies on the dynamics of 3D inverted flags in the $Re = 100-500$ flow regime. In particular, they examined the effects of the aspect ratio of the flag and the incidence angle of the oncoming flow. For zero incidence angle, they found that, as the dimensionless

<footnote>† Throughout the present paper, including the Introduction, the mass ratio $\mu$ is defined as fluid mass to structure mass. Some other authors use the inverse definition.</footnote>
<figure><img src="image_2.png" /></figure>

Figure 1: Two-dimensional inverted flag in axial flow: free at the upstream and clamped at the downstream end. The flag is assumed to be infinitely wide (i.e. $H/L \to \infty$); hence, the fluid flow surrounding the flag is treated as two-dimensional.

flow velocity is increased (or the bending rigidity parameter is decreased) beyond the deflected mode, the flag may undergo flapping motion about the deflected shape with the free end pointing downstream. They called this behaviour the “deflected-flapping” mode. They also found that, as the mass ratio or aspect ratio is increased, the flag starts flapping at a lower dimensionless flow velocity, which also persists for a wider range of flow velocities. Bistability and hysteresis are present in the range between the stretched-straight and flapping modes and between the flapping and deflected modes. For small non-zero incidence angles, they found that the flapping motion becomes asymmetric about the static equilibrium.

Shoele & Mittal (2016a) investigated interactions between a 2D inverted piezo-electric flag and axial flow. The majority of their simulations were conducted at $Re = 200$. They found that the inverted flag is dynamically nearly insensitive to the mass ratio for the range of dimensionless flow velocities where low- or high-amplitude symmetric flapping occurs. This finding agrees with the observations made by Kim et al. (2013), but contradicts the computational findings by Tang et al. (2015). However, it was shown that at sufficiently high flow velocities the flag bends completely, pointing downstream, and flapping becomes asymmetric. The dynamical behaviour shows a strong correlation with the mass ratio in this flow regime. This latter configuration, in fact, behaves very much like a conventional flag in axial flow – where the fluid flows from the clamped end towards the free end – whose dynamical behaviour is known to be strongly dependent on the mass ratio; see, for example, Kornecki et al. (1976), Tang & Païdoussis (2007), Connell & Yue (2007), Eloy et al. (2008) and Howell et al. (2009). Finally, they attributed the mechanism of flapping of light flags to vortex-induced vibration (VIV) and the lock-in phenomenon.

Yu et al. (2017) conducted experiments with an inverted flag in a relatively confined water channel ($W/L = 2$, $W$ being the width of the test-section) in the $Re = 4660 - 7600$ flow regime. The inverted flag had an aspect ratio $AR = 3$ and mass ratio $\mu \simeq 143$. The initially straight flag buckled in the $Re < 4940$ flow regime, becoming slightly deflected ($A/L < 0.3$). Beyond the buckling regime, while increasing (decreasing) the Reynolds number (the dimensionless bending stiffness), three other modes were encountered: the biased mode, the flapping mode, and the deflected mode. The biased mode corresponds to moderate-amplitude flapping ($0.3 \le A/L \le 1.2$) about the buckled shape in the $5100 \le Re \le 6000$ flow regime. The other two modes correspond, respectively, to large-amplitude flapping ($A/L \sim 1.8$) in the $6000 < Re \le 7122$ flow regime, and large-amplitude buckling beyond $Re = 7400$. The flapping motion was found to undergo
a hysteresis loop when the Reynolds number was decreased from $Re = 7600$ down to $Re = 4660$. The existence of large-amplitude flapping and deflected modes, as well as the hysteresis phenomenon, had previously been reported by Kim et al. (2013).

More recently, Gurugubelli & Jaiman (2019) computationally explored the mechanism for large-amplitude flapping of 3D inverted flags, assuming spanwise periodicity while ignoring the side-edge effects. They employed the arbitrary Lagrangian-Eulerian technique which allows the deforming mesh used for solving the fluid system to adapt to the deforming body's Lagrangian mesh. Most simulations were conducted for a flag with $\mu = 1$ and at $Re = 3 \times 10^4$. In particular, they showed that flapping does not arise below a critical Re, presumably because vortex shedding does not occur in the low-Re flow regime ($Re \lesssim 40$; refer to Blevins 1990 and Païdoussis et al. 2010). They found that the flag undergoes divergence leading to large static deformation instead, in that flow regime. However, most recently, Goza et al. (2018) showed that for sufficiently heavy flags, i.e. $\mu \ll 1$, large-amplitude flapping arises even for $Re < 50$, but the mechanism of instability is not classical VIV. Gurugubelli & Jaiman (2019) also found that the aspect ratio of the flag ($0.25 \le \mathcal{A}R \le 2$) has a marginal effect on large-amplitude flapping. Finally, they attributed the periodic flapping motion to a complex interplay between the unsteady shedding of leading-edge vortices and the structural dynamics of the flexible flag; the trailing-edge vortices were found to have an insignificant effect on the flapping.

In an elaborate study, Sader et al. (2016a) employed a combination of mathematical theory, scaling analysis and experimental measurement to understand the underlying physical mechanisms of flapping of inverted flags. They concluded that the large-amplitude flapping motion of an inverted flag is a vortex-induced vibration, in contrast to flapping of a conventional flag, which is a ‘self-excited’ oscillation. They provided some observations and analyses to support this hypothesis, e.g. that flapping occurred over the approximate range of Strouhal numbers $St \in [0.1 \ 0.2]$ (the Strouhal number being independent of the Reynolds number) and flapping occurred above a minimum amplitude, accompanied by flow separation and periodic vortex shedding. The Strouhal number was calculated by measuring the peak-to-peak amplitude of flapping and considering the synchronization of vortex shedding and flag motion, as observed in experiments (refer also to Kim et al. 2013). Adopting steady-state aerodynamic theory in the equations obtained by Kornecki et al. (1976), they showed that the onset of flapping is due to a static divergence instability and is thus independent of the mass ratio. In addition, they showed analytically (using lifting-line theory and the vortex-lattice method) and experimentally that the critical dimensionless flow velocity at which the zero-deflection equilibrium loses stability decreases as the aspect ratio is increased.

In a subsequent paper, Sader et al. (2016b) examined theoretically the stability of low aspect ratio (or slender) flags and cylindrical beams (or ‘rods’ as they call them) in uniform steady axial flow. Employing the Bollay (1939) quadratic normal force coefficient (i.e. $C_N = 2 \sin^2 \alpha$ for rigid, rectangular, very low aspect ratio plates, $\alpha$ being the incidence angle, or angle of attack), they developed a theory to examine the static stability of slender flexible and flexibly supported rigid inverted flags. They found that at low flow velocities the zero-deflection equilibrium is the only possible stable state. However, at a critical flow velocity, a saddle-node bifurcation occurs, from which a new deflected equilibrium state emanates and grows in amplitude as the flow velocity is increased. Their experimental measurements for low aspect ratio (down to $\mathcal{A}R = 0.033$) inverted flags agree well with their theoretical results, both qualitatively and quantitatively.

Sader et al. (2016b) also extended the theory to the case of inverted flexible cylinders and produced some results by employing experimental system parameters in Rinaldi & Païdoussis (2012). The normal force coefficient was taken as $C_D = 1.1$ (Taylor 1952; Hoerner 1965) for the theoretical calculations. They found that an inverted cylinder is stable at low flow velocities, but undergoes a saddle-node bifurcation at higher flow velocities, leading to a large-amplitude deflected equilibrium state, similarly to low aspect ratio inverted flags. This differs from Rinaldi & Païdoussis’s observations of a ‘weak’ flutter-like instability at low flow velocities, followed by
divergence at higher flow velocities. Moreover, the transition from zero-deflection equilibrium to large-amplitude deflected equilibrium should normally occur abruptly, which is again different from Rinaldi & Païdoussis's observations of a steadily growing bow with increasing flow. On the other hand, the very recent nonlinear model developed by Abdelbaki et al. (2018) for inverted cylinders in confined axial flow corroborates qualitatively several aspects of the dynamics observed by Rinaldi & Païdoussis (2012). All of the above call for fresh experimental studies with inverted flexible cylinders. Interestingly, the theoretical predictions by Sader et al. (2016b) for the saddle-node bifurcation point agree well with Rinaldi & Païdoussis's experimental measurements of the critical flow velocity for static divergence.

In a series of studies, Tavallaeinejad et al. (2020, 2018) developed nonlinear theoretical models to explore the static and dynamic responses of flexible, low aspect ratio (mostly $\mathcal{R} < 1$) inverted flags in axial flow. They used a Hamiltonian framework to derive the equations of motion. The ‘reactive’ fluid force – the fluid flow momentum change in the transverse direction due to the immersed flag’s motion – was formulated using the large-amplitude elongated body theory of Lighthill (1971), while the nonlinear or ‘vortex’ lift – due to strong vortex formation along the lateral edges of the flag (Hoerner & Borst 1985) – was approximated by the nonlinear wing theory proposed by Bollay (1939). They cast Bollay’s implicit equation for the normal force coefficient of a finite aspect ratio plate-like wing into a closed-form equation similar to that given by Polhamus (1966).

The nonlinear static response of the system was examined in Tavallaeinejad et al. (2018) by considering steady fluid-related forces. It was found that a low aspect ratio flag undergoes a static divergence via a subcritical pitchfork bifurcation followed by a saddle-node bifurcation. However, beyond a critical aspect ratio, the subcritical bifurcation was replaced by a supercritical bifurcation. Also, increasing the aspect ratio was found to have a destabilising effect, which is consistent with the experimental/analytical observations of Sader et al. (2016a) and the computational results of Tang et al. (2015). For initially curved flags (i.e. imperfect flags), the pitchfork bifurcation was replaced by a saddle-node bifurcation. By increasing the aspect ratio, the bistable range between the two saddle-node bifurcations was progressively diminished and eventually a supercritical response appeared.

Tavallaeinejad et al. (2020) explored the nonlinear dynamics of low aspect ratio ($\mathcal{R} \le 0.5$), relatively heavy ($\mu \le 0.4$) inverted flags in axial flow. For a not too slender flag, they found the following typical sequence of dynamical states, as the flow velocity is increased: (i) stable at the position of rest, (ii) flapping about the position of rest via a supercritical Hopf bifurcation (i.e. symmetric flapping mode), (iii) flapping about a deflected equilibrium (i.e. asymmetric flapping or deflected-flapping mode) via two saddle-node bifurcations, and (iv) highly deflected state on one side via a Hopf bifurcation. They also found that as $\mathcal{R}$ was increased from $\mathcal{R} = 0.1$ to $\mathcal{R} = 0.5$, the behaviour changed from a fully static response (as reported in Tavallaeinejad et al. 2018) to the complex sequence of dynamical behaviour outlined above. In addition, it was shown that the fluid flow with a non-zero incidence angle would break the symmetry of the response from essentially zero flow. Increasing the incidence angle was found to change the dynamics dramatically. At a sufficiently large incidence angle, the bistable regime vanishes altogether. Increasing the mass ratio from $\mu = 0.1$ to $\mu = 0.4$ was found to lower the critical flow velocity for flapping and to increase the flapping amplitude. The latter differs from the experimental observations by Kim et al. (2013) (with only two data points, however) but agrees with the computational results of Tang et al. (2015). Finally, increasing the skin friction coefficient, approximated by boundary layer solutions (Anderson 2017), was found to lower the critical flow velocities for the secondary bifurcations, as well as the flapping amplitude.

Finally, inspired by the planform of tree leaves, Fan et al. (2019) explored the effects of the plan-view shape of the inverted flag on its stability and dynamics. Wind tunnel tests were carried out with trapezoidal sheets of $\mathcal{R} = 1.25$, taper ratio $0 \le \Delta/H \le 1.75$ ($\Delta$ and $H$ being the height
<table><caption>Table 1: Taxonomy of recent studies on inverted flags. IBM and ALE are the abbreviations for ‘immersed boundary method’ and ‘arbitrary Lagrangian-Eulerian’ scheme, respectively; N/A stands for “not applicable”.</caption><thead><tr><th>Study</th><th>Approach</th><th>R</th><th>μ</th><th>Re</th></tr></thead><tbody><tr><td rowspan="2">Kim et al. (2013)</td><td>Experimental (air)</td><td>1.0 – 1.3</td><td>0.303 – 0.4</td><td>~ 10<sup>4</sup> – 10<sup>5</sup></td></tr><tr><td>Experimental (water)</td><td>1.3 – 2.0</td><td>166.7 – 250</td><td>~ 10<sup>4</sup> – 10<sup>5</sup></td></tr><tr><td>Ryu et al. (2015)</td><td>Computational (IBM)</td><td>→ ∞</td><td>—</td><td>≤ 250</td></tr><tr><td>Gilmanov et al. (2015)</td><td>Computational (IBM)</td><td>1.0</td><td>0.31</td><td>~ 10<sup>5</sup></td></tr><tr><td>Gurugubelli & Jaiman (2015)</td><td>Computational (ALE)</td><td>→ ∞</td><td>0.1 – 10</td><td>1000</td></tr><tr><td>Tang et al. (2015)</td><td>Computational (IBM)</td><td>0.5 – 2.0</td><td>0.2 – 1.0</td><td>100 – 500</td></tr><tr><td>Shoele & Mittal (2016a)</td><td>Computational (IBM)</td><td>→ ∞</td><td>0.1 – 10</td><td>25 – 800</td></tr><tr><td rowspan="2">Sader et al. (2016a)</td><td>Experimental (air)</td><td>0.1-7.5</td><td>~ 0.38</td><td>~ 10<sup>4</sup> – 10<sup>5</sup></td></tr><tr><td>Analytical</td><td>≥ 1</td><td>steady flow</td><td>inviscid flow</td></tr><tr><td rowspan="2">Sader et al. (2016b)</td><td>Experimental (air)</td><td>0.033 – 0.13</td><td>~ 0.38</td><td>~ 10<sup>4</sup> – 10<sup>5</sup></td></tr><tr><td>Analytical</td><td>→ 0</td><td>steady flow</td><td>inviscid flow</td></tr><tr><td>Yu et al. (2017)</td><td>Experimental (water)</td><td>3.0</td><td>143</td><td>~ 4 – 7 × 10<sup>3</sup></td></tr><tr><td>Goza et al. (2018)</td><td>Computational (IBM)</td><td>→ ∞</td><td>0.02-20</td><td>20, 200</td></tr><tr><td>Gurugubelli & Jaiman (2019)</td><td>Computational (ALE)</td><td>0.25 – 2.0</td><td>1.0</td><td>0.1 – 50, 3 × 10<sup>4</sup></td></tr><tr><td>Tavallaeinejad et al. (2018)</td><td>Analytical</td><td>0.067 – 1.7</td><td>steady flow</td><td>inviscid flow</td></tr><tr><td>Tavallaeinejad et al. (2020)</td><td>Analytical</td><td>0.1 – 0.5</td><td>0.1 – 0.4</td><td>N/A</td></tr><tr><td rowspan="2">Fan et al. (2019)</td><td>Experimental (air)</td><td>1.25</td><td>0.4, 0.61</td><td>~ 10<sup>4</sup></td></tr><tr><td>Analytical</td><td>≥ 1</td><td>steady flow</td><td>inviscid flow</td></tr><tr><td>Present paper</td><td>Analytical</td><td>→ ∞</td><td>0.1 – 5.0</td><td>inviscid flow</td></tr></tbody></table>

at the free and clamped ends, respectively) and two different thicknesses, yielding $\mu = 0.40$ and $\mu = 0.61$ (presumably for the rectangular sheets), in $Re \approx 10^4 - 10^5$ flow. They found that as the taper ratio (or shape parameter) $\Delta/H$ is increased, the critical flow velocities for both the onset and cessation of flapping and the flow speed range over which flapping occurs decrease. Interestingly, the ratio of the two critical flow velocities (i.e. cessation to onset) for all trapezoidal shapes tested was found to be about 4. In addition, they developed an extended form of the analytical model developed in Sader et al. (2016a) for predicting the divergence instability of an inverted flag of arbitrary shape, and good agreement with experimental results was achieved. Moreover, they studied experimentally and analytically a new configuration of inverted flags, where a rectangular sheet is attached to a slender elastic strip (to mimic a ‘petiole’ or the stalk that attaches the leaf to the stem) which is clamped at its other end. Similar global dynamical behaviour as for the flag with no petiole was observed, namely: undeflected equilibrium, flapping and deflected equilibrium. It was also reported that it is the petiole which provides the primary source of elasticity, while the elastic sheet undergoes rigid-body motion.

A summary and taxonomy of previous studies is presented in Table 1.†

The aim of the present paper is to investigate analytically the nonlinear dynamics and post-critical behaviour of inverted flags in axial flow. In particular, we aim at deriving an analytical model for the dynamics of high aspect ratio, heavy inverted flags by leveraging on inviscid flow

<footnote>† This table strictly aims at reporting the range of system parameters investigated in each study and does not assess the validity of their findings. The present paper explores the effect of mass ratio for $\mu \in [0.1, 5.0]$, but the results may not match all features present in the experiments of Kim et al. (2013).</footnote>
aerodynamic theory. Despite recent advances in computational fluid dynamics (CFD) methods, which make them practically viable for use in fully coupled fluid-structure interaction solvers, analytical approaches are still attractive for solving complex FSI problems. Analytical solutions are obtained in a faster fashion and thus are superior for parametric studies to, for example, computational models based on direct numerical simulation (DNS) or large-eddy simulation (LES) of the fluid flow. Moreover, analytical models often give valuable insight into the underlying mechanisms of the dynamics of the system.

Thus, a continuum *analytical* representation of fluid forces is obtained and detailed based on the quasi-steady version of the unsteady aerodynamic theory of Theodorsen (1949) to formulate the reactive force, and the leading-edge suction analogy of Polhamus (1966) to deal with the separated flow at the flag leading edge, with no account taken of periodicity of fluid dynamic forces due to leading-edge vortex shedding. The analytical model will be used to predict, for example, the onset of instabilities, frequency, and amplitude of flapping, as well as the physical mechanisms involved in the transition from one regime to another.

The rest of the paper is organized as follows. In § 2, the extended Hamilton’s principle is employed to derive the equations of motion from the kinetic and potential energies of the system, as well as the virtual work due to fluid-related forces. The resulting partial-integro-differential equations are then discretized spatially by means of the Galerkin method and they are recast as a set of nonlinear ordinary differential equations (ODEs). In § 3, the linearised form of the equations of motion is obtained, and the linear dynamics of the system is explored by solving an eigenvalue problem. The nonlinear response of the system is examined in detail in § 4, by utilizing direct time integration and bifurcation analysis. The sensitivity of the dynamics to the mass ratio and the angle of attack is also explored. The results are summarised in the form of bifurcation diagrams to identify the transition between various regimes. It is shown that lighter flags are more likely to exhibit hysteresis in their dynamical response. In § 5, the effectiveness of the model is demonstrated by comparing the obtained results with those available in the literature, highlighting that the proposed model captures both the qualitative behaviour and the detailed dynamics observed in previous studies.

## 2. Analytical modelling

In this section, a nonlinear model is developed for the inverted flag problem. First, the inverted flag problem is defined mathematically. The kinetic and potential energies of the inverted flag, idealized as a cantilevered elastic beam as shown in figure 2, are derived according to the nonlinear Euler-Bernoulli beam theory. Next, expressions for fluid forces are formulated separately, and the associated virtual work is derived. Finally, the nonlinear integro-partial-differential equations governing the rotation of the mid-plane are derived via the extended Hamilton’s principle. These equations are non-dimensionalised and then discretised via the Galerkin method.

### 2.1. *Definitions and preliminaries*

The system under consideration is shown schematically in figure 1, consisting of a vertical cantilevered thin plate subjected to an inviscid axial flow impinging on its free end. The flag is assumed to be infinitely wide (i.e. $\mathcal{R} \to \infty$); hence, the fluid flow surrounding the flag is treated as being two-dimensional. The spanwise variation of the fluid-related loads is neglected; the flag is therefore idealized as a cantilevered Euler-Bernoulli beam. The inverted flag is of infinite span, finite chord $L$, thickness $h$, and flexural rigidity $D$. The mechanical properties of the flag, i.e., the mass density, Poisson ratio, Young’s modulus, and internal damping coefficient are represented by $\rho_p$, $v_p$, $E$, and $\eta$, respectively; $\rho_f$ stands for the density of the fluid flowing with mean flow velocity $U$. Two coordinate systems are defined: (i) a right-handed rectangular Cartesian reference system ($O; XZ$), with the $X(\mathbf{e}_X)$ and $Z(\mathbf{e}_Z)$ axes being in the axial and transverse directions,
<figure><img src="image_3.png" /></figure>

Figure 2: (Left) The flag is idealised as an inextensible (i.e. $e = 0$) cantilevered beam. An arbitrary infinitesimal element $dx$ of the beam undergoes longitudinal $u(X,t)$ and transverse $w(X,t)$ displacements, related to each other by the rotation angle $\psi(X,t)$. (Right) Joukowski's conformal mapping of the plate into a circle of radius $L/4$ (i.e. $\bar{X}\bar{Z}$-plane); $q_r$ and $q_{\theta}$ are the radial and circumferential velocities on the circle in the $\bar{X}\bar{Z}$-plane, respectively. The impermeability condition is satisfied by considering a set of sources and sinks of strength $\sigma$ on the top and bottom halves of the circle, respectively. The induced velocity due to a single bound vortex of strength $\Gamma$ and its image (i.e. the wake vortex of strength $-\Gamma$) at an arbitrary point on the circle satisfies the Kutta condition.

respectively; (ii) a curvilinear coordinate system ($O; xZ$), with the $x$-axis being along the flag centreline from its clamped end towards the free one. The unit vectors $e_x$ and $e_z$ are tangential and normal to the centreline, respectively.

The rotation angle $\psi(X,t)$, the curvature $\kappa(X,t)$, and the axial strain $e(X,t)$ are related to the longitudinal and transverse motions, $u(X,t)$ and $w(X,t)$ respectively, of a generic point at a distance $Z$ from the mid-plane on the cross-section by

$$ \sin \psi(X,t) = \frac{\partial_X w(X,t)}{1 + e(X,t)}, \quad \cos \psi(X,t) = \frac{1 + \partial_X u(X,t)}{1 + e(X,t)}, \quad \kappa(X,t) = \frac{\partial_X \psi(X,t)}{1 + e(X,t)}, \quad (2.1) $$

where $\partial_X \equiv \partial/\partial X$ denotes the first spatial derivative.† The inextensibility assumption reads $e(X,t) = 0$, which yields $(1 + \partial_X u(X,t))^2 + (\partial_X w(X,t))^2 = 1$ (Farokhi et al. 2018), indicating that the length of the flag centreline remains constant during vibrations. This condition allows us to express all physical quantities in terms of the curvilinear coordinate $(x,t)$ (refer to Nayfeh & Pai 2008 or Lopes et al. 2002). Additionally, it reduces the dependent variables to one, and $w(x,t)$ and $u(x,t)$ may be expressed in terms of the rotation angle of the cross-section $\psi(x,t)$. Hence, velocity components of an element may be formulated as

$$ \begin{aligned} \dot{w}(x,t) &= \int_{0}^{s} \dot{\psi}(s,t) \cos \psi(s,t) \, ds, \\ \dot{u}(x,t) &= -\int_{0}^{s} \dot{\psi}(s,t) \sin \psi(s,t) \, ds, \end{aligned} \quad (2.2) $$

where the overdot denotes the derivative with respect to time.

Considering a set of longitudinal and transverse displacements (i.e. $u(x,t)$ and $w(x,t)$, respectively) of the element $dx$ shown in figure 2, away from its undeflected static equilibrium, the relative velocity between the structure and the incident flow in the tangential and normal

<footnote>† In the remainder of this paper, $\partial_\bullet$ is defined as the first derivative with respect to $\bullet$, $\partial_{\bullet\bullet}$ denotes the second derivative, and so on.</footnote>
directions may be written as

$$ 
\begin{aligned}
\mathbf{V}_{\text{rel}}(x, t) &= V_{\tau}(x, t)\mathbf{e}_x + V_{\text{n}}(x, t)\mathbf{e}_z \\
&= [(\dot{u} + U) \cos \psi + \dot{w} \sin \psi]\mathbf{e}_x + [-(\dot{u} + U) \sin \psi + \dot{w} \cos \psi]\mathbf{e}_z.
\end{aligned}
\quad (2.3) $$

In this paper, the governing equations are derived in a Hamiltonian framework. The extended Hamilton's principle may be expressed as

$$ \int_{t_1}^{t_2} (\delta\mathcal{T}(t) - \delta\mathcal{V}(t) + \delta\mathcal{W}_f(t) + \delta\mathcal{W}_d(t)) dt = 0, \quad (2.4) $$

with $\mathcal{T}(t)$ and $\mathcal{V}(t)$ denoting the kinetic and potential energies of the system, respectively, and $\delta\mathcal{W}_f(t)$ and $\delta\mathcal{W}_d(t)$ the virtual work by the fluid forces acting on the inverted flag and damping, respectively; here $\delta$ denotes the variational operator.

## 2.2. Kinetic and potential energies

In this study, the centreline rotation of the inverted flag $\psi(x, t)$, is taken as the primary variable for describing the flag motion. This can be achieved with the aid of equations (2.1) and taking into account the inextensibility assumption. Derivation of the equation of motion in terms of the rotation angle allows predicting large-amplitude deformations even when the tip rotation exceeds $\pi/2$ (refer to Tavallaeinejad et al. 2018). Thus, in what follows, the kinetic and potential energies as well as the virtual work of fluid-related forces and damping are formulated as functions of $\psi(x, t)$.

The kinetic energy of the system for translational and rotational motions of the inverted flag may be written as

$$ 
\begin{aligned}
\mathcal{T}(t) = & \frac{1}{2} \rho_{\text{p}} h \int_{0}^{L} \left[ \left( \int_{0}^{x} \sin \psi(s, t) \dot{\psi}(s, t) \, ds \right)^{2} + \left( \int_{0}^{x} \cos \psi(s, t) \dot{\psi}(s, t) \, ds \right)^{2} \right] dx \\
& + \frac{1}{24} \rho_{\text{p}} h^{3} \int_{0}^{L} \dot{\psi}(x, t)^{2} \, dx.
\end{aligned}
\quad (2.5) $$

The potential strain energy of the inextensible inverted flag in terms of the rotation angle is

$$ \mathcal{V}(t) = \frac{1}{2} D \int_{0}^{L} (\partial_{x} \psi(x, t))^{2} dx, \quad (2.6) $$

in which $D = Eh^3/(12(1 - v_p^2))$ is the plane-strain flexural rigidity for plates of large aspect ratios. The virtual work done by the Kelvin–Voigt structural damping† may be formulated as

$$ \delta W_d(t) = -D\eta \int_0^L \partial_{xx} \dot{\psi}(x, t) \delta\psi \, dx, \quad (2.7) $$

where $\eta$ denotes the viscoelastic damping coefficient. The virtual work done by the fluid-related force normal to the flag is given by

$$ \delta W_F(t) = \int_{0}^{L} (F_N(x, t) \cos \psi(x, t) \delta w - F_N(x, t) \sin \psi(x, t) \delta u) \, dx, \quad (2.8) $$

in which $F_N$ represents the normal force acting on the flag, and is derived in the following section.

† Material damping is introduced in the model since it is an inherent property of all materials; it provides a way for modelling the dissipation of the energy pumped into the flag from the fluid flow. The structural damping model utilized in the proposed model is more appropriate than a generic linear viscous damping model.
## 2.3. Aerodynamic model

The inviscid pressure-related forces are modelled in this section. The flow is assumed to be inviscid and two-dimensional around a thin flat plate. The 2D quasi-steady thin airfoil theory, involving large angles of attack, is employed to formulate the reactive fluid forces corresponding to the relative motion of the inverted flag with respect to the incident flow; for details the reader is referred to Bisplinghoff et al. (2013), Ramesh et al. (2013) and Yan et al. (2014). The influence of leading edge separation is dealt with separately. First, the quasi-steady forces are derived using a velocity potential approach in a way that the solution is divided into non-circulatory and circulatory contributions; each part is obtained using a Joukowski conformal transformation. Next, the separated flow at the leading edge is modelled utilizing the leading edge suction analogy of Polhamus (1966).

We consider a two-dimensional perturbation velocity potential due to the motion of the plate $\phi(x, z, t)$, which satisfies the Laplace equation (i.e. $\nabla^2\phi = 0$, with $\nabla$ being a first order spatial differential operator). Following Theodorsen's approach, the perturbation velocity potential may be divided into two parts: $\phi = \phi^{\text{nc}} + \phi^{\text{c}}$, where $\phi^{\text{c}}$ is the circulatory component, fulfilling the Kutta condition, and $\phi^{\text{nc}}$ denotes the non-circulatory component which satisfies the impermeability boundary condition.

The pressure difference between the lower and upper surfaces may be formulated using the unsteady Bernoulli equation (refer to Ramesh et al. 2013 and Yan et al. 2014),

$$ \Delta P(x,t) = P_u(x,t) - P_l(x,t) = -2\rho[\dot{\phi}_u - ((\dot{u} + U)\cos\psi + \dot{w}\sin\psi)\partial_x\phi_u]. \quad (2.9) $$

Following Theodorsen's approach, the pressure difference over the flag is decomposed into two contributions, namely non-circulatory and circulatory components. Hence, the total pressure may be written as

$$ \Delta P(x,t) = \Delta P^{\text{nc}}(x,t) + \Delta P^{\text{c}}(x,t). \quad (2.10) $$

The non-circulatory (or added mass) component, represented by $\Delta P^{\text{nc}}(x,t)$, contributes to the perturbation pressure due to the displacement of the surrounding fluid (Bisplinghoff et al. 2013). The circulatory contribution (vortex-induced pressure), denoted by $\Delta P^{\text{c}}(x,t)$, on the other hand, dominates the perturbation pressure in the vicinity of the trailing edge due to shedding of vorticity (Staubli & Rockwell 1989).

### 2.3.1. Non-circulatory contribution

In order to find the solution of the two-dimensional incompressible potential flow, the plate is mapped into a circle with its centre at the origin and radius $r = L/4$, through the following standard conformal mapping illustrated in figure 2:

$$ x + iz = \bar{X} + \frac{L}{2} + i\bar{Z} + \frac{L^2}{16(\bar{X} + i\bar{Z})}. \quad (2.11) $$

Given that $\bar{X} = L/(4 \cos \theta)$ and $\bar{Z} = L/(4 \sin \theta)$, the relation between points on the flag and the circle may be found using (2.11) as $x = (1 + \cos \theta)/2$.

Considering a continuous distribution of source-sink of strength $\sigma(\theta, t)$ and $-\sigma(\theta, t)$ over the upper and lower surfaces of the sheet (or on its circle counterpart in figure 2), which satisfies the impermeability boundary condition, and after several manipulations, the strength per unit length of the source sheet can be found as $\sigma(\theta, t) = 4V_n(\theta, t) \sin\theta$ (refer to Bisplinghoff et al. 2013). Considering an infinitesimal element of the source sheet (L/4) $d\varphi$ on the circle at $\theta = \varphi$, along with its sink counterpart at $\theta = -\varphi$, and integrating over the half circle, the non-circulatory
component of the resultant tangential velocity may be formulated as

$$q_{\theta}^{\text{nc}}(\theta, t) = \frac{2}{\pi} \int_{0}^{\pi} V_{n}(\varphi, t) \frac{\sin^2 \varphi}{\cos \varphi - \cos \theta} d\varphi. \qquad (2.12)$$

Assuming that $\phi_u^{\text{nc}}(L/4, 0, t) = 0$ at the leading edge (i.e. $\theta = 0$), the non-circulatory perturbation velocity potential in the $\bar{X}\bar{Z}$ plane at $r = L/4$ may be found,

$$\phi_u^{\text{nc}}(\theta, t) = -\frac{L}{2\pi} \int_{\theta}^{0} \int_{0}^{\pi} \frac{V_n(\varphi, t) \sin^2 \varphi}{\cos \varphi - \cos \theta} d\varphi \, d\theta. \qquad (2.13)$$

Substitution of (2.13) into (2.9), changing the variables by making use of $x = (1 + \cos \theta)/2$, and using the chain rule, the pressure distribution over the flag due to the non-circulatory contribution of the flow in the $\bar{X}\bar{Z}$ plane at $r = L/4$ may be written as

$$\begin{aligned} \Delta P^{\text{nc}}(\theta, t) = & -\rho \left[ \frac{L}{2} \partial_t \left( \int_{\theta}^{0} \int_{0}^{\pi} \frac{V_n(\varphi, t) \sin^2 \varphi}{\cos \varphi - \cos \theta} d\varphi \, d\theta \right) \right. \\ & \left. - \frac{1}{\sin \theta} ((\dot{u} + U) \cos \psi + \dot{w} \sin \psi) \int_{0}^{\pi} \frac{V_n(\varphi, t) \sin^2 \varphi}{\cos \varphi - \cos \theta} d\varphi \right]. \end{aligned} \qquad (2.14)$$

### 2.3.2. Circulatory contribution

The Kutta condition enforces that the flow leaves smoothly the trailing edge of the inverted flag and hence avoids an infinite velocity at the trailing edge (i.e. $\theta = \pi$); see figure 3. It is known that the non-circulatory solution by itself is unable to satisfy the Kutta condition. Theodorsen (1949) fulfilled the condition by considering bound vortices, along with a wake of counter-vortices continuously shedding from the trailing edge at the freestream velocity. The induced velocity due to the bound vortices and their ‘images’ cancels the tangential velocity at the trailing edge due to the non-circulatory component of the flow. Hence,

$$q_{\theta}^{\text{nc}}(L/4, \theta = \pi, t) + q_{\theta}^{c}(L/4, \theta = \pi, t) = 0. \qquad (2.15)$$

Subsequently, Theodorsen obtained the induced velocity due to a single bound vortex of strength $\Gamma$ and its image (i.e. the wake vortex of strength $-\Gamma$) at an arbitrary point on the circle and formulated the circulatory component of the tangential velocity as

$$q_{\theta}^{c}(\theta, \bar{X}, t) = -\frac{\Gamma}{2\pi} \left( \frac{\bar{X}^2 - L^2/4}{\bar{X}^2 + L^2/4 + \bar{X}(L/2)\cos\theta} \right). \qquad (2.16)$$

Applying the Kutta condition at the trailing edge (i.e. $\theta = \pi$) and the quasi-steady assumption,† which neglects only the influence of the wake vortices on the flow (refer to (Bisplinghoff et al. 2013), the pressure distribution originating from the circulatory contribution of the flow may be written as‡

$$\Delta P^{c}(\theta, t) = -2\rho \left( -\dot{\phi}_{u}^{c}(\theta, t) + \frac{2V_{\tau}}{L \sin \theta} \partial_{\theta} \phi_{u}^{c}(\theta, t) \right) = -\frac{\rho V_{\tau} q_{\theta}^{\text{nc}}(\pi, t)}{\sin \theta}. \qquad (2.17)$$

It ought to be noted that expression (2.17) ignores the effects of flow history in the wake behind the trailing edge, and instead considers the forces generated by the instantaneous deformation, velocity, and acceleration of the flag; the rationale for this is detailed in the following.

Direct numerical simulations by Gurugubelli & Jaiman (2019) for an inverted flag with a splitter

† This is an approximation of Theodorsen’s theory, in which $C(k) \to 1$, with $C(k)$ being Theodorsen’s function based on Hankel functions of the second kind, with $k$ denoting the reduced frequency. The quasi-steady theory can be used in the time domain for slow harmonic oscillations with small $k$, or slowly varying motion that is not harmonic (Hodges & Pierce 2011).

‡ To see the full derivation, the interested reader is referred to Bisplinghoff et al. (2013).
plate at the trailing edge reveal that, despite the total suppression of vortex shedding from the trailing edge, the large-amplitude flapping persists, with only a minor change in the frequency and amplitude of the oscillations. This suggests that the large-amplitude flapping of inverted flags may not require wake-flag synchronization. In other words, the effects of wake flow unsteadiness on the dynamical behaviour of the system seems to be insignificant, suggesting that a quasi-steady modelling of fluid forces may be acceptable.

Nevertheless, future work ought to employ a nonlinear analytical model accounting for the unsteady wake function (in the time domain) behind the trailing edge, to further investigate the global dynamical behaviour of the inverted-flag system; this is not a trivial task and it is beyond the scope of this paper.

### 2.3.3. Leading-edge suction force

The 2D quasi-steady aerodynamic model employed in this paper is based on inviscid potential flow theory which, by definition, ignores any viscous effects such as flow separation or viscous stresses acting on the inverted flag. Accounting for the effects of flow separation from the leading edge is of particular importance for modelling the dynamics of inverted flags. As observed by Yu et al. (2017), the flow separates from the leading edge, giving rise to a stable attached vortex above the leading edge up to the angle of attack of about 22°, and a normal force stall at larger angles of attack where reattachment fails to occur. The leading edge flow separation reduces the pressure over the upper surface of the plate and its impact becomes typically stronger at high incidence angles on thin sharp-edged plates. In fact, potential flow theories alone tend to underestimate the normal force for thin sharp-edged plates at high angles, and they need to be modified either by employing analytical models, such as the Polhamus suction analogy, or by employing fully numerical schemes to model the actual flow field.

In the present work, in order to account for the effects of the leading edge vortices, the nonlinear vortex force associated with the separated flow from the leading edge is formulated via the Polhamus analogy (Polhamus 1966) as a correction to the quasi-steady pressure derived earlier. This model formulates a nonlinear relation between the potential force and the angle of attack, which acts normal to the flag, and thus modifies the normal force, as shown in figure 3. Thus, the effect of leading edge vortices can be taken into account via superimposing linearly the nonlinear vortex lift and the quasi-steady forces due to the attached flow (Sane 2003).†

In order to implement the method outlined above, the leading edge suction needs to be expressed in terms of the quasi-steady force obtained earlier. Hence, the Blasius formula (Yan et al. 2014) is employed to determine the leading edge suction. Integration of the pressure around a contour at the leading edge (LE), and application of Bernoulli's equation yields

$$F_x - iF_z = -\frac{i\rho L}{4} \int_{LE} \left( \frac{q_{\bar{X}} - iq_{\bar{Z}}}{1 - 2e^{-i\theta}} \right)^2 \sin\theta \, d\theta, \quad (2.18)$$

where $F_x$ and $F_z$ represent the aerodynamic forces acting on an immersed body. Assuming that the contribution of the non-circulatory component in the induced velocity near the leading edge outweighs the other contributions and recalling that $q_{\bar{X}} = -q_{\theta}^{nc} \sin \theta$ and $q_{\bar{Z}} = q_{\theta}^{nc} \cos \theta$, the leading edge suction force at $r = L/4$ is

$$F_p(t) = \frac{\pi \rho L}{16} (q_{\theta}^{nc}(0, t))^2. \quad (2.19)$$

Although the suction analogy does not model the actual flow field, it has been shown to be

<footnote>† A more comprehensive approach may be to determine all contributions together, via a direct utilization of the Navier-Stokes equations, for instance. The approach employed here simplifies the analysis considerably and has been shown to give acceptable results for a wide range of problems; for example, refer to the work of Pedersen & Żbikowski (2006) on insect-like flapping wings in hover and Azuma & Okamoto (2005) on deflected wings, among others.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Flow around the 2D inverted flag placed in an inviscid fluid. The Kutta condition enforces a smooth, tangential flow leaving the trailing edge. However, according to Polhamus's leading edge suction analogy, the leading edge suction vanishes due to separation at the leading edge, and it forms a leading edge vortex. This results in a normal force component of equal magnitude to the leading edge suction. Effectively, the leading edge suction force is rotated by 90° onto the low-pressure side of the flag.</figcaption></figure>

remarkably accurate, especially at high angles of incidence. For instance, this model has been used successfully in predicting the nonlinear vortex lift of delta wings (Purvis 1981, DeVoria & Mohseni 2017) and in the theoretical modelling of flapping flight of insects (Pedersen & Żbikowski 2006).

It ought to be emphasized that the leading-edge suction model used in the present paper fails to include the effects of vortex-shedding from the leading edge. Thus, the time-variant fluid forces acting on the flag model may not be truly representative of those in reality, which in turn may result in inaccurate predictions of the dynamical behaviour (stability, amplitude and frequency of oscillations). Nevertheless, how the vortex shedding load affects both the local and global dynamics of the flag is not yet known and is yet to be understood. Accounting for such a load in the analytical model is a very challenging task and is not attempted in this paper.

### 2.3.4. Total aerodynamic force

As mentioned in the last section, the nonlinear vortex lift can be superimposed linearly on the quasi-steady pressure. Hence, the total quasi-steady normal force acting on the inverted flag is obtained as

$$F_N(\theta, t) = -\Delta P(\theta, t) - 2F_p(t) \delta_D(\cos \theta - 1), \quad (2.20)$$

with the total quasi-steady pressure formulated as

$$\Delta P(\theta, t) = \rho \left[ \frac{L}{2} \partial_t \left( \int_0^\theta q_\theta^{nc}(\vartheta, t) \, d\vartheta \right) + \frac{V_\tau}{\sin \theta} (q_\theta^{nc}(\theta, t) - q_\theta^{nc}(\pi, t)) \right] \quad (2.21)$$

where $\delta_D(x - L)$ denotes the Dirac delta function. It is noted that dealing with a non-homogenous, non-classical boundary condition is a more challenging task; hence, the end shear due to the force that results from the Polhamus effect at the leading edge is transferred from the boundary condition into the equation of motion via the Dirac delta function.

Substituting (2.21) into (2.8) and recalling that $x = (1 + \cos \theta)/2$ yields the virtual work done by the fluid-related force normal to the flag.
## 2.4. Equation of motion

Substituting equations. (2.5)-(2.8) into the generalized Hamilton's principle given in (2.4) and performing integration by parts yields the nonlinear governing equation of motion in terms of the rotation angle $\psi(x, t)$:

$$ \frac{1}{12} \rho_p h^3 \ddot{\psi}(x,t) - D \partial_{xx}(1 + \eta \partial_t) \psi(x,t) \\ + \rho_p h \left[ \cos \psi(x,t) \partial_{tt} \left( \int_x^L \int_0^s \sin \psi(\xi, t) \, d\xi \, ds \right) - \sin \psi(x,t) \partial_{tt} \left( \int_x^L \int_0^s \cos \psi(\xi, t) \, d\xi \, ds \right) \right] \\ + \cos \psi(x,t) \int_L^x F_N(s,t) \cos \psi(s,t) \, ds + \sin \psi(x,t) \int_L^x F_N(s,t) \sin \psi(s,t) \, ds = 0. \quad (2.22) $$

Utilizing the dimensionless parameters

$$ x^* = \frac{x}{L}, \quad t^* = \frac{t}{\tau}, \quad \beta = \frac{1}{12} \left(\frac{h}{L}\right)^2, \quad \eta_s = \frac{\eta}{\tau}, \quad f^* = f\tau, \quad \mu = \frac{\rho_f L}{\rho_p h}, \quad \Pi = \sqrt{\frac{\rho_f L}{D}} LU, \quad (2.23) $$

in which $f$ is the frequency of flapping and $\tau = \sqrt{\rho_p h/D L^2}$, one obtains the following nonlinear dimensionless equation of motion, in which the asterisk notation for $x^*$ and $t^*$ has been dropped throughout for simplicity:

$$ \beta \ddot{\psi}(x,t) - (\partial_{xx}(1 + \eta_s \partial_t))\psi(x,t) \\ + \cos \psi(x,t) \partial_{tt} \left( \int_x^1 \int_0^s \sin \psi(\xi, t) \, d\xi \, ds \right) - \sin \psi(x,t) \partial_{tt} \left( \int_x^1 \int_0^s \cos \psi(\xi, t) \, d\xi \, ds \right) \\ + \cos \psi(x,t) \int_L^x C_f(s,t) \cos \psi(s,t) \, ds + \sin \psi(x,t) \int_L^x C_f(s,t) \sin \psi(s,t) \, ds = 0, \quad (2.24) $$

where

$$ C_f = \frac{\mu}{\pi} \int_{\theta}^{0} \int_{0}^{\pi} \left( (\ddot{w} - \dot{u}\dot{\psi}) \cos \psi - (\ddot{u} + \dot{w}\dot{\psi}) \sin \psi - \frac{\Pi}{\sqrt{\mu}} \dot{\psi} \cos \psi \right) f(\varphi, \vartheta) \, d\varphi \, d\vartheta \\ - \frac{2}{\pi \sin \theta} \left[ (\dot{u} \cos \psi + \dot{w} \sin \psi) \left( \int_{0}^{\pi} (\mu(\dot{w} \cos \psi - \dot{u} \sin \psi) - \Pi \sqrt{\mu} \sin \psi) g(\varphi, \theta) \, d\varphi \right) \right. \\ \left. + \Pi \cos \psi \int_{0}^{\pi} (\sqrt{\mu}(\dot{w} \cos \psi - \dot{u} \sin \psi) - \Pi \sin \psi) g(\varphi, \theta) \, d\varphi \right] - \frac{1}{16\pi} I_1 |I_1| \delta_D(s-1), $$

with

$$ I_1 = \int_0^\pi (1 + \cos \varphi)(\Pi \sqrt{\mu} \dot{w} \cos \psi - (\Pi \sqrt{\mu} \dot{u} + 1) \sin \psi) \, d\varphi $$

and

$$ f(\varphi, \vartheta) = \frac{\sin^2 \varphi}{\cos \varphi - \cos \vartheta}, \quad g(\varphi, \theta) = \frac{(1 - \cos \varphi)(1 + \cos \theta)}{\cos \varphi - \cos \theta}, $$

in which $\psi = \psi(\xi, t)$, $w = w(\xi, t)$ and $u = u(\xi, t)$ for simplicity. It is recalled that $C_f = C_f(s, t)$ via $s = (1 + \cos \theta)/2$ and $\xi = (1 + \cos \varphi)/2$. The absolute operator is used in $C_f$ to modify the force, accounting for the inherent symmetry of deformations with respect to the longitudinal axis.

The clamped free boundary conditions read $\psi(0, t) = 0$ and $\partial_x \psi(L, t) = 0$.
### 3. Linear stability analysis

In this section, typical linear dynamics of the system, the mechanism of instability, and the sensitivity of the first two critical flow velocities to the mass ratio are explored. A linearised partial differential equation is first derived in the vicinity of the position at rest and is then transformed into an ordinary differential equation (ODE) via the Galerkin method, which may subsequently be transformed to an eigenvalue problem. The stability analysis is performed by considering the eigenvalues or equivalently eigenfrequencies of the system obtained by solving the eigenvalue problem for given system parameters.

#### 3.1. Linear equation of motion

The linearised form of the nonlinear equation presented in § 2.4 is expressed in terms of the transverse displacement $w(x, t)$. The resulting equation differs from that in Sader et al. (2016a) in that the time-dependent terms in the fluid model adopted here have been retained. Assuming a small and purely transverse deflection of the inverted flag (i.e. $w(x, t) \ll 1$, $u(x, t) \approx 0$), equation (2.24) may be linearised as

$$ 
\begin{aligned}
& \ddot{w}(x,t) + (1 + \eta_s \partial_t) \partial_{xxxx} w(x,t) \\
& + \frac{1}{\pi} \int_{\theta}^{0} \int_{0}^{\pi} (\mu \ddot{w}(\xi, t) - \Pi \sqrt{\mu} \partial_{\xi} \dot{w}(\xi, t)) f(\varphi, \vartheta) \, d\varphi \, d\vartheta \\
& + \frac{2\Pi}{\pi \sin \theta} \int_{0}^{\pi} (\Pi \partial_{\xi} w(\xi, t) - \sqrt{\mu} \dot{w}(\xi, t)) g(\varphi, \theta) \, d\varphi = 0,
\end{aligned} 
\quad (3.1) 
$$

where $x = (1 + \cos \theta)/2$ and $\xi = (1 + \cos \varphi)/2$. The boundary conditions for the clamped-free flag are expressed as

$$ w(0, t) = 0, \quad \partial_x w(0, t) = 0 \quad \partial_{xx} w(1, t) = 0, \quad \partial_{xxx} w(1, t) = 0. \quad (3.2) $$

The series expansion

$$ w(x, t) = \sum_{i=1}^{M} W_i(x) p_i(t) \quad (3.3) $$

is introduced as an approximate solution for the linearised equation of motion of the inverted flag, where $p_i(t)$ are the time-dependent generalized coordinates, and $W_i(x)$ are the corresponding mode shapes of a cantilevered beam in vacuo; $M$ is the number of modes utilized in the solution.

It is noted that direct application of conventional eigenfunctions in the expansion (3.3) hinders the development of the exact evaluation of the non-circulatory velocity potential due to the presence of singular terms. In order to overcome this problem, an expansion of Chebyshev polynomials of the first kind with a sufficiently large number of terms is utilized to approximate each mode shape of the sheet. Polynomial representation of the mode shapes allows the integration in (3.1) to be performed exactly over $\varphi$ with the aid of Glauert's integral given in Eq. (A 1).

Applying the Galerkin technique, equation (3.1) reduces to the well-known second-order form $M\ddot{p} + C\dot{p} + Kp = 0$ (see Eq. (A 7)). Defining the following partitioned vector and matrices of order $2M$:

$$ r = \begin{Bmatrix} \dot{p} \\ p \end{Bmatrix}, \quad B = \begin{bmatrix} 0 & M \\ M & C \end{bmatrix}, \quad E = \begin{bmatrix} -M & 0 \\ 0 & K \end{bmatrix}, \quad (3.4) $$

the problem reduced into the first-order form $B\dot{r} + Er = 0$. Seeking periodic solutions in the form of $r(t) = \tilde{r} \exp(i\omega t)$, where $\tilde{r}$ is the modal amplitude and $\omega$ is the complex frequency, one can obtain the following eigenvalue problem $[i\omega B + E]\tilde{r} = 0$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Argand diagram for an inverted flag with $\mu = 1.0$, showing the evolution of the first three modes with the dimensionless flow velocity $\Pi$.</figcaption></figure>

<caption>Figure 5: Stability map for a 2D inverted flag; [—]: the onset of divergence, [—]: the onset of flutter.</caption>

## 3.2. Numerical results

The dynamics of the linear system is examined by obtaining the eigenfrequencies of the system as a function of the dimensionless flow velocity. The evolution of the generally complex eigenfrequencies of the system with varying flow velocity is commonly displayed as an Argand diagram. In the Argand diagram, the ordinate and the abscissa are associated with the imaginary and real parts of the eigenfrequencies, i.e. $Im(\omega)$ and $Re(\omega)$, respectively. It is recalled that $Im(\omega)$ corresponds to the damping ratio, and $Re(\omega)$ represents the dimensionless frequency of the oscillation. Instability materializes by the crossing of the eigenfrequency locus from the positive to the negative half-plane in the Argand diagram (i.e. when $Im(\omega) < 0$); refer to Païdoussis (2014) for more details.

In this study, six modes are employed for numerical solutions (i.e. $M = 6$); this is sufficient for obtaining converged solutions. The other system parameters utilized in the simulations are as follows: $L = 10$ cm, $h = 0.1$ mm, $\rho = 1200$ kg m$^{-3}$, $\rho_f = 1.2$ kg m$^{-3}$, $D = 2454$ N cm$^2$ and $\eta_s = 0.0002$.

Figure 4 is the Argand diagram for an inverted flag with $\mu = 1.0$, showing the evolution of the lowest three eigenfrequencies of the system. The linear response can be interpreted as follows.

(i) Free motion of the inverted flag is damped at low flow velocities; this is evident from the positive imaginary part of the eigenfrequencies.

(ii) At a sufficiently high flow velocity, the frequency of the first mode becomes purely imaginary, and then bifurcates on the $Im(\omega)$-axis. One of the solution branches associated with this mode eventually vanishes altogether (i.e. $Re(\omega) = 0$, $Im(\omega) = 0$) at $\Pi = 1.36$ and crosses from the positive to the negative half-plane, indicating the onset of static divergence (buckling) via a pitchfork bifurcation.

(iii) At a higher flow velocity ($\Pi = 7.58$), the negative and positive branches of the first mode loci coincide and leave the imaginary axis, while $Im(\omega) < 0$ (shown in the inset in this figure). This indicates the onset of Païdoussis-type coupled-mode flutter (Done & Simpson 1977,
Païdoussis 2014). This bifurcation emanates directly from a divergent state; hence, the frequency of oscillation is zero at the onset of flapping.

(iv) The eigenfrequency of the first mode becomes purely imaginary once again as the flow velocity is increased further, at $Π = 14.74$, giving rise to a second buckling of the inverted flag in the first mode.

(v) The frequency of the second mode diminishes, and it vanishes circa $Π = 15$.

Let us now consider how well linear theory predicts the dynamics of the system as observed in the experiments. For a physical system, the inverted flag loses its straight, undeformed configuration and buckles similarly to a column under axial load. Thereafter, the flag exhibits flapping motion. Increasing the flow velocity further, the flapping motion gradually diminishes in magnitude, and the flag displays a buckled shape once again, but with a larger amplitude of deformation (compared to the initial buckling) as the fluid dynamic forces grow with increasing flow velocity. Thus, up to this point, the linear model succeeds in predicting, qualitatively at least, the dynamics of the physical system.

It is worth mentioning that stability of the stretched-straight state of the inverted flag up to the occurrence of divergence is not affected by the mass ratio parameter. Thus, for static divergence, the mass ratio does not change the value of the critical flow velocity. This is evident from the stability map shown in figure 5, which highlights the sensitivity of the foregoing bifurcations to the mass ratio parameter. The onset of flapping, however, is sensitive to the mass ratio. More specifically, the Païdoussis-type coupled-mode flutter occurs at considerably higher values of $Π$ as $μ$ is decreased, in the range $μ ∈ [0.1, 0.6]$. For larger values of $μ$, the change becomes much less significant, and the curve of the boundary of stability displays a near-plateau.

The divergence instability mechanism was originally suggested by the Kim et al. (2013) experiments (where they argued that changing the mass ratio does not affect the onset of flapping), and by Sader et al. (2016a) theoretically via a linear stability analysis. Sader et al. (2016a) argued that the divergence instability of the inverted flag is a steady process, with steady aerodynamic forces being involved. Then, unsteadiness in the flow comes into play and leads to flapping motion. This is also predicted by the nonlinear model in the next section.

Before closing this section, an important note should be made. Since linear theory is valid only up to the onset of the first instability, in this case divergence, the system must be analysed by the nonlinear model to ascertain the reliability of the predicted post-divergence behaviour. Thus, the loss of the stability of the original equilibrium predicted by the linear model is valid. However, other states may be encountered once nonlinear effects are taken into account. The reason for this is that, in the linear model, it is assumed that motions are small, remaining in the proximity of the undeflected equilibrium. This is not the case for larger values of the flow velocity, where the inverted flag diverges considerably, away from original equilibrium.

## 4. Post-divergence behaviour

The nonlinear aspects of the system dynamics, such as limit-cycle frequencies and amplitudes, as well as the transition between the various dynamical states are studied numerically. The nonlinear integro-partial-differential model of equation (2.24) is first discretized in space, and then recast in a system of first-order ODEs. The resultant set of ODEs is solved using a pseudo-arclength continuation technique. The Gear backward differentiation formula (BDF) is also used to obtain the time histories of oscillation.

### 4.1. *Nonlinear discretised model*

The spatial discretisation of the nonlinear equation of motion (2.24) is made in a similar fashion to that described in § 3.1 by employing the Galerkin method. Hence, a linear combination of
suitable comparison functions $\Psi_i(x)$ is selected to approximate the rotation of the flag

$$ \psi(x, t) = \sum_{i=1}^{M} \Psi_i(x) q_i(t), \quad (4.1) $$

where $q_i(t)$ represent the unknown time-dependent generalized coordinates; $M$ denotes the number of modes used.

As discussed in § 3.1, Chebyshev polynomials of the first kind with a sufficient number of terms are utilised to approximate the comparison functions $\Psi_i(x)$ (see Eq. (A 8)). Applying Galerkin's technique, i.e. inserting Eq. (4.1) into Eq. (2.24), multiplying it by the corresponding comparison functions, and integrating the resulting equation with respect to the spatial coordinate from 0 to 1 results in a set of $M$ second-order ODEs with trigonometric nonlinearities:

$$
\begin{aligned}
& \beta \sum_{j=1}^{M} \left( \int_{0}^{1} \Psi_{i} \Psi_{j} \, dx \right) \ddot{q}_{j} - \eta_{s} \sum_{j=1}^{M} \left( \int_{0}^{1} \Psi_{i} \Psi_{j}'' \, dx \right) \dot{q}_{j} - \sum_{j=1}^{M} \left( \int_{0}^{1} \Psi_{i} \Psi_{j}'' \, dx \right) q_{j} \\
& \qquad + \int_{0}^{1} \Psi_{i}(x) \left[ \sin\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \int_{0}^{s} \partial_{tt} \left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) \sin\left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) d\xi \, ds \right] dx \\
& \qquad + \int_{0}^{1} \Psi_{i}(x) \left[ \cos\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \int_{0}^{s} \partial_{tt} \left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) \cos\left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) d\xi \, ds \right] dx \\
& \qquad + \int_{0}^{1} \Psi_{i}(x) \left[ \sin\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \int_{0}^{s} \left(\partial_{t} \left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right)\right)^{2} \cos\left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) d\xi \, ds \right] dx \\
& \qquad - \int_{0}^{1} \Psi_{i}(x) \left[ \cos\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \int_{0}^{s} \left(\partial_{t} \left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right)\right)^{2} \sin\left(\sum_{j=1}^{M} \Psi_{j}(\xi) q_{j}(t)\right) d\xi \, ds \right] dx \\
& \qquad + \int_{0}^{1} \Psi_{i}(x) \left[ \cos\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \Delta p(s, t) \cos\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) ds \right] dx \\
& \qquad + \int_{0}^{1} \Psi_{i}(x) \left[ \sin\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) \int_{1}^{x} \Delta p(s, t) \sin\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) ds \right] dx \\
& \qquad - \left( f_{p}(t) \sin\left(\sum_{j=1}^{M} \Psi_{j}(1) q_{j}(t)\right) \right) \int_{0}^{1} \Psi_{i}(x) \sin\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) dx \\
& \qquad - \left( f_{p}(t) \cos\left(\sum_{j=1}^{M} \Psi_{j}(1) q_{j}(t)\right) \right) \int_{0}^{1} \Psi_{i}(x) \cos\left(\sum_{j=1}^{M} \Psi_{j}(x) q_{j}(t)\right) dx = 0, \quad (4.2)
\end{aligned}
$$

with $\Delta p$ and $f_p$ denoting the dimensionless counterparts of $\Delta P$ and $F_p$, respectively; the prime denotes a spatial derivative.

It is evident from (4.2) that sine and cosine terms have been retained intact in the discretisation scheme and the numerical integration. Consequently, the number of terms in the discretised system and, in turn, the computational cost, increase significantly. Nevertheless, this ensures that reliable predictions are made when the flag deformation becomes very large, at high flow velocities (Tavallaeinejad et al. 2018). In this study, six generalized coordinates are retained in the series expansion of the rotation angle, resulting in 6 second-order ODEs. These equations are then recast into state-space form, i.e. a set of twelve first-order ODEs, and are solved using the pseudo-arclength continuation technique implemented in AUTO (Doedel et al. 1997) to construct
the bifurcation diagrams. This technique enables us to perform continuation of both stable and unstable branches of the solution. Direct time integration is also performed via Gear's backward differentiation formula (BDF) (Gear 1971) yielding the time-varying generalised coordinates $q_i$. Time histories of the amplitude of oscillation are used to plot phase-plane portraits, Fast Fourier transforms (FFT), and (flapping) flag shapes at points of interest in the parameter space. In some cases, where the change in the state of the system is sharp, the stable periodic solution emerging from a saddle-node bifurcation is numerically traced via direct time integration, using the calculated generalised coordinates at a given flow velocity as initial conditions for the next flow velocity. In order to plot the flag shapes in various states, the transverse and longitudinal motions (i.e. $w(x, t)$ and $u(x, t)$, respectively) are retrieved, making use of equation (2.2).

## 4.2. Numerical results

The numerical results from the nonlinear model are shown in the form of bifurcation diagrams, in which the ordinate is the maximum tip rotation of the inverted flag and the abscissa is the dimensionless flow velocity $\Pi$. The system parameters in § 3.2 are again utilized in this section. The convention in all bifurcation diagrams is as follows: continuous and dashed lines represent stable and unstable branches of the solution, respectively. A stable solution represents either a buckled state corresponding to static equilibrium, or the amplitude of oscillation for flapping. Each colour used in the bifurcation diagrams corresponds to a specific regime of dynamical behaviour in the considered range of flow velocities.

Figure 6 gives the bifurcation diagram for a two-dimensional inverted flag with $\mu = 1.0$, demonstrating four distinct states/regimes, which can be summarized as follows.

(i) Stretched-straight state: as also expected from the linear analysis, the trivial solution associated with the undeflected static equilibrium is stable at small flow velocities. Small perturbations generate motions which ultimately die out, and the inverted flag remains stable at its undeflected equilibrium, as illustrated in figure 7(a). This regime was first identified by Kim et al. (2013) and its stability was characterized to some extent by Kim et al. (2013), Gurugubelli & Jaiman (2015) and Sader et al. (2016a).

(ii) Buckled state: the undeflected static equilibrium loses stability at point A where $\Pi = 1.36$ through a supercritical pitchfork bifurcation, in line with the linear model. This leads to divergence (buckling) in the first mode, the amplitude of which increases with flow velocity (see figure 7(b)). The existence of the buckled state was first reported by Gurugubelli & Jaiman (2015), observed experimentally by Yu et al. (2017), and subsequently demonstrated to be an equilibrium of the inverted flag system by Goza et al. (2018).

(iii) Deformed-flapping state: the deflected solution loses stability via a supercritical Hopf bifurcation at point B where $\Pi = 1.70$, leading to “deformed-flapping motions”, that is periodic oscillation around the buckled state. The response is predominantly in the first mode of a cantilevered beam. This motion lasts over the range of $\Pi \in [1.70 \ 1.82]$, and the amplitude of oscillation increases with $\Pi$, as illustrated in figures 7(c)-(d). Additional features of the system dynamics are given in figure 8, displaying the time history, phase-plane portrait and FFT plot at $\Pi = 1.80$, all indicating a periodic motion; the dimensionless frequency of oscillation is $f^* = 0.33$. The appearance of deformed flapping was identified by Yu et al. (2017), and its emergence was reported to be via a supercritical Hopf bifurcation by Goza et al. (2018).

(iv) Large-amplitude flapping state: increasing the flow velocity further causes a transition from the deformed-flapping regime to flapping around the undeflected equilibrium, the amplitude of which increases with the flow velocity. The transition is accompanied by a jump in the amplitude of flapping via two saddle-node bifurcations at points C and D, at $\Pi = 1.82$ and $\Pi = 1.65$, respectively. This forms a region wherein the response of the system is attracted by either the stable limit-cycle around the deflected equilibrium or that around the undeflected equilibrium. Hence, the behaviour of the system is indeed subcritical. The existence of this bistable zone
<figure><img src="image_6.png" /></figure>

Figure 6: Bifurcation diagram for a two-dimensional inverted flag with $\mu = 1.0$, indicating the minimum and maximum values of the leading edge rotation over a cycle of steady-state oscillation of period $T^*$. [—] Undeflected stable equilibrium; [---] deflected stable equilibrium corresponding to buckling of the inverted flag to either side; [—] stable periodic solutions (oscillation around deflected equilibria) corresponding to deformed-flapping of the inverted flag to either side; [---] stable periodic solution (oscillations around undeflected equilibrium). This colour scheme is also utilized in figures 7-10.

<figure><img src="image_7.png" /></figure>

Figure 7: Shapes of the two-dimensional inverted flag studied in figure 6: (a) stretched-straight state ($\Pi = 1.35$), (b) buckled state ($\Pi = 1.40$), (c) at the onset of deformed-flapping mode ($\Pi = 1.71$), (d) deformed-flapping mode ($\Pi = 1.75$), and (e) large-amplitude flapping ($\Pi = 1.65$). Colour scheme as in figure 6.

implies the existence of an unstable repelling limit cycle (not shown in the bifurcation diagram as AUTO was unable to find it) bounded by the two saddle-node bifurcations. Figure 9 shows the time history, phase-plane portrait and FFT plot for $\Pi = 1.65$, all displaying a periodic motion; the dimensionless frequency of oscillation is $f^* = 0.52$. This regime has been characterized in a number of references (see, for instance, Gurugubelli & Jaiman 2015, Shoele & Mittal 2016a and Goza et al. 2018).

The qualitative dynamical behaviour of the inverted flag with increasing flow velocity is illustrated in figure 7, showing successively the stretched-straight, buckled, deformed-flapping, and the large-amplitude flapping states.

Figure 10 shows the sensitivity of $f^*$ to the dimensionless flow velocity $\Pi$. As depicted, $f^*$ sharply decreases with increasing $\Pi$ in the deformed-flapping regime, while it increases only
<figure><img src="image_8.png" /><figcaption>Figure 8: Response of the inverted flag, with $\mu = 1.0$ at $\Pi = 1.80$, displaying a periodic motion around the deflected equilibrium; (a) time history, (b) phase portrait, and (c) the frequency content of the response, obtained from FFT. Colour scheme as in figure 6.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Response of the inverted flag with $\mu = 1.0$ at $\Pi = 1.65$, displaying a periodic motion around the undeflected equilibrium; (a) time history, (b) phase-plane portrait, and (c) the frequency content of the response, obtained from FFT. Colour scheme as in figure 6.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10: Dimensionless frequency $f^*$ versus dimensionless flow velocity $\Pi$ for an inverted flag with $\mu = 1.0$. Colour scheme as in figure 6.</figcaption></figure>
slightly in the large-amplitude flapping regime as $\Pi$ is increased. This latter aspect is in good agreement with the experimental measurements by Sader et al. (2016a), which confirms that the large-amplitude flapping frequency varies only slightly with flow velocity.

The model developed in the present study is believed to simulate fairly accurately the effect of fluid flow on stability of the system, since it is capable of predicting reasonably well the dynamical response of the system up to moderate angles of attack. It is nevertheless realised that the flag exhibits highly curved shapes in the course of large-amplitude flapping (see figure 7), signalling that the limit of applicability of inviscid flow theory may well have been surpassed.

At very large angles of attack, massive separation occurs at the leading edge of the flag, resulting in a sharp drop in the normal force, which in turn causes a sudden oscillation decrease. This phenomenon is not captured by inviscid flow theory which, in principle, cannot model such viscous effects. Interestingly, nevertheless, the ‘qualitative route’ (i.e. stability, buckling, deformed-flapping, and large-amplitude flapping), as well as the bifurcation values for $\Pi$ are in very good agreement with the fully coupled computational study of Goza et al. (2018), who conducted simulations using an immersed boundary method. For example, for a two-dimensional inverted flag with $\mu = 0.2$ (which is the reciprocal of the dimensionless mass ratio used by Goza et al., denoted by $M_\rho = 5.0$), the present model predicts the onset of divergence (pitchfork bifurcation) at $\Pi_A = 1.36$, the onset of deformed flapping (Hopf bifurcation) at $\Pi_B = 1.71$, and the onset of large-amplitude flapping (saddle-node bifurcation) at $\Pi_D = 1.80$, while the computational model of Goza et al. (2018) predicts $\Pi_A = 1.42$, $\Pi_B = 1.55$, and $\Pi_D = 1.78$, respectively. A more complete comparison between the results of the present model and those of Goza et al. (2018) is provided in § 5. It should be noted that the simulations by Goza et al. (2018) were performed in the $Re = 200$ flow regime, while the results obtained via the present model correspond to a higher flow regime, $Re \sim O(10^5)$.

An additional comment is useful here. The buckled and deformed-flapping regimes predicted by the proposed model and the simulations of Goza et al. (2018) have not been observed in the experiments of Kim et al. (2013). Nevertheless, these two regimes are experimentally identified in the measurements of Yu et al. (2017). † Additionally, the existence of the (unstable) deflected equilibrium was observed experimentally by Sader et al. (2016b) via the introduction of additional damping.

The numerical results in figure 6 offer a plausible explanation for the above-mentioned discrepancy. It is known that there is a potential for hysteresis in the onset of flapping, resulting in a bistable region. For a physical flag, the existence of this phenomenon gives rise to the possibility of a sudden jump from quiescent undeflected equilibrium to large-amplitude flapping motion, without encountering other states. Vortex shedding, overshoots during flow-speed increments, or incoming flow turbulence may provide a sufficiently strong disturbance for the jump to occur.

The absence of unsteady vortex shedding in the model proposed in this paper could well explain the discrepancy in the predicted sequence of displayed regimes and experimental observations. However, it should be noted that in the low-$Re$ simulations by Goza et al. (2018) and Ryu et al. (2015) ($Re = 200$), as well as high-$Re$ simulations by Gurugubelli & Jaiman (2019) ($Re = 30000$), vortex shedding is present, yet the deformed-flapping regime is observed in these studies over a range of flow speeds, in contrast to experimental observations of Kim et al. (2013). Indeed, the relationship between experimentally observed regime transitions and the predictions of low- and high-$Re$ simulations is an interesting avenue for further experimental and theoretical research.

The numerical results discussed here suggest that ‘fluid-elastic instabilities’ may be the underlying mechanisms for the rich dynamical behaviour displayed by two-dimensional, strictly

<footnote>† It is noted that Yu et al. (2017) used a channel with a half-width dimension equal to the flag length (and hence, close to the flapping amplitude). Thus, the flow was partially blocked, which may have affected the dynamics.</footnote>
<figure><img src="image_11.png" /></figure>

Figure 11: Bifurcation diagram of an inverted flag subjected to steady-fluid forces showing the leading edge rotation angle versus dimensionless flow velocity. [—] Undeflected stable equilibrium and [---] deflected static equilibrium solutions corresponding to buckling of the inverted flag to either side. No secondary bifurcation (dynamic or static) materializes. Colour scheme as in figure 6.

speaking *heavy* inverted flags. Thus, the flapping motion may be viewed as a 'self-excited oscillation', that is flutter. The same mechanism is responsible for flapping of conventional flags in axial flow. More specifically, the proposed model is developed based on self-excited (or movement-induced) vibration principles and does not include flow instability and its consequences, such as vortex-shedding-induced fluctuating loads. Nevertheless, it does capture qualitatively the bifurcations, as well as some detailed intricate dynamics (including the large-amplitude flapping regime), observed in previous studies for *small* mass ratio flags.

4.3. Influence of mass ratio

As discussed in the foregoing, the undeflected equilibrium of the inverted flag loses stability first via divergence, and the post-buckling response of the system remains static prior to the occurrence of a Hopf bifurcation. As seen from equation (2.24), the mass ratio parameter $\mu$ is present only in the time-dependent terms. Thus, it is expected that $\mu$ has no impact on the critical point for divergence or the existence of the deflected equilibria. This conclusion agrees well with experimental/theoretical findings of Sader et al. (2016a) for the onset of divergence. The existence of a deflected equilibrium following the onset of divergence was realized in their experiments by introducing additional damping by touching the flag with a slender rigid pole.

The mass ratio $\mu$, on the other hand, does affect the dynamical behaviour of the system and significantly affects the stability of the associated solution branches. Setting the time-dependent fluid forces to zero (all terms multiplied by $\mu$ in equation (2.24) disappear) results in governing equations equivalent to a model for an inverted flag subject only to steady-flow forces. As can be seen from figure 11, the deflected static solutions branching off from the pitchfork bifurcation point remain stable with increasing flow velocity. On the other hand, considering all time-dependent fluid forces (with $\mu \neq 0$) and performing a stability analysis leads to the detection of Hopf bifurcations on each deflected static solution, giving rise to periodic solutions, i.e. flutter.

Figure 12 shows how $\mu$ affects the nonlinear response of the system. It is seen that the stability of the trivial solution, the critical flow velocities for divergence and virtually the Hopf bifurcation
<figure><img src="image_12.png" /><figcaption>Figure 12: Bifurcation diagrams for a two-dimensional inverted flag with different mass ratios: (a) $\mu = 0.1$, (b) $\mu = 0.5$, (c) $\mu = 2.0$, and (d) $\mu = 5.0$, indicating the minimum and maximum values of the leading edge rotation over a cycle of steady-state oscillation of period $T^*$. [—] Stable static solution (undeflected stable equilibria); [—] stable static solutions (deflected static equilibria) corresponding to buckling of the inverted flag to either side; [—] stable periodic solutions (oscillation around deflected equilibria) corresponding to deformed-flapping of the inverted flag to either side; [—] stable periodic solution (oscillations around undeflected equilibrium). Colour scheme as in figure 6.</figcaption></figure>

are not affected by $\mu$. The post-critical behaviour, on the other hand, is significantly altered. For instance, for $\mu = 0.5$, the critical values of the flow velocity are $\Pi_A = 1.36$, $\Pi_B = 1.70$, $\Pi_C = 1.83$, and $\Pi_D = 1.75$, while for $\mu = 5.0$ these values are $\Pi_A = 1.36$, $\Pi_B = 1.70$, $\Pi_C = 1.80$ and $\Pi_D = 1.40$.

Two additional points should be made here about the effects of the mass ratio parameter on the dynamics of the system. First, the critical flow velocity for the Hopf bifurcation, $\Pi_B$, obtained via the nonlinear model is considerably lower than that obtained via the linearised model. This may be explained by considering the fact that a linear model assumes that an instability emanates from the position of rest where the deflections are vanishingly small. However, in reality, once a divergence has developed, deflections away from the equilibrium state develop and grow with increasing flow
velocity. As the deflection becomes larger, the nonlinear fluid-related forces become stronger. This
causes a Hopf bifurcation, leading to a limit cycle oscillation, to occur at appreciably different,
and in this case lower flow velocities, compared to the predictions of the linear model.

Secondly, as seen from figure 12, for sufficiently large μ (i.e., for a light flag or a heavily
loaded flag), the transition to large-amplitude flapping becomes strongly subcritical. An unstable
periodic solution (not shown here as AUTO was unable to find it) which originates from the
periodic solution around the deflected equilibrium (i.e. at Π<sub>C</sub>) goes backward and folds at a
flow velocity (i.e. Π<sub>D</sub>) even lower than that for the onset of divergence.† The unstable solution
then becomes stable, generating the large-amplitude periodic solution. The existence of this
solution would give rise to the possibility of a sudden jump from static undeflected equilibrium to
large-amplitude flapping motion, without encountering other states. More specifically, considering
the case of μ = 5.0 (figure 12(d)), two stable solutions coexist for Π ∈ [1.10 1.36]: the trivial
solution and a large-amplitude periodic one. For Π = 1.20, for instance, small perturbations
about the undeflected equilibrium would die out and the inverted flag would come to rest at the
original static equilibrium. Perturbations of sufficiently large amplitude, on the other hand, could
lead the system to enter into the large-amplitude flapping state, without experiencing the usual
route outlined in § 4.2. Interestingly, the results in figure 12 highlight that light flags are more
likely to exhibit hysteresis in the onset of large-amplitude flapping.

Figure 13 shows the variation of the reduced frequency as a function of the dimensionless flow velocity for different values of mass ratio. Here, to be able to compare the numerical values for flapping frequency with those obtained by Goza et al. (2018) and Shoele & Mittal (2016a), a fluid time scale, defined as $\tau_f = L/U$, is used to non-dimensionalize the flapping frequency $f$. Thus, the reduced frequency is defined as $f_R = fL/U = f^* \sqrt{\mu/\Pi}$.‡ As seen in the figure, for all cases studied, $f_R$ decreases monotonically with increasing $\Pi$. For heavier flags (and thus lower mass ratios), the reduced frequency decreases approximately linearly with $\Pi$.

As also seen in figure 13, the reduced frequency decreases as μ is decreased. For example, for
the heaviest flag with μ = 0.1, the reduced frequency is f<sub>R</sub> = 0.09 at Π = 2.0, while for the
lightest flag with μ = 5.0, f<sub>R</sub> = 0.43. The reduction of flapping frequency with mass ratio is in
line with Gurugubelli & Jaiman's predictions. Strikingly, the results in figure 13 indicate that the
dimensionless flapping frequency lies in the range 0.1 ≲ f<sub>R</sub> ≲ 0.2, for μ = 0.3, which agrees
well with predictions by Shoele & Mittal (2016a) and Goza et al. (2018), and measurements
by Kim et al. (2013). This suggests that the absence of unsteady vortex shedding in the model may
have a minimal impact on the dynamical behaviour of heavy flags. For lighter flags (or heavier
fluid forcing), on the other hand, the model overestimates the flapping frequency as compared
to the measurements of Yu et al. (2017) who reported that the dimensionless flapping frequency
reaches to f<sub>R</sub> = 0.18, with the shedding frequency of the vortices being twice the flapping
frequency of the inverted flag. Thus, the unsteady vortex shedding could well be central for light
flags by virtue of injecting another time scale into the problem, and therefore may not be ignored.

The reduced frequency may be considered as a measure of the unsteadiness of the flow.
It is generally acceptable to use quasi-steady-flow theory for a body oscillating with a low
value of reduced frequency in fluid flow (Hodges & Pierce 2011). Small values of the reduced
frequency allow sufficient time for vortex formation and shedding from the trailing edge (Shoele
& Mittal 2016a). Thus, the results shown in figure 13 suggest that the quasi-steady flow theory
may be more appropriate for predicting the dynamics of heavier flags.

<footnote>† Refer to figure 6 for definitions of Π<sub>C</sub> and Π<sub>D</sub>.</footnote>

<footnote>‡ The notation adopted here can be linked to that used in the first footnote in § 2.3.2 as f<sub>R</sub> = k/π.</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 13: Sensitivity of the large-amplitude flapping reduced frequency to the dimensionless flow velocity Π for inverted flags of different mass ratios: [▲] μ = 0.1; [■] μ = 0.5; [●] μ = 2.0; and [♦] μ = 5.0.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14: Inverted flag in axial flow with a small initial angle of attack ψ₀</figcaption></figure>

### 4.4. Influence of angle of attack

It is of interest to explore the dynamical behaviour of inverted flags when they are subject to oncoming flow with a non-zero angle of attack ψ₀, as illustrated in figure 14. Note that ψ₀ is constant with respect to temporal and spatial coordinates.

The projection of the undisturbed flow velocity onto the X and Z directions leads to $U_X = U \cdot e_X = U \cos \psi_0$ and $U_Z = U \cdot e_Z = U \sin \psi_0$, respectively. Consequently, the normal and tangential components of the relative flow velocity in equation (2.3), and in turn the fluid-related forces, should be revised in order to take into account both components. Hence,

$$V_n(x, t) = (\dot{w} - U \sin \psi_0) \cos \psi - (\dot{u} + U \cos \psi_0) \sin \psi, \quad (4.3a)$$

$$V_{\tau}(x, t) = (\dot{u} + U \cos \psi_0) \cos \psi + (\dot{w} - U \sin \psi_0) \sin \psi. \quad (4.3b)$$

Inserting (4.3a) into (2.12) and substituting the resulting expression together with (4.3b) into (2.21) and (2.19) gives the total normal force acting on the flag. On the structural dynamics side, the
<figure><img src="image_15.png" /><figcaption>Figure 15: Nonlinear response of a two-dimensional inverted flag with $\mu = 1.0$ and $\psi_0 = 0.1$. (a) Bifurcation diagram indicating the minimum and maximum values of the leading edge rotation over a cycle of steady-state oscillation of period $T^*$: [—] Stable static solution; [---] unstable and [—] stable periodic solution (asymmetric oscillations). (b) Dimensionless frequency $f^*$ versus dimensionless flow velocity $\Pi$.</figcaption></figure>

equations remain unchanged. The resulting equation of motion is solved using the techniques discussed in § 4.1.

The bifurcation diagram and the flapping frequency for a 2D inverted flag with $\mu = 1.0$, and $\psi_0 = 0.1 \approx 5.7^\circ$ are plotted in figure 15. Even such a seemingly small angle of attack gives rise to asymmetric fluid loading on the inverted flag which consequently breaks the symmetry of the nonlinear response. More specifically, the supercritical pitchfork bifurcation ‘point’ (point A in figure 6) is replaced by a gradual transition from small to large deflections. In contrast to the case of $\psi_0 = 0$, the flag bends continuously more with increasing flow velocity. This behaviour is followed by a Païdoussis-type Hopf bifurcation, giving rise to a stable limit cycle corresponding to asymmetric flapping around the deflected equilibrium. This solution folds at the first saddle-node bifurcation encountered (at $\Pi = 1.8$) and becomes unstable. By tracing the unstable solution, a second saddle-node bifurcation arises at $\Pi = 1.63$, where the response of the system folds once more. Following the second saddle-node bifurcation, the solution becomes stable, which is physically manifested as asymmetric large-amplitude flapping around the origin.

As seen in figure 15(b), except for a sharp increase of $f^*$ within a short range in the bistable regime, the behaviour is generally similar to that for the higher mass ratio $\mu = 1.0$ and $\psi_0 = 0$ (see figure 10).

The bifurcation diagrams shown in figure 16 for $\psi_0 = 0.2$ illustrate the effects of the resulting symmetry breaking on the response of the system for various mass ratios (cf. figure 12 for $\psi_0 = 0$). Interestingly, the broken symmetry as a result of a non-zero initial angle of attack enables us to obtain numerically the unstable solution branch connecting the two periodic solution branches – a task which was not feasible for $\psi_0 = 0$; see § 4.2 and § 4.3 and figures 6 and 12.

As discussed in § 4.2, the periodic solution around the deflected equilibrium becomes unstable through a saddle-node bifurcation and recovers stability via another saddle-node bifurcation at a lower flow velocity. Different scenarios can be expected for flags of different mass ratios. While the sequence of progressive bending → asymmetric deformed-flapping → asymmetric large-amplitude flapping is more plausible for the inverted flag with a low mass ratio, e.g. $\mu = 0.1$, the
<figure><img src="image_16.png" /><figcaption>Figure 16: Bifurcation diagrams for $\psi_0 = 0.2$ and various mass ratios: (a) $\mu = 0.1$, (b) $\mu = 0.5$, (c) $\mu = 2.0$, and (d) $\mu = 5.0$, indicating the minimum and maximum values of the leading edge rotation over a cycle of steady-state oscillation of period $T^*$. [---] Stable static solution; [---] unstable and [---] stable periodic solutions (asymmetric oscillations).</figcaption></figure>

flag with the higher mass ratio $\mu = 5.0$ may go directly from the buckled shape to large-amplitude flapping due to the appearance of the limit cycle in the proximity of the static solution. This behaviour was also observed for flags with zero initial angle of attack, where an abrupt jump may occur from a static equilibrium to large-amplitude flapping (see § 4.3). However, a non-zero $\psi_0$ seems to facilitate the direct transition of light flags from the buckled shape to large-amplitude flapping.

The effect of the mass ratio on the flapping frequency of inverted flags with an initial incidence angle $\psi_0 = 0.2$ is shown in figure 17. Here, the reduced frequency for both small- and large-amplitude flapping is shown. The branch in the low-$\Pi$ region corresponds to small-amplitude flapping, while that appearing in the high-$\Pi$ region corresponds to large-amplitude flapping. It is seen that the reduced frequency generally decreases as the flow velocity is increased, for both small- and large-amplitude flapping. However, in all the cases shown, a jump occurs in the value of the reduced frequency, in the region where the asymmetric small-amplitude flapping is
<figure><img src="image_17.png" /><figcaption>Figure 17: Reduced frequency of the asymmetric flapping motion as a function of the dimensionless flow velocity $\Pi$ for $\psi_0 = 0.2$ and different mass ratios: [▲] $\mu = 0.1$, [■] $\mu = 0.5$, [●] $\mu = 2.0$, and [♦] $\mu = 5.0$. [---] Unstable and [—] stable periodic solutions.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 18: Bifurcation diagram with $\mu = 1.0$ and various initial angles of attack $\psi_0$.</figcaption></figure>

transformed into large-amplitude flapping. Similarly to what was found for the flag with zero $\psi_0$ (see figure 13), the reduced frequency decreases with decreasing mass ratio.

The sensitivity of the nonlinear response of the system to the initial angle of attack $\psi_0$ is shown in figure 18. The present numerical results are in fairly good agreement with the observations provided by Cossé et al. (2014) who investigated experimentally the effect of flag orientation to the impinging flow. In particular, they reported that the transition for cases with $\psi_0 \neq 0$ was gradual, as opposed to the case with $\psi_0 = 0$, where an abrupt jump to large-amplitude flapping was observed. Likewise, the numerical results show that the transition between the two flapping motions becomes smoother as $\psi_0$ is increased. Furthermore, varying $\psi_0$ affects the onset of flapping significantly. For inverted flags with larger initial angles of attack, the asymmetric flapping motion around the buckled shape occurs at lower flow velocities. The
amplitude of large-amplitude flapping also decreases correspondingly. This also agrees well with measurements by Cossé et al. (2014) which reveal that, as the initial angle of attack is increased, the transition from static equilibrium to the flapping mode occurs at lower flow velocities and with smaller amplitudes of oscillation. Moreover, as seen in figure 18, the inverted flag undergoes the asymmetric deformed-flapping motion at different critical flow velocities, which is consistent with computational findings by Shoele & Mittal (2016b) and, again, the experimental observations of Cossé et al. (2014).

These findings highlight the ability of inverted flags to perform a flapping motion even when the oncoming flow is not parallel to their longitudinal axis. This is of practical importance for designing robust energy harvesters in actual ambient conditions where the wind direction is not necessarily parallel to the neutral plane of the inverted flag (Orrego et al. 2017).

**5. Comparison with simulations of Goza et al. (2018)**

Goza et al. (2018) examined the dynamical behaviour of two-dimensional inverted flags for four different values of the mass ratio, defined as $M_ρ = μ^{-1}$, at $Re = 20$ and $Re = 200$. They performed computational investigations by solving strongly coupled 2D Navier-Stokes equations and the geometrically nonlinear Euler-Bernoulli beam equation. The simulation results were summarised in the form of bifurcation diagrams in which the peak values of tip displacement of the flag were plotted as a function of the dimensionless bending stiffness of the flag, denoted as $K_B = Π^{-2}$. Here, the numerical results from the present model for the dynamics of an inverted flag with $μ = 0.2$ are compared with those obtained by Goza et al. (2018) for $M_ρ = 5.0$ and at the above two Reynolds numbers.

Figure 19 shows the peak values† of the tip deflection of the 2D inverted flag over a cycle of steady-state oscillation of period $T^*$ as a function of the dimensionless flow velocity. The solid lines show our model predictions, while symbols show the results from Goza et al. (2018). A good agreement on the global dynamical behaviour between the two sets of results may be observed. In particular, the type of instabilities and their sequence are identical across the two sets of results. Interestingly, also the amplitudes for different dynamical regimes predicted by the present model are comparable to those computed by Goza et al. (2018) at $Re = 20$.

A quantitative comparison is made in Table 2 where thresholds of divergence, deformed-flapping motion, and large-amplitude flapping motion predicted by the present model and the computational model of Goza et al. (2018) are given. As seen from the table, excellent agreement was achieved. The maximum relative error is about 8%.

Moreover, large-amplitude flapping frequencies obtained by the two models are in the same ballpark. For example, for $μ = 0.2$ at $Π = 2.0$, $f_R = 0.13$ is obtained by our model, while $f_R \simeq 0.09$ is obtained by Goza et al. (2018) at $Re = 200$. The discrepancy may be explained by the fact that viscous flow effects have not systematically been considered in our model. These effects perhaps can produce a distributed compressive load along the plate which would tend to reduce the effective stiffness of the fluid-structure system and thus to decrease the flapping frequency.

An additional point ought to be made here concerning the Reynolds-number independence in the dynamics of inverted flags. Interestingly, the qualitative comparison between the results obtained via the proposed model and the simulations of Goza et al. (2018) at low Reynolds numbers (i.e. $Re = 200$ and $Re = 20$) suggests that the overall dynamical behaviour and characteristic features of inverted flags is insensitive to the Reynolds number. This has been shown by Shoele & Mittal (2016a) who performed direct numerical simulations at small Reynolds

<footnote>† Two sets of values for each $Π$ can be seen for a range of flow velocities on either positive or negative side of the bifurcation diagram, as indicated by [♦]. This is perhaps because the local maxima and minima that appeared at the peak values of displacement (where tip transverse velocity goes thorough zero) have been included in the bifurcation diagram (Goza et al. 2018, figure 14(a)).</footnote>
<figure><img src="image_19.png" /><figcaption>Figure 19: Bifurcation diagrams for inverted flags with $\mu = 0.2$ over a cycle of period $T^*$ steady-state oscillation. [—] stretched-straight state; [—] buckled state; [—] deformed-flapping regime; and [—] large-amplitude flapping around the origin. Lines show the results obtained by the present model, while symbols [♦], [♦], and [♦] show the corresponding results obtained by Goza et al. (2018). (a) $Re = 200$, and (b) $Re = 20$.</figcaption></figure>

<table><caption>Table 2: Dimensionless flow velocity for various bifurcation points: comparison between the present model and simulations by Goza et al. (2018);† simulations for $Re = 200$, ‡ simulations for $Re = 20$. Values extracted from (Goza et al. 2018, figure 3).</caption><thead><tr><th rowspan="2">Studies</th><th>Buckling</th><th>Deformed-flapping</th><th>Large-amplitude flapping</th></tr><tr><th>(Pitchfork bifurcation)<br>$\Pi_A$</th><th>(Hopf bifurcation)<br>$\Pi_B$</th><th>(Saddle-node bifurcation)<br>$\Pi_C$</th></tr></thead><tbody><tr><td>Present study</td><td>1.36</td><td>1.71</td><td>1.80</td></tr><tr><td rowspan="2">Goza et al. (2018)</td><td>1.42†</td><td>1.55†</td><td>1.78†</td></tr><tr><td>1.48‡</td><td>1.61‡</td><td>1.85‡</td></tr></tbody></table>

numbers and reported that the frequency and the amplitude of large-amplitude flapping remain almost unchanged in the range $100 \le Re < 1000$. Tang et al. (2015) explored the Reynolds number effect in the range $100 \le Re \le 500$ and found that the flag inertia dominates over the viscous effects with increasing $Re$, leading to fairly similar results for $Re > 100$. Furthermore, a comparison between the high Reynolds number simulations of Gurugubelli & Jaiman (2019) at $Re = 30000$ and those by others at small Reynolds numbers (i.e., Tang et al. 2015, Ryu et al. 2015, Gurugubelli & Jaiman 2015, Gilmanov et al. 2015, Shoele & Mittal 2016a) shows that the key features of the dynamics of the inverted flag do not change significantly with $Re$.

## 6. Conclusions

The main concern of the present paper is the mechanism for flapping of infinitely high aspect ratio heavy inverted flags in axial flow. It is shown that it may well be a ‘fluidelastic instability’. This means that flapping of two-dimensional inverted flags may be viewed as a self-excited
(or movement-induced) vibration, similar to flutter of aircraft wings or of cantilevered pipes conveying fluid.

Our numerical results show that two-dimensional inverted flags undergo multiple bifurcations as the flow velocity is increased. Physically speaking, the flags exhibit four dynamical states (or regimes) as the flow velocity is increased: (i) stretched-straight state, (ii) buckled state, (iii) deflected-flapping state, and (iv) large-amplitude flapping state. Our findings suggest that the system dynamics depends on the mass ratio as well as on the initial angle of attack. It was shown that the mass ratio parameter does not affect the stability of the stretched-straight state and the onset of divergence, as they are static phenomena; however, it controls the possibility of direct transition from undeflected equilibrium to large-amplitude flapping motion and it affects the amplitude of large-amplitude flapping. Furthermore, inverted flags of larger mass ratios are more likely to exhibit hysteresis, and thus more prone to undergo flapping motion at lower flow velocities due to the presence of a subcritical periodic solution. At a certain critical flow velocity (via a saddle-node bifurcation), flow-induced disturbances may result in spontaneous large-amplitude flapping of inverted flags of sufficiently large mass ratio.

In addition, it was shown that the symmetry of the system response breaks down when the initial angle of attack is not zero. As this angle is increased, the system loses stability at lower flow velocities, while the amplitude of large-amplitude flapping decreases. In addition, smoother transitions ensue between the dynamical regimes. These results were found to be in qualitative good agreement with the experimental observations of Cossé et al. (2014). Where possible, the predictions based on the present model were also compared with experimental observations of Sader et al. (2016a), as well as computational results of Goza et al. (2018). In general, very good agreement was achieved. In particular, the present analytical model predicts a similar sequence of instabilities to that obtained computationally by Goza et al.. The values for critical flow velocities, corresponding to the initiation of each regime, were also comparable. Interestingly, the qualitative comparison made in § 5 suggests that the global dynamics of inverted flags is independent of the Reynolds number. However, the comparison for some aspects of the post-critical dynamics, such as amplitudes of buckling/flapping, is less satisfactory. The discrepancy can probably be explained by the fact that the analytical model does not accurately account for the effects of massive flow separation from the leading edge of the inverted flag. These effects were considered to some extent by employing the leading-edge suction analogy proposed by Polhamus (1966) for the calculation of the vortex or nonlinear lift, but clearly not sufficiently well.

As evidenced from the review of previous work in the Introduction, many studies on the dynamics of inverted flags proposed that vortex-induced vibration (VIV) is the underlying mechanism for flapping (Sader et al. 2016a, for example).† Nevertheless, the model proposed in the present paper, though totally ignoring vortex shedding excitation, does predict many of the salient features of the dynamics of heavy inverted flags. Figure 19 bears witness to that, convincingly. Furthermore, flapping has been found to arise in very low Reynolds number flows (Re = 20) by Goza et al. (2018), where classical vortex shedding was not observed. Therefore, although not casting doubts that VIV is important, or even determining, in the generation of flapping in the case of light flags, the work in the present paper shows that self-excitation is the principal underlying mechanism for heavy flags; even so, we recognize that further work should be done to establish this unequivocally.

The findings of the present paper are based on several assumptions, such as that the flow over the flag is a two-dimensional quasi-steady flow. These assumptions might limit the capability of the model to accurately predict the dynamics of the flag, for example, when very large flapping

<footnote>† There is also the possibility of a combined vortex-shedding (VIV) and self-excited (SEV) fluidelastic excitation mechanism, or VIV-SEV, in a similar way as found for slender prismatic bodies with bluff cross-section and sufficiently long afterbody by Nemes et al. (2012) and Mannini et al. (2014).</footnote>
motion occurs. Nevertheless, the findings of this paper will hopefully motivate further research on the analytical modelling of inverted flags subject to axial flow. Such modelling can provide deep insights into, for example, the underlying mechanisms for the various instabilities, and the impact of different parameters on the global as well as local dynamics and stability of the system. They can also be used to corroborate experimental and computational studies.

Finally, it ought to be mentioned also that the findings of this study should be of practical importance for the design of small-scale energy harvesting systems.

## Acknowledgments

The authors gratefully acknowledge the financial support by the Natural Sciences and Engineering Research Council of Canada, the Solution Mining Research Institute (SMRI) and Pipeline Research Council International (PRCI). The last author also appreciates the financial support by Concordia University through a Start-Up grant.

## Appendix A. Analytical solution to the aeroelastic response

The pressure distribution along the inverted flag depends on the normal component of relative velocity via the singular integral formulated in equation (2.12), which in turn is a function of flag deformation. This integral may be evaluated using Glauert's principal value integral given as

$$ \frac{1}{\pi} \int_{0}^{\pi} \frac{\cos n\varphi}{\cos \varphi - \cos \vartheta} d\varphi = \frac{\sin n\vartheta}{\sin \vartheta}, \quad n = 0, 1, 2, \dots \quad (A1) $$

To enable the use of (A1), polynomial representation of exact solutions of a linear uniform Euler-Bernoulli cantilevered beam is utilized, using Chebyshev polynomial expansions of the first kind, which are used in the discretization procedure. It is noted that convergence and accuracy of the solution obtained via the Galerkin method depend on the number and type of trial functions in the expansion series. Hence, a sufficiently large number of terms in the Chebyshev polynomial expansions are considered to guarantee that each polynomial-type comparison function mimics its associated exact counterpart, even up to the fourth spatial derivative. The polynomial eigenfunctions are written as

$$ W_i(x) = \sum_{j=1}^{J_i} a_{ij} x^j, \quad (A2) $$

with $a_{ij}$ being the constants of the $j$th term associated with the $i$th eigenfunction to be approximated; $J_i$ represents the maximum power needed for approximation of the $i$th eigenfunction. Recalling that $x = \frac{1}{2}(1 + \cos \varphi)$, and making use of the binomial theorem as well as the multiple angle formula of the sine and cosine functions, the transverse displacement of the inverted flag may be formulated as

$$ w(\varphi, t) = \sum_{i=1}^{M} p_i(t) \sum_{j=0}^{J_i} a_{i,j} \frac{1}{2^j} \sum_{m=0}^{j} \binom{j}{m} \cos^m \varphi. \quad (A3) $$

The perturbation pressure may be obtained analytically upon inserting (A3) into (3.1), and applying Glauert's integral equality (A1) to perform integrations over $\varphi$. The resultant expression
may be decomposed into three terms, associated with three effects of the pressure force in (3.1):

$$
\begin{align}
\Delta P(\theta, t) &= \mu \sum_{i=1}^{M} \left( \ddot{p}_i(t) \int_{\theta}^{0} (\Sigma 1_i - \Sigma 2_i) \, d\vartheta \right) + \Pi^2 \sum_{i=1}^{M} p_i(t) (\Sigma 3_i - \Sigma 4_i - \Sigma 6_i) \nonumber \\
&\quad - \Pi \sqrt{\mu} \sum_{i=1}^{M} \left[ \dot{p}_i(t) \left( \int_{\theta}^{0} (\Sigma 3_i - \Sigma 4_i) \, d\vartheta + (\Sigma 1_i - \Sigma 2_i - \Sigma 5_i) \right) \right] \tag{A 4} \\
&\equiv \mu \sum_{i=1}^{M} \ddot{p}_i(t) P^{\mathcal{M}}(\theta) + \Pi^2 \sum_{i=1}^{M} p_i(t) P^{\mathcal{K}}(\theta) - \Pi \sqrt{\mu} \sum_{i=1}^{M} \dot{p}_i(t) P^{\mathcal{C}}(\theta), \nonumber
\end{align}
$$

with $P^M(\theta)$, $P^K(\theta)$, and $P^K(\theta)$ being the components of $\Delta P(\theta, t)$ obtained by retaining terms associated with $\ddot{p}_i(t)$, $\dot{p}_i(t)$, and $p_i(t)$ respectively. The values of $\Sigma 1_i - \Sigma 6_i$ are defined as

$$
\begin{align}
\Sigma 1_i &= \sum_{j=1}^{J_i} a_{i,j} \frac{1}{2^j} \sum_{m=0}^{j} \binom{j}{m} 2^{1-m} \sum_{n=0}^{\frac{m-1}{2}} \binom{m}{n} \Sigma 1_n, \notag \\
\Sigma 2_i &= \sum_{j=1}^{J_i} a_{i,j} \frac{1}{2^j} \sum_{m=0}^{j} \binom{j}{m} \frac{1}{2^{m+1}} \sum_{n=0}^{\frac{m-1}{2}} \binom{m+2}{n} \Sigma 2_n, \notag \\
\Sigma 3_i &= \sum_{j=1}^{J_i} a_{i,j} \frac{j}{2^{j-1}} \sum_{m=0}^{j} \binom{j-1}{m} 2^{1-m} \sum_{n=0}^{\frac{m-1}{2}} \binom{m}{n} \Sigma 1_n, \tag{A 5} \\
\Sigma 4_i &= \sum_{j=1}^{J_i} a_{i,j} \frac{j}{2^{j-1}} \sum_{m=0}^{j} \binom{j-1}{m} \frac{1}{2^{m+1}} \sum_{n=0}^{\frac{m-1}{2}} \binom{m+2}{n} \Sigma 1_n, \notag \\
\Sigma 5_i &= \sum_{j=0}^{J_i} a_{i,j} \frac{1}{2^j} \int_0^\pi (\cos \varphi + 1)^j (1 - \cos \varphi) \, d\varphi, \notag
\end{align}
$$

$$
\Sigma 6_i = \sum_{j=0}^{J_i} a_{i,j} \frac{j}{2^{j-1}} \int_{0}^{\pi} (\cos \varphi + 1)^{j-1} (1 - \cos \varphi) \, d\varphi,
$$

where

$$
\begin{align}
\Sigma 1_n &= \sum_{r=0}^{m-2n} \binom{m-2n}{r} \cos^r \vartheta \sin^{k-1} \vartheta \sin\left(k \frac{\pi}{2}\right), \notag \\
\Sigma 2_n &= \sum_{r=0}^{m-2n+2} \binom{m-2n+2}{r} \cos^r \vartheta \sin^{k+1} \vartheta \sin\left((k+2) \frac{\pi}{2}\right), \tag{A 6}
\end{align}
$$

and $k = m - 2n - r$.

Hence, integrations over $\varphi$ are performed exactly and the pressure distribution over the inverted
<figure><img src="image_20.png" /><figcaption>Figure A.1: First six in vacuo comparison functions of a cantilevered beam utilized in the Galerkin discretization of the nonlinear system. First [—], second [—], third [—], fourth [—], fifth [—] and sixth [—] comparison function.</figcaption></figure>

flag may be evaluated analytically. Applying the Galerkin technique to (3.1) yields

$$ 
\begin{align*}
M_{ij} &= \int_{0}^{\pi} \frac{1}{2} (\mu P_{i}^{M}(\theta) + W_{i}(\theta)) W_{j}(\theta) \sin \theta \, d\theta, \\
C_{ij} &= -\Pi \sqrt{\mu} \int_{0}^{\pi} \frac{1}{2} P_{i}^{C}(\theta) W_{j}(\theta) \sin \theta \, d\theta + \eta_{s} \int_{0}^{1} \partial_{xxxx} W_{i}(x) W_{j}(x) \, dx, \tag{A 7} \\
K_{ij} &= \mu \int_{0}^{\pi} \frac{1}{2} P_{i}^{K}(\theta) W_{j}(\theta) \sin \theta \, d\theta + \int_{0}^{1} \partial_{xxxx} W_{i}(x) W_{j}(x) \, dx.
\end{align*}
$$

For the nonlinear analysis, the following comparison functions are defined by their polynomial successors to discretize the nonlinear equation of motion (2.24):

$$ \Psi_i(x) = \sin(\lambda_i x) + \sinh(\lambda_i x) + \frac{\cos(\lambda_i) + \cosh(\lambda_i)}{\sin(\lambda_i) + \sinh(\lambda_i)} (\cos(\lambda_i x) - \cosh(\lambda_i x)), \quad (A 8) $$

recalling that $\lambda_i$ is the ith root of the transcendental equation $1 + \cosh(x) \cos(x) = 0$ (Bishop & Johnson 2011). These functions are suitable since they satisfy the clamped-free boundary conditions, and also convenient, as they are already known (note that $\psi(x, t) = \partial_x w(x, t)$ holds for the linear equation of motion of the flag; hence $\Psi_i(x) = \partial_x W_i(x)$ are chosen). Similarly to equation (A 2), polynomial representations of these comparison functions are constructed with the aid of a Chebyshev polynomial expansion of the first kind with a sufficiently large number of terms (i.e. $\Psi_i(x) = \sum_{j=1}^{J_i} a_{ij} j x^{j-1}$). The first six comparison functions are shown in figure A.1.

It is clear from (2.3) that the normal component of the relative velocity is a function of trigonometric expressions in terms of $\psi(x, t)$. For the problem at hand, retaining these trigonometric terms intact hinders the evaluation of the exact solution of the non-circulatory velocity potential, by virtue of the fact that the integral is singular and hence its principal value needs to be calculated. In order to overcome this difficulty, the normal component of the relative velocity is expressed in a succession of Taylor's expansions. This results in a polynomial representation of the source/sink
sheet strength in terms of the rotation angle. Hence,

$$ V_n(x,t) \approx -\Pi \sqrt{\mu} \left[ \partial_t \left( \int_0^x \sum_{r=1}^R \frac{(-1)^r}{(2r)!} \psi(s,t)^{2r} \, ds \right) + 1 \right] \sum_{r=1}^R \frac{(-1)^r}{(2r+1)!} \psi(x,t)^{2r+1} \\ + \Pi \sqrt{\mu} \sum_{r=1}^R \frac{(-1)^r}{(2r)!} \psi(x,t)^{2r} \partial_t \left( \int_0^x \sum_{r=1}^R \frac{(-1)^r}{(2r+1)!} \psi(s,t)^{2r+1} \, ds \right), \quad (A9) $$

where R is the order of truncation. Similarly to the procedure detailed earlier, inserting $\Psi_i(x,t) = \sum_{i=1}^M q_i(t) (\sum_{j=1}^{J_i} a_{ij} j x^{j-1})$ into (A9), the velocity potential given in (2.13) can be integrated based on Glauert's integral. This enables the analytical evaluation of the pressure distribution over the flag, in a similar fashion to the discretization technique detailed for the linear model. In this paper, terms up to third order are retained in the Taylor expansions in equation (A9), resulting in

$$ V_n(x,t) \approx \Pi \sqrt{\mu} \left[ \psi(x,t) \int_0^x \dot{\psi}(s,t) \psi(s,t) \, ds + \left(1 - \frac{1}{2}\psi(x,t)^2\right) \int_0^x \dot{\psi}(s,t) \, ds \right. \\ \left. - \frac{1}{2} \int_0^x \dot{\psi}(s,t) \psi(s,t)^2 \, ds \right] - \psi(x,t) + \frac{1}{6} \psi(x,t)^3. \quad (A10) $$

## REFERENCES

* ABDELSBAKI, A., PAÏDOUSSIS, M. & MISRA, A. 2018 A nonlinear model for a free-clamped cylinder subjected to confined axial flow. *Journal of Fluids and Structures* **80**, 390–404.

* ANDERSON, J. D. J. 2017 *Fundamentals of Aerodynamics*, sixth edn. New York, NY: McGraw-Hill Education.

* AZUMA, A. & OKAMOTO, M. 2005 Theoretical study on two-dimensional aerodynamic characteristics of unsteady wings. *Journal of Theoretical Biology* **234** (1), 67–78.

* BISHOP, R. E. D. & JOHNSON, D. C. 2011 *The Mechanics of Vibration*. Cambridge University Press.

* BISPLINGHOFF, R. L., ASHLEY, H. & HALFMAN, R. L. 2013 *Aeroelasticity*. Courier Corporation.

* BLEVINS, R. D. 1990 *Flow-Induced Vibration*, 2nd edn. New York, NY (USA); Van Nostrand Reinhold Co., Inc.

* BOLLAY, W. 1939 A non-linear wing theory and its application to rectangular wings of small aspect ratio. *ZAMM-Journal of Applied Mathematics and Mechanics* **19** (1), 21–35.

* CONNELL, B. S. & YUE, D. K. 2007 Flapping dynamics of a flag in a uniform stream. *Journal of Fluid Mechanics* **581**, 33–67.

* COSSÉ, J., SADER, J., KIM, D., CERDEIRA, C. H. & GHARIB, M. 2014 The effect of aspect ratio and angle of attack on the transition regions of the inverted flag instability. In *ASME 2014 Pressure Vessels and Piping Conference*. American Society of Mechanical Engineers.

* DEVORIA, A. C. & MOHSENI, K. 2017 A vortex model for forces and moments on low-aspect-ratio wings in side-slip with experimental validation. *Proceedings of the Royal Society (London)* A **473** (2198), 20160760.

* DOEDEL, E., CHAMPNEYS, A., FAIRGRIEVE, T., KUZNETSOV, Y., SANDSTEDE, B. & WANG, X. 1997 *Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont)*. AUTO97, Concordia University, Montreal, QC, Canada.

* DONE, G. T. & SIMPSON, A. 1977 Dynamic instability of certain conservative and non-conservative systems. *Journal of Mechanical Engineering Science* **19** (6), 251–263.

* DOWELL, E. H. 1975 *Aeroelasticity of Plates and Shells*. Neatherlands: Noordhoff International Publishing.

* ELOY, C., LAGRANGE, R., SOUILLIEZ, C. & SCHOUVEILER, L. 2008 Aeroelastic instability of cantilevered flexible plates in uniform flow. *Journal of Fluid Mechanics* **611**, 97–106.

* FAN, B., HUERTAS-CERDEIRA, C., COSSÉ, J., SADER, J. E. & GHARIB, M. 2019 Effect of morphology on the large-amplitude flapping dynamics of an inverted flag in a uniform flow. *Journal of Fluid Mechanics* **874**, 526–547.

* FAROKHI, H., PAÏDOUSSIS, M. P. & MISRA, A. K. 2018 Nonlinear behaviour of cantilevered carbon
nanotube resonators based on a new nonlinear electrostatic load model. *Journal of Sound and Vibration* **419**, 604–629.

GEAR, C. W. 1971 The automatic integration of ordinary differential equations. *Communications of the ACM* **14** (3), 176–179.

GILMANOV, A. , LE, T. B. & SOTIROPOULOS, F. 2015 A numerical approach for simulating fluid structure interaction of flexible thin shells undergoing arbitrarily large deformations in complex domains. *Journal of Computational Physics* **300**, 814–843.

GOSSELIN, F. , DE LANGRE, E. & MACHADO-ALMEIDA, B. A. 2010 Drag reduction of flexible plates by reconfiguration. *Journal of Fluid Mechanics* **650**, 319–341.

GOZA, A. , COLONIUS, T. & SADER, J. E. 2018 Global modes and nonlinear analysis of inverted-flag flapping. *Journal of Fluid Mechanics* **857**, 312–344.

GUO, C. Q. & PAÏDOUSSIS, M. P. 2000 Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow. *Journal of Applied Mechanics* **67** (1), 171–176.

GURUGUBELLI, P. & JAIMAN, R. 2019 Large amplitude flapping of an inverted elastic foil in uniform flow with spanwise periodicity. *Journal of Fluids and Structures* **9**, 139–163.

GURUGUBELLI, P. S. & JAIMAN, R. K. 2015 Self-induced flapping dynamics of a flexible inverted foil in a uniform flow. *Journal of Fluid Mechanics* **781**, 657–694.

HODGES, D. H. & PIERCE, G. A. 2011 *Introduction to Structural Dynamics and Aeroelasticity*. Cambridge University Press.

HOERNER, S. F. 1965 *Fluid-Dynamic Drag*. Bricktown, NJ: Published by the author.

HOERNER, S. F. & BORST, H. V. 1985 *Fluid-Dynamic Lift: Practical Information on Aerodynamic and Hydrodynamic Lift*, 2nd edn. Bricktown, NJ: Mrs. Liselotte A. Hoerner.

HOWELL, R. , LUCEY, A. , CARPENTER, P. & PITMAN, M. W. 2009 Interaction between a cantilevered-free flexible plate and ideal flow. *Journal of Fluids and Structures* **25** (3), 544–566.

HUANG, L. 1995a *Flutter of cantilevered plates in axial flow*. *Journal of Fluids and Structures* **9** (2), 127–147.

HUANG, L. 1995b *Mechanical modeling of palatal snoring*. *The Journal of the Acoustical Society of America* **97** (6), 3642–3648.

KIM, D. , COSSÉ, J. , HUERTAS CERDEIRA, C. & GHARIB, M. 2013 Flapping dynamics of an inverted flag. *Journal of Fluid Mechanics* **736**.

KORNECKI, A. , DOWELL, E. H. & O'BRIEN, J. E. 1976 On the aeroelastic instability of two-dimensional panels in uniform incompressible flow. *Journal of Sound and Vibration* **47** (2), 163–178.

LIGHTHILL, M. J. 1971 Large-amplitude elongated-body theory of fish locomotion. *Proceedings of the Royal Society of London B* **179** (1055), 125–138.

LOPES, J.-L. , PAÏDOUSSIS, M. P. & SEMLER, C. 2002 Linear and nonlinear dynamics of cantilevered cylinders in axial flow. Part 2: The equations of motion. *Journal of Fluids and Structures* **16** (6), 715–737.

MANNINI, C. , MARRA, A. & BARTOLI, G. 2014 VIV-galloping instability of rectangular cylinders: Review and new experiments. *Journal of Wind Engineering and Industrial Aerodynamics* **132**, 109–124.

MEI, C. , ABDEL-MOTAGALY, K. & CHEN, R. 1999 Review of nonlinear panel flutter at supersonic and hypersonic speeds. *Applied Mechanics Reviews* **52** (10), 321–332.

NAYFEH, A. H. & PAI, P. F. 2008 *Linear and Nonlinear Structural Mechanics*. John Wiley & Sons.

NEMES, A. , ZHAO, J. , JACONO, D. L. & SHERIDAN, J. 2012 The interaction between flow-induced vibration mechanisms of a square cylinder with varying angles of attack. *Journal of Fluid Mechanics* **710**, 102–130.

ORREGO, S. , SHOELE, K. , RUAS, A. , DORAN, K. , CAGGIANO, B. , MITTAL, R. & KANG, S. H. 2017 Harvesting ambient wind energy with an inverted piezoelectric flag. *Applied Energy* **194**, 212–222.

PAÏDOUSSIS, M. P. 2014 *Fluid-Structure Interactions: Slender Structures and Axial Flow. Volume 1*, 2nd edition. Oxford: Academic Press.

PAÏDOUSSIS, M. P. 2016 *Fluid-Structure Interactions: Slender Structures and Axial Flow. Volume 2*, 2nd edition. Oxford: Academic Press.

PAÏDOUSSIS, M. P. , PRICE, S. J. & DE LANGRE, E. 2010 *Fluid-Structure Interactions: Cross-Flow-Induced Instabilities*. Cambridge University Press.

PEDERSEN, C. B. & ŽBIKOWSKI, R. 2006 An individual-Polhamus aerodynamic model of insect-like flapping wings in hover. *WIT Transactions on State-of-the-art in Science and Engineering* **4**.
POLHAMUS, E. C. 1966 A concept of the vortex lift of sharp-edge delta wings based on a leading-edge-suction analogy. *NASA Technical Note TN D-3767*.

PURVIS, J. W. 1981 Analytical prediction of vortex lift. *Journal of Aircraft* **18** (4), 225-230.

RAMESH, K., GOPALARATHNAM, A., EDWARDS, J. R., OL, M. V. & GRANLUND, K. 2013 An unsteady airfoil theory applied to pitching motions validated against experiment and computation. *Theoretical and Computational Fluid Dynamics* **27** (6), 843-864.

RINALDI, S. & PAÏDOUSSIS, M. P. 2012 Theory and experiments on the dynamics of a free-clamped cylinder in confined axial air-flow. *Journal of Fluids and Structures* **28**, 167-179.

RYU, J., PARK, S. G., KIM, B. & SUNG, H. J. 2015 Flapping dynamics of an inverted flag in a uniform flow. *Journal of Fluids and Structures* **57**, 159-169.

SADER, J. E., COSSÉ, J., KIM, D., FAN, B. & GHARIB, M. 2016a Large-amplitude flapping of an inverted flag in a uniform steady flow – a vortex-induced vibration. *Journal of Fluid Mechanics* **793**, 524-555.

SADER, J. E., HUERTAS-CERDEIRA, C. & GHARIB, M. 2016b Stability of slender inverted flags and rods in uniform steady flow. *Journal of Fluid Mechanics* **809**, 873-894.

SANE, S. P. 2003 The aerodynamics of insect flight. *Journal of Experimental Biology* **206** (23), 4191-4208.

SHELLEY, M. J. & ZHANG, J. 2011 Flapping and bending bodies interacting with fluid flows. *Annual Review of Fluid Mechanics* **43**, 449-465.

SHOELE, K. & MITTAL, R. 2016a Energy harvesting by flow-induced flutter in a simple model of an inverted piezoelectric flag. *Journal of Fluid Mechanics* **790**, 582-606.

SHOELE, K. & MITTAL, R. 2016b Flutter instability of a thin flexible plate in a channel. *Journal of Fluid Mechanics* **786**, 29-46.

STAUBLI, T. & ROCKWELL, D. 1989 Pressure fluctuations on an oscillating trailing edge. *Journal of Fluid Mechanics* **203**, 307-346.

TANG, C., LIU, N.-S. & LU, X.-Y. 2015 Dynamics of an inverted flexible plate in a uniform flow. *Physics of Fluids* **27** (7).

TANG, L. & PAÏDOUSSIS, M. P. 2007 On the instability and the post-critical behaviour of two-dimensional cantilevered flexible plates in axial flow. *Journal of Sound and Vibration* **305** (1-2), 97-115.

TAVALLAEINEJAD, M., LEGRAND, M. & PAÏDOUSSIS, M. P. 2020 Nonlinear dynamics of slender inverted flags in uniform steady flows. *Journal of Sound and Vibration* **467**, 115048.

TAVALLAEINEJAD, M., PAÏDOUSSIS, M. P. & LEGRAND, M. 2018 Nonlinear static response of low-aspect-ratio inverted flags subjected to a steady flow. *Journal of Fluids and Structures* **83**, 413-428.

TAYLOR, G. I. 1952 Analysis of the swimming of long and narrow animals. *Proceedings of the Royal Society (London)* A **214** (1117), 158-183.

THEODORSEN, T. 1949 General theory of aerodynamic instability and the mechanism of flutter. *NACA Technical Report, No. 496*.

YAN, Z., TAHA, H. E. & HAJJ, M. R. 2014 Geometrically-exact unsteady model for airfoils undergoing large amplitude maneuvers. *Aerospace Science and Technology* **39**, 293-306.

YU, Y., LIU, Y. & CHEN, Y. 2017 Vortex dynamics behind a self-oscillating inverted flag placed in a channel flow: Time-resolved particle image velocimetry measurements. *Physics of Fluids* **29** (12), 125104(1-15).