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

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 materializes, 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.


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
Email address for correspondence: michael.paidoussis@mcgill.ca (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 D f L= p h; f and 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 UL= 10 4 10 5 , 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 D 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 D 250). For Re < 50, the stretched-straight and flapping modes did not arise. For 100 Re 250, however, the stretchedstraight, 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 A D H=L D 1 and 0:31 flapping in a high Reynolds number axial flow (Re 10 5 ). They obtained excellent agreement with the experimental results reported by Kim et al. (2013), particularly on the time history of the flag freeend 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 D 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 Throughout the present paper, including the Introduction, the mass ratio is defined as fluid mass to structure mass. Some other authors use the inverse definition. 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 ! 1); 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 D 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 D 2, W being the width of the test-section) in the Re D 4660 7600 flow regime. The inverted flag had an aspect ratio A D 3 and mass ratio ' 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 A=L 1:2) about the buckled shape in the 5100 Re 6000 flow regime. The other two modes correspond, respectively, to large-amplitude flapping (A=L 1:8) in the 6000 < Re 7122 flow regime, and large-amplitude buckling beyond Re D 7400. The flapping motion was found to undergo a hysteresis loop when the Reynolds number was decreased from Re D 7600 down to Re D 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 D 1 and at Re D 3 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 . 40; refer to Blevins 1990 andPaï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. 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 A 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 2 OE0: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 D 2 sin 2˛f or rigid, rectangular, very low aspect ratio plates,˛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 A D 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 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 saddlenode 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. (2020Tavallaeinejad et al. ( , 2018 developed nonlinear theoretical models to explore the static and dynamic responses of flexible, low aspect ratio (mostly A < 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 (A 0:5), relatively heavy ( 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 A was increased from A D 0:1 to A D 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 D 0:1 to D 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 A D 1:25, taper ratio 0 =H 1:75 ( and H being the height Study Approach A Re Kim et al. (2013) Experimental (air) 1:0 1:3 0:303 0:4 10 4 10 5 Experimental (water) 1:3 2:0 166:7 250 10 4 10 5 Ryu et al. (2015) Computational (  at the free and clamped ends, respectively) and two different thicknesses, yielding D 0:40 and D 0:61 (presumably for the rectangular sheets), in Re 10 4 10 5 flow. They found that as the taper ratio (or shape parameter) =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 postcritical 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 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 2 OE0:1 5:0, but the results may not match all features present in the experiments of Kim et al. (2013). 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.

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.

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. A ! 1); 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 p , p , E, and , respectively; 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 .OI XZ/, with the X.e X / and Z.e Z / axes being in the axial and transverse directions, Figure 2: (Left) The flag is idealised as an inextensible (i.e. e D 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 .X; t/. (Right) Joukowski's conformal mapping of the plate into a circle of radius L=4 (i.e. N X N Z-plane); q r and q are the radial and circumferential velocities on the circle in the N X N Z-plane, respectively. The impermeability condition is satisfied by considering a set of sources and sinks of strength on the top and bottom halves of the circle, respectively. The induced velocity due to a single bound vortex of strength and its image (i.e. the wake vortex of strength ) at an arbitrary point on the circle satisfies the Kutta condition.
respectively; (ii) a curvilinear coordinate system .OI 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 .X; t/, the curvature .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 .X; t/ D @ X w.X; t/ 1 C e.X; t/ ; cos .X; t/ D 1 C @ X u.X; t/ 1 C e.X; t/ ; .X; t/ D @ X .x; t/ 1 C e.X; t/ ; (2.1) where @ X @=@X denotes the first spatial derivative. The inextensibility assumption reads e.X; t/ D 0, which yields .1 C @ X u.X; t// 2 C .@ X w.X; t// 2 D 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 2008or 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 .x; t/. Hence, velocity components of an element may be formulated as 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 In the remainder of this paper, @ is defined as the first derivative with respect to , @ denotes the second derivative, and so on. directions may be written as In this paper, the governing equations are derived in a Hamiltonian framework. The extended Hamilton's principle may be expressed as with T .t/ and V.t/ denoting the kinetic and potential energies of the system, respectively, and ıW f .t/ and ıW d .t/ the virtual work by the fluid forces acting on the inverted flag and damping, respectively; here ı denotes the variational operator.

Kinetic and potential energies
In this study, the centreline rotation of the inverted flag .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 =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 .x; t/.
The kinetic energy of the system for translational and rotational motions of the inverted flag may be written as P .x; t/ 2 dx: (2.5) The potential strain energy of the inextensible inverted flag in terms of the rotation angle is in which D D Eh 3 =.12.1 2 p // 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 where denotes the viscoelastic damping coefficient. The virtual work done by the fluid-related force normal to the flag is given by 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.

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 .x; z; t/, which satisfies the Laplace equation (i.e. r 2 D 0, with r being a first order spatial differential operator). Following Theodorsen's approach, the perturbation velocity potential may be divided into two parts: D nc C c , where c is the circulatory component, fulfilling the Kutta condition, and 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 andYan et al. 2014), 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 P .x; t/ D P nc .x; t/ C P c .x; t/: (2.10) The non-circulatory (or added mass) component, represented by P 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 P 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).

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 D L=4, through the following standard conformal mapping illustrated in figure 2: Given that N X D L=.4 cos / and N Z D L=.4 sin /, the relation between points on the flag and the circle may be found using (2.11) as x D .1 C cos /=2.
Considering a continuous distribution of source-sink of strength .; t/ and .; 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 .; t/ D 4V n .; t/ sin (refer to Bisplinghoff et al. 2013).
Considering an infinitesimal element of the source sheet .L=4/ d' on the circle at D ', along with its sink counterpart at D ', and integrating over the half circle, the non-circulatory component of the resultant tangential velocity may be formulated as Assuming that nc u .L=4; 0; t/ D 0 at the leading edge (i.e. D 0), the non-circulatory perturbation velocity potential in the N X N Z plane at r D L=4 may be found, Substitution of (2.13) into (2.9), changing the variables by making use of x D .1 C cos /=2, and using the chain rule, the pressure distribution over the flag due to the non-circulatory contribution of the flow in the N X N Z plane at r D L=4 may be written as (2.14)

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. D ); 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, Subsequently, Theodorsen obtained the induced velocity due to a single bound vortex of strength and its image (i.e. the wake vortex of strength ) at an arbitrary point on the circle and formulated the circulatory component of the tangential velocity as Applying the Kutta condition at the trailing edge (i.e. D ) 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 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/ ! 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.

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 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 N X D q nc sin and q N Z D q nc cos , the leading edge suction force at r D L=4 is Although the suction analogy does not model the actual flow field, it has been shown to be 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)   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 & Zbikowski 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.

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 (2.20) with the total quasi-steady pressure formulated as where ı 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 D .1 C cos /=2 yields the virtual work done by the fluid-related force normal to the flag.

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 .x; t/: Utilizing the dimensionless parameters in which f is the frequency of flapping and D p p h=DL 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:

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.

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/ 1, u.x; t/ 0), equation (2.24) may be linearized as The series expansion 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 ' 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 R p C C P p C K p D 0 (see Eq. (A 7)). Defining the following partitioned vector and matrices of order 2M : the problem reduced into the first-order form B P r C Er D 0. Seeking periodic solutions in the form of r.t/ D Q r exp.i!t/, where Q r is the modal amplitude and ! is the complex frequency, one can obtain the following eigenvalue problem

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.!/ and Re.!/, respectively. It is recalled that Im.!/ corresponds to the damping ratio, and Re.!/ 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.!/ < 0); refer to Païdoussis (2014) for more details.
In this study, six modes are employed for numerical solutions (i.e. M D 6); this is sufficient for obtaining converged solutions. The other system parameters utilized in the simulations are as follows: L D 10 cm, h D 0:1 mm, D 1200 kg m 3 , f D 1:2 kg m 3 , D D 2454 N cm 2 and s D 0:0002. Figure 4 is the Argand diagram for an inverted flag with D 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.!/-axis. One of the solution branches associated with this mode eventually vanishes altogether (i.e. Re.!/ D 0, Im.!/ D 0) at … D 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 (… D 7:58), the negative and positive branches of the first mode loci coincide and leave the imaginary axis, while Im.!/ < 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 form 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 … D 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 … D 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 2 OE0: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.

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 pseudoarclength continuation technique. The Gear backward differentiation formula (BDF) is also used to obtain the time histories of oscillation.

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 ‰ i .x/ is selected to approximate the rotation of the flag (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 ‰ 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: with p and f p denoting the dimensionless counterparts of 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).

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 …. 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 D 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 … D 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 … D 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 … 2 OE1:70 1:82, and the amplitude of oscillation increases with …, 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 … D 1:80, all indicating a periodic motion; the dimensionless frequency of oscillation is f D 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 … D 1:82 and … D 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 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 … D 1:65, all displaying a periodic motion; the dimensionless frequency of oscillation is f D 0:52. This regime has been characterized in a number of references (see, for instance, Gurugubelli & Jaiman 2015, Shoele & Mittal 2016aand 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 …. As depicted, f sharply decreases with increasing … in the deformed-flapping regime, while it increases only  slightly in the large-amplitude flapping regime as … 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 … 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 D 0:2 (which is the reciprocal of the dimensionless mass ratio used by Goza et al.,  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 D 200), as well as high-Re simulations by Gurugubelli & Jaiman (2019) (Re D 30 000), 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 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. 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.

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 is present only in the time-dependent terms. Thus, it is expected that 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 , 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 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 ¤ 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 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 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, … 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 form the periodic solution around the deflected equilibrium (i.e. at … C ) goes backward and folds at a flow velocity (i.e. … D ) 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 D 5:0 ( figure 12(d)), two stable solutions coexist for … 2 OE1:10 1:36: the trivial solution and a large-amplitude periodic one. For … D 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 f D L=U , is used to non-dimensionalize the flapping frequency f . Thus, the reduced frequency is defined as f R D f L=U D f p =…. As seen in the figure, for all cases studied, f R decreases monotonically with increasing …. For heavier flags (and thus lower mass ratios), the reduced frequency decreases approximately linearly with ….
As also seen in figure 13, the reduced frequency decreases as is decreased. For example, for the heaviest flag with D 0:1, the reduced frequency is f R D 0:09 at … D 2:0, while for the lightest flag with D 5:0, f R D 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 R . 0:2, for D 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 R D 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.
Refer to figure 6 for definitions of … C and … D . The notation adopted here can be linked to that used in the first footnote in ÷ 2.3.2 as f R D k= .

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 0 , as illustrated in figure 14. Note that 0 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 D U e X D U cos 0 and U Z D U e Z D U sin 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/ D . P w U sin 0 / cos . P u C U cos 0 / sin ; (4.3a) V .x; t/ D . P u C U cos 0 / cos C . P w U sin 0 / sin : (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 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 D 1:0, and 0 D 0:1 5:7 ı 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 0 D 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 … D 1:8) and becomes unstable. By tracing the unstable solution, a second saddle-node bifurcation arises at … D 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 D 1:0 and 0 D 0 (see figure 10).
The bifurcation diagrams shown in figure 16 for 0 D 0:2 illustrate the effects of the resulting symmetry breaking on the response of the system for various mass ratios (cf. figure 12 for 0 D 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 0 D 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 largeamplitude flapping is more plausible for the inverted flag with a low mass ratio, e.g. D 0:1, the flag with the higher mass ratio D 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 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 0 D 0:2 is shown in figure 17. Here, the reduced frequency for both small-and largeamplitude flapping is shown. The branch in the low-… region corresponds to small-amplitude flapping, while that appearing in the high-… 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  transformed into large-amplitude flapping. Similarly to what was found for the flag with zero 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 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 0 ¤ 0 was gradual, as opposed to the case with 0 D 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 0 is increased. Furthermore, varying 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).  Goza et al. (2018) for M D 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 D 20.

Comparison
A quantitative comparison is made in Table 2 where thresholds of divergence, deformedflapping 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 D 0:2 at … D 2:0, f R D 0:13 is obtained by our model, while f R ' 0:09 is obtained by Goza et al. (2018) at Re D 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 D 200 and Re D 20) suggests that the overall dynamical behaviour and characteristic features of inverted flags is insensitive to the Reynolds number. This has been 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)).   shown by Shoele & Mittal (2016a) who performed direct numerical simulations at small Reynolds numbers and reported that the frequency and the amplitude of large-amplitude flapping remain almost unchanged in the range 100 Re < 1000. Tang et al. (2015) explored the Reynolds number effect in the range 100 Re 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 D 30 000 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.

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 qualitatively 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 D 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 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). the model to accurately predict the dynamics of the flag, for example, when very large flapping 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.
.cos ' C 1/ j 1 .1 cos '/ d'; where †1 n D m 2n X rD0 m 2n r ! cos r # sin k 1 # sin k 2 ; †2 n D m 2nC2 X rD0 m 2n C 2 r ! cos r # sin kC1 # sin .k C 2/ 2 ; (A 6) and k D m 2n r. Hence, integrations over ' are performed exactly and the pressure distribution over the inverted For the nonlinear analysis, the following comparison functions are defined by their polynomial successors to discretize the nonlinear equation of motion (2.24): recalling that i is the ith root of the transcendental equation 1 C cosh.x/ cos.x/ D 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 .x; t/ D @ x w.x; t/ holds for the linear equation of motion of the flag; hence ‰ i .x/ D @ 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. ‰ i .x/ D P J i j D1 a ij jx 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 .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, . 1/ r .2r/Š .s; t/ 2r ds C 1 . 1/ r .2r C 1/Š .s; t/ 2rC1 ds ; where R is the order of truncation. Similarly to the procedure detailed earlier, inserting ‰ i .x; t/ D P M iD1 q i .t/ P J i j D1 a ij jx j 1 into (A 9), 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 (A 9)