Self-induced flapping dynamics of a flexible inverted foil in a uniform flow

We present a numerical study on the self-induced flapping dynamics of an inverted flexible foil in a uniform flow. A high-order coupled fluid–structure solver based on fully coupled Navier–Stokes and nonlinear structural dynamic equations has been employed. Unlike a conventional flexible foil flapping where the leading edge is clamped, the inverted elastic foil is fixed at the trailing edge and the leading edge is allowed to oscillate freely. We investigate the evolution of flapping instability of an inverted foil as a function of the non-dimensional bending rigidity, $K_{B}$ , Reynolds number, $\mathit{Re}$ , and structure-to-fluid mass ratio, $m^{\ast }$ , and identify three distinct stability regimes, namely (i) fixed-point stable, (ii) deformed steady and (iii) unsteady flapping state. With the aid of a simplified analytical model, we show that the fixed-point stable regime loses its stability by static-divergence instability. The transition from the deformed steady state to the unsteady flapping regime is marked by a flow separation at the leading edge. We also show that an inverted foil is more vulnerable to static divergence than a conventional foil. Three distinct unsteady flapping modes have been observed as a function of decreasing $K_{B}$ : (i) inverted limit-cycle oscillations, (ii) deformed flapping and (iii) flipped flapping. We characterize the transition to the deformed-flapping regime through a quasistatic equilibrium analysis between the structural restoring and the fluid forces. We further examine the effects of $m^{\ast }$ on the post-critical flapping dynamics at a fixed $\mathit{Re}=1000$ . Finally, we present the net work done by the fluid and the bending strain energy developed in a flexible foil due to the flapping motion. For small $m^{\ast }$ , we demonstrate that the flapping of an inverted flexible foil can generate $O(10^{3})$ times more strain energy in comparison to a conventional flexible foil flapping, which has a profound impact on energy harvesting devices.

The fluid-elastic instability with such a configuration can be easily seen by blowing air over a thin piece of paper or a leaf and in waving of flag-like structures. A flexible foil clamped at the leading edge and with the trailing edge free to oscillate, from here on referred to as a conventional foil, is stabilized by the viscous drag and structural stiffness, while the pressure force tends to produce the opposite effect. Applications of this phenomenon include the implementation of new surgical methods (Huang 1995), increasing the speed of paper printing (Watanabe et al. 2002a,b), nuclear plate assemblies (Guo & Paidoussis 2000), flow control devices (Lucey 1998;Jaiman, Loth & Dutton 2004) and botanical applications (Huang, Rominger & Nepf 2011). A system of flexible foils has been also proposed as a means of harvesting energy, which can be utilized to generate electric energy (Allen & Smits 2001;Tang, Paidoussis & Jiang 2009). In this view, Michelin & Doare (2013) recently proposed a theoretical framework to study the optimal position of piezoelectric patches on an elastic foil to maximize the efficiency.
The flapping of a conventional flexible foil has been extensively studied using various theoretical, experimental and numerical methods.  and Shelley & Zhang (2011) have provided an extensive literature review on this topic. Early studies on the flapping instability (Dowell 1966(Dowell , 1970 were based on linear instability analysis using the Galerkin decomposition in space to predict the critical velocity above which the flexible foil exhibits flapping motion. These studies largely focused on supersonic flow and were motivated by aerospace applications. Kornecki, Dowell & O'Brien (1976) have theoretically investigated the flapping instability for subsonic flows by introducing a circulatory fluid loading to satisfy the Kutta condition at the trailing edge. Extending this theoretical analysis, Argentina & Mahadevan (2004) included the effects of viscous drag induced tension and finite size in their theoretical model. To analyse the influence of trailing and leading edge boundary conditions on the flapping instability, Guo & Paidoussis (2000) performed a linear stability analysis in Fourier space.
More recently, Shelley, Vandenberghe & Zhang (2005), with the aid of a simplified analytical model and water-tunnel experiments, have shown that a flexible foil in a heavy fluid can experience strong added-mass effects, which can stabilize the foil. By employing a three-dimensional potential flow with two-dimensional foil deflections, Eloy, Souilliez & Schouveiler (2007) studied the effects of aspect ratio on the stability and dynamics of a flexible foil. In this study, the authors observed that flexible foils with small aspect ratios tend to be more stable compared with foils with larger aspect ratios. Eloy, Kofman & Schouveiler (2012) further analysed the hysteresis effects observed experimentally and attributed these effects to the surface defects of the flexible plate. In addition to linear/nonlinear stability analyses, Tang & Paidoussis (2007) and Michelin, Llewellyn Smith & Glover (2008) have presented reduced-order models based on lumped vortex and unsteady point vortex models respectively to predict the post-critical behaviour of the flexible foil.
Since the experimental study by Taneda (1968), the flapping instability of a conventional foil has been the focus for a large number of experimental studies (Datta & Gottenberg 1975;Huang 1995;Zhang et al. 2000;Allen & Smits 2001;Watanabe et al. 2002b;Shelley et al. 2005;Eloy et al. 2008Eloy et al. , 2012. Zhang et al. (2000) utilized a soap-film flow to replicate a two-dimensional Navier-Stokes fluid and examined the flapping instability of a conventional foil. Two distinct dynamic regimes called a stretched-straight mode and a regular limit-cycle oscillation were observed. In this experiment, the authors noted a large region of bi-stability for a sufficiently large initial perturbation. Watanabe et al. (2002b) studied the influence of Flapping dynamics of a flexible inverted foil 659 bending rigidity and aspect ratio on the flapping instability for improving the speed of paper printing. In this study, the authors reported that the flapping amplitudes and frequencies increase with the flow velocity, and no boundary layer separation is observed along the foil. To understand the origin of the difference between the stability analysis and the experimental observations, Eloy et al. (2008) compared the stability analysis with the experimental observations. They observed that the three-dimensional effects become significant with increase in the aspect ratio of foils.
In parallel to the experimental studies, a large number of numerical studies have been performed on a conventional foil set-up using various methods. In particular, Zhu & Peskin (2002 simulated the soap-film experiment of Zhang et al. (2000) using the immersed boundary method. Although the Reynolds number of the direct numerical simulations was O(10 2 ) times lower than the experiment, the simulations predicted the dynamic regimes observed experimentally (Zhang et al. 2000). In this numerical study, it was found that the flexible foil did not experience a flapping instability when the mass of the flexible foil was neglected. The soap-film experiment of Zhang et al. (2000) was simulated by Farnell, David & Barton (2004) using a loosely coupled fluid-structure solver. Alben (2009) reported a numerical method to simulate the flapping dynamics in the limit of infinite Reynolds number. In this numerical method, a fluid flow was modelled as an inviscid flow consisting of a thin vortex sheet, and the viscous effects were captured through a thin vortex sheet. Using this numerical method, Alben (2009) analysed the transition from regular limit-cycle flapping to chaotic flapping. Connell & Yue (2007) and Liu, Jaiman & Gurugubelli (2014) presented a comprehensive numerical analysis on the effects of the mass ratio, m * , and Reynolds number, Re, on the flapping dynamics; here, Connell & Yue (2007) used a strongly coupled finite difference formulation and Liu et al. (2014) developed a numerical formulation based on the high-order finite-element method using arbitrary Eulerian-Lagrangian (ALE) coordinates. In both of these studies, three distinct regimes were identified: (i) a steady fixed-point response for small m * , (ii) a periodic limit-cycle flapping for intermediate m * and (iii) a chaotic flapping regime found for large m * . These results showed a good agreement with the analytical findings of critical mass ratio.
In the present work, we gain our motivation from the ability of self-sustained flapping of an elastic foil to extract energy and generate electric power (Tang et al. 2009;Akcabay & Young 2012) from the vast and untapped source of energy available in the form of ocean currents and tidal flows. A large number of energy harvesting models (Allen & Smits 2001;Tang et al. 2009;Akcabay & Young 2012;Michelin & Doare 2013) have been proposed based on the idea of converting fluid kinetic energy into strain energy through the flapping and bending process of a conventional foil. This strain energy in turn can be converted to electric potential using capacitive, conductive or piezoelectric methods. To optimize the piezoelectric energy harvesting, Michelin & Doare (2013) performed coupled piezoelectric fluid-structure simulations of a conventional foil. In this study, the authors reported that the efficiency can be up to 10-12 % for a structure-to-fluid mass ratio of 20, it decreases to less than 1 % for m * ≈ 1 and becomes relatively small for m * 1. This leads to the requirement to develop new energy harvesting devices with an alternative configuration to extract the maximum possible kinetic energy from the fluid flow even for m * 1, which typically represents the case of a flexible foil flapping in flowing water.
All of the aforementioned literature corresponds to a conventional flexible foil, whose dynamics is comparatively well understood. The objective of this paper is to study the flapping mechanism of a flexible foil whose trailing edge is clamped and whose leading edge is free to oscillate; hereafter this will be referred to as an inverted foil. An inverted foil induces a remarkably different flapping dynamics with respect to the conventional flapping. Guo & Paidoussis (2000) employed a linear stability theory and observed that an inverted flexible foil becomes unstable for non-zero flow velocities when the viscous damping effects are neglected. This observation clearly contradicts the physical observations, where an inverted flexible foil can exhibit stability for non-zero fluid velocity. Buchak, Eloy & Reis (2010) performed an interesting experimental study by placing a book in a wind tunnel and analysed the clapping interaction among a stack of papers. This experimental study closely represented the dynamics of an inverted foil flapping. More recently, Kim et al. (2013) investigated the dynamics of an inverted foil in both wind and water tunnels. They reported the existence of three flapping regimes: two quasisteady modes and a limit-cycle flapping mode. They observed that an inverted foil was more prone to flapping instability compared with a conventional foil, i.e. the critical velocity of an inverted foil was lower than that of a conventional foil. However, the physical phenomenon underlying this observation requires some attention through nonlinear viscous simulations of fluid-structure interactions.
Here, we perform two-dimensional numerical simulations to assess the flapping dynamics of an inverted flexible foil configuration. For this purpose, we adopt a high-order fluid-structure interaction solver based on the combined field with explicit interface (CFEI) formulation proposed by Liu et al. (2014). This solver captures the nonlinearities of the problem coming from the Navier-Stokes equation and the geometrically nonlinear structural dynamics. Through a series of coupled numerical simulations, we investigate the evolution of flapping instability in an inverted foil as a function of the non-dimensional bending rigidity K B , m * and Re. We characterize the underlying physical mechanism for the onset of flapping instability and present stability phase diagrams of K B versus Re and m * through direct numerical simulation (DNS) results and simplified analytical solutions. Furthermore, we conduct a systematic parametric study to understand the behaviour of the post-critical flapping dynamics and its dependence on the relevant non-dimensional parameters. As a function of decreasing K B , we observe three distinct flapping regimes: inverted limit cycle, deformed flapping and flipped flapping. We visualize the evolution of vortex patterns as a function of K B for these flapping response regimes. We finally analyse the ability of inverted foil to convert the available fluid kinetic energy to structural strain energy during the flapping and bending process.
The content of the paper is structured as follows. Section 2 describes the problem of a flexible plate in a uniform axial flow and its governing equations. In § 3, we provide a brief description of the coupled fluid-structure solver used in this study. Section 4 presents a mesh convergence study and numerical verification of the coupled fluidstructure formulation. Finally, § 5 covers the detailed results and analysis of inverted foil, namely the onset of flapping instability, the effects of bending rigidity and mass ratio, the transition mechanism to the unsteady deformed-flapping mode, the regimes of flapping dynamics, the wake topology and the net energy transfer.

Problem statement
We consider a two-dimensional thin flexible foil Ω s interacting with an incompressible uniform axial flow Ω f (t). As shown in figure 1, the trailing edge (TE) of the flexible foil is clamped and the leading edge (LE) is free to perform flapping motion. Here, L and h represent the length and thickness of the foil respectively, U 0 is the magnitude of the uniform axial flow exciting the flapping instability, ρ f and ρ s are the densities of the fluid and structure, µ f denotes the fluid dynamic viscosity, α represents the leading edge angle of attack, E is Young's modulus and ν is the Poisson's ratio of the foil.
The Navier-Stokes equations governing an incompressible flow in an arbitrary Lagrangian-Eulerian reference frame are where u f = u f (x, t) and w = w(x, t) represent the fluid and mesh velocities defined for each spatial point x ∈ Ω f (t) respectively, f f is the body force applied on the fluid and σ f is the Cauchy stress tensor for a Newtonian fluid, written as where p is the fluid pressure, I denotes the second-order identity tensor and T represents the fluid viscous stress tensor. The deformation of the flexible foil is governed by the structural equation where u s = u s (z, t) is the structural velocity defined for each material point z ∈ Ω s , f s represents the external forces applied on the solid and σ s denotes the first Piola-Kirchhoff stress tensor. In this study, the structural stress tensor is modelled using the Saint Venant-Kirchhoff model (Antman 1995). The coupled system required to satisfy the velocity and traction continuity conditions along the fluid-structure interface is as follows: where n f and n s are respectively the outward normals to the deformed fluid and the undeformed solid interface boundaries, Γ 0 denotes the fluid-structure interface Γ (t) at time t = 0 and ϕ s is the displacement function that maps each Lagrangian point z ∈ Ω s to its deformed position at time t. Here, γ is any part of interface Γ 0 and ϕ s (γ , t) represents the corresponding fluid part over the interface Γ (t) at time t.

Coupled fluid-structure formulation
We employ a high-order variational method stable for very low mass ratios based on the CFEI formulation. In Liu et al. (2014), we have proved the numerical stability of the formulation for low structure-to-fluid mass ratio through energy estimates and demonstrated the ability of the scheme to simulate flapping dynamics. In this formulation, the solid positions and arbitrary Lagrangian-Eulerian mesh velocities are decoupled from the remaining variables (i.e. fluid velocity, pressure and structural velocity) and are handled explicitly. Consequently, the CFEI formulation only requires a linear system of equations to be solved for each time step. This step significantly decreases the size of the matrix that needs to be solved at each time step.
The weak form of the Navier-Stokes equations (2.1) and (2.2) can be written as Here, ∂ t denotes the partial time derivative operator ∂(·)/∂t, φ f ∈ H 1 (Ω f (t)) and q ∈ L 2 (Ω f (t)) are the smooth test functions in the Sobolev space for the fluid velocity and pressure respectively, Γ f n (t) represents the fluid Neumann boundary along which σ f (x, t) · n f = T f . The weak form of the structural dynamics equation (2.4) can be written as where φ s ∈ H 1 (Ω s ) denotes a smooth test function for the structural velocity. Similarly to the fluid domain, Γ s n represents the solid Neumann boundary and σ s (z, t) · n s = T s . The weak form of the traction continuity condition along the fluid-structure interface (2.5) is given by (3.4) In the above equation, one may observe that u f and u s are defined on different domains Ω f (t) and Ω s respectively, and the condition can be realized by considering a conforming mesh along the interface Γ 0 with ϕ s (z, t) being the position vector of the deformed solid. We next enforce the weak form of the traction continuity condition in (3.4) to combine (3.1)-(3.3) to construct the weak form of the coupled fluid-structure. For the weak form (Temam 2001) of the coupled fluid-structure formulation, let us introduce the finite-dimensional trial solution function space S along with the corresponding test function space V . The weak form of the formulation can be written as follows: (3.6) and the velocity continuity condition is enforced in the function space (Liu et al. 2014). As mentioned earlier, the main feature of the CFEI scheme is to explicitly determine the solid positions and the ALE-mesh displacements at the start of each time step. This feature enables us to decouple the ALE-mesh field from the remaining fluid-structure variables (u f , p and u s ). The solid position, ϕ s,n , for the nth time step can be determined using the second-order Adam-Bashforth method, The ALE-mesh nodes on the fluid domain Ω f (t) can be updated using a steady pseudo-elastic material model (Stein, Tezduyar & Benney 2003) given by where σ m is the stress experienced by the ALE mesh due to the strain induced by the interface deformation. Assuming that the ALE mesh behaves as a linearly elastic material, the stress experienced can be written as where η f denotes the ALE-mesh nodal displacement satisfying the boundary conditions (3.10) Here, Γ f (t) is the fluid domain boundary and Γ f (t)\Γ (t) denotes the non-interface fluid boundary; τ m is a mesh stiffness variable chosen as a function of the element size to limit the distortion of the small elements located in the immediate vicinity of the fluid-structure interface. The mesh stiffness variable τ m has been defined as τ m = (max i |T i | − min i |T i |)/|T j |, where T j represents the jth element on the mesh T.
The weak variational form in (3.6) is discretized in space using P n /P n−1 /P n isoparametric finite elements for the fluid velocity, pressure and solid velocity respectively. In this paper, we consider the stable P 2 /P 1 /P 2 isoparametric finiteelement meshes, which satisfy the inf-sup condition for well-posedness.

Numerical verification and convergence study
4.1. Numerical verification To verify the coupled solver, we simulate a conventional flexible foil flapping in a uniform axial flow and compare our findings with the theoretical model and DNS results (Connell & Yue 2007). The flapping dynamics of the conventional foil is strongly influenced by three critical non-dimensional parameters, namely m * , Re and K B (Shelley et al. 2005;Connell & Yue 2007), which are defined as Here, B represents the flexural rigidity defined as B = Eh 3 /12(1 − ν 2 ). The reduced velocity U r = √ m * /K B can also be used to quantify the dynamic effects of flapping. We consider a two-dimensional computational domain of size [22L × 10L], as shown in figure 2(a). In this computational domain, a flexible foil of length L and thickness h is placed along the x-axis with its leading edge clamped and centred at (0, 0). In figure 2(a), Γ To select a suitable mesh for the simulation of the conventional foil, a mesh convergence study is performed for two different meshes M1 and M2 consisting of 15 699 and 30 368 P 2 /P 1 high-order triangular elements. The critical non-dimensional parameters considered for this study are m * = 0.1, K B = 0.0001 and Re = 1000. A constant time-step size of t = 0.001 is employed with a characteristic time scale of L/U 0 . Figures 2(b) and 2(c) show a typical high-order undeformed fluid mesh and a closeup view of the boundary layer mesh around the flexible foil respectively. Table 1 summarizes the trailing edge tip-displacement statistics for the meshes M1 and M2. The values within the brackets represent the percentage difference in the numerical solutions with respect to the mesh M2. The maximum difference between the numerical solutions calculated for the meshes M1 and M2 is less than 1 %. Therefore, the mesh M1 will be used in the numerical verification study.
First, we verify the ability of the coupled fluid-structure solver to predict the transition boundary from a stable mode to a flapping mode for a conventional flexible foil. An expression for the critical K B below which an infinitely wide conventional flexible foil loses its stability can be written as Here, C T denotes the non-dimensional viscous tension coefficient which can be approximated to 1.328 with the aid of the laminar boundary layer theory (Connell & Yue 2007). The above expression can be obtained by solving the quadratic dispersion relation given by Connell & Yue (2007) for Im(ω) < 0. For low m * , Connell & Yue (2007) and Shelley et al. (2005) have shown that a conventional flexible foil loses its stability when the fundamental mode k x = 2π becomes unstable. As part of this verification, we compare the critical K B observed from our DNS results with the analytical solutions of (4.2) for a range of K B and Re at a fixed m * = 0.1. Figure 3 summarizes our DNS results for Re = 500, 1000 and 2000 in a stability phase diagram of a conventional foil. In this figure, the curve (--) represents the neutral curve between the stable and fundamental flapping modes, which is given by (4.2) for k x = 2π and m * = 0.1. The curve (---) represents the transition between the fundamental mode k x = 2π and the higher mode k x = 3π given by (4.2) for m * = 0.1. The stability phase diagram shows good agreement between our DNS and the analytical solution for Re 1000.
To prove the existence of the fundamental mode k x = 2π and the transition from the  fundamental mode k x = 2π to the higher mode k x = 3π, we plot the full-body profiles of the flexible foil when the trailing edge is at its extrema and when it is crossing the mean position in figure 4 for K B = 3.27 × 10 −4 and 1.8 × 10 −4 . From figure 4(a) at K B = 3.27 × 10 −4 , it can be inferred that the flapping modes range between k x = 2π and 5π/2. Since k x = 2π represents the least stable mode of k x ∈ [2π, 5π/2], the use of k x = 2π seems to be appropriate. For K B = 1.8 × 10 −4 , as shown in figure 4(b), the flexible foil exhibits flapping modes in the range k x ∈ [3π, 7π/2], and the explanation given for k x = 2π can be extended to k x = 3π. Next, we assess the accuracy of the coupled solver in simulating the flapping dynamics of a conventional foil for m * = 0.075, Re = 1000 and K B = 0.0001. In table 2, the key flapping properties for these cases are compared with the DNS results of Connell & Yue (2007) (2007) at Re = 1000 and K B = 0.0001.
solution (Connell & Yue 2007) and the current study. Despite different underlying structural formulations, fluid-structure coupling techniques and discretization methods, it can be observed that there is reasonable agreement among the flapping properties for both the mass ratios.

Mesh convergence of inverted foil
To examine the instability and post-critical nonlinear dynamics of an inverted flexible foil, we consider two computational domain configurations. The first computational domain will be used to study the initial development of instability in the inverted flexible foil. This computational domain is symmetrical about the foil centreline and is similar to the one presented earlier for the conventional foil, except that the trailing edge of the foil is clamped and the leading edge is free to oscillate. The converged finite-element mesh M1 defined for the conventional foil will be employed to simulate the evolution of flapping instability. Figure 5(a) shows a typical schematic of the second computational domain, which will be used to simulate the post-critical flapping dynamics of the inverted foil. In contrast to the first computational domain, this domain is asymmetrical about the foil centreline i.e. H u = H l . The size of the computational domain is [24L × 10L], with H u = 5.25L and H l = 4.75L. The inverted flexible foil is neutrally stable in its undeformed state and it can bend on either side of its centreline. In order to have a numerically repeatable solution, we have considered a non-symmetrical computational domain. To study the influence of domain size for K B = 0.2, Re = 1000 and m * = 0.1, two different domain widths, 10L and 40L, have been selected for assessment of the flapping response dynamics. The difference in the maximum tip displacement δ max has been found to be less than 1 %. The maximum blockage ratio, which is defined as δ max /(H u + H l ), has been observed to be approximately 8.5 % and 2.125 % for the two domains respectively. We next perform a detailed mesh convergence study for the asymmetrical computational domain using meshes M1, M2 and M3 consisting of 11 643, 23 197 and 46 104 P 2 /P 1 triangular elements. A typical undeformed high-order finite-element fluid mesh M3 and the corresponding central block consisting of the boundary layer mesh can be seen in figure 5(b). Figure 5(c) shows a closeup view of the deformed mesh around the flexible foil for three time instants. For the purpose of mesh convergence study, we consider a flexible foil with m * = 0.1 and K B = 10 −4 interacting with fluid flow at U 0 = 1 corresponding to Re = 1000. A constant time-step size t = 0.001 has been employed. Table 3 summarizes the results of the mesh independence test for the inverted foil. The values within the brackets indicate the percentage difference in the numerical solutions with respect to the mesh M3. The r.m.s. value of the cross-stream tip displacement is calculated using δ rms where δ * is the non-dimensional tip displacement defined as δ/L andδ * denotes the non-dimensional mean tip displacement. A similar definition has been used to calculate C rms l . Therefore, the mesh M2 will be employed to simulate the post-critical flapping dynamics of the inverted flexible foil.  this study. We first examine the evolution of the flapping instability as a function of K B for a fixed m * = 0.1 and Re = 1000 by gradually decreasing K B from 0.8 to 0.2. Figure 6 shows the time history of the cross-stream leading edge displacements for K B = 0.56, 0.55, 0.5475, 0.545 and 0.5425. From this figure, we can identify three distinct stability regimes: (i) fixed-point stable, (ii) deformed steady state and (iii) unsteady flapping regime. In the first regime, for K B 0.56, the flexible foil remains steady in its initial configuration with negligibly small leading edge cross-stream displacements. In the second regime, for 0.56 < K B < 0.5425, the flexible foil deforms slowly to achieve a steady state. This steady state deformation increases with a decrease in K B . In the unsteady flapping regime, for K B 0.5425, the flexible foil loses its stability, to perform regular sinusoidal oscillations with a constant frequency and amplitude.

Static divergence
The transition from the fixed-point stable to the deformed steady state regime can possibly be attributed to the static-divergence instability (Bisplinghoff, Ashley & Halfman 1957;Fung 1969). We consider a simplified model of the static divergence by introducing a small normal impulsive force P at the leading edge. As a result of this impulsive force, the flexible foil deforms with a leading edge cross-stream displacement δ y and an angle of attack α. The inverted flexible foil will experience the static-divergence instability when the aerodynamic moment acting on it exceeds the elastic restoration moment (Fung 1969). The aerodynamic moment can be evaluated by the lift force acting at the aerodynamic centre of the inverted flexible foil. The elastic restoration moment of the inverted flexible foil can be idealized as a spring connected at the leading edge. Therefore, the condition for the static divergence will be kδ y L < 1 2 ρ f U 2 0 C l (α, Re)Lx ac , where k is the spring constant of an elastically mounted foil, C l is the lift coefficient expressed as function of Re and α, and x ac denotes the distance between the aerodynamic centre and the trailing edge. From the DNS results corresponding to the deformed steady state regime, we construct an empirical formulation for C l as a function of Re and α. Figure 7(a) shows the relationship between C l and α obtained from the DNS simulations for a typical Re = 400. From this figure, it can be noticed that C l is directly proportional to α and can be represented as where (∂C l /∂α)(Re) denotes the slope of the curve in figure 7(a). Figure 7(b) shows the variation of ∂C l /∂α for a range of Re, and the curve (--) represents the best-fit curve given by the equation where a = −5.524 × 10 −7 , b = 1.139 × 10 −3 and c = 2.492. Therefore, (5.1) will be On substituting x ac = 3/4L, k = 3EI/L 3 , δ y = PL 3 /(3EI) and α = PL 2 /(2EI) (Crandall, Dahl & Lardner 1995) in (5.4) and simplifying we obtain EI ρ f U 2 0 L 3 < 3 16 (aRe 2 + bRe + c). (5.5) The expression on the left-hand side represents K B from (4.1), and the critical K B for transition from the fixed-point stable to the deformed steady state regime will be (K B ) st = 3 16 (aRe 2 + bRe + c).
(5.6) 5.1.2. Stability regimes Figure 8 summarizes the stability regimes of an inverted flexible foil for a range of Reynolds numbers. In this figure, the curve (--) depicts (K B ) st from (5.6) and the curve (---) represents the transition from the deformed steady state to the unsteady flapping regime given by the empirical relation We will further discuss the effects of mass ratio in § 5.3.

Conventional foil versus inverted foil
To understand why a conventional foil configuration does not experience a static divergence, we perform two sets of simulations on a fixed rigid flat plate at Re = 1000. In the first case, we rotate the rigid plate about the trailing edge to provide a small angle of attack, 5 • at the leading edge. For the second case, a small angle of attack, −5 • , is provided at the leading edge by rotating about the leading edge. A symmetrical computational domain with H u = H l = 5L has been used for this study. Figure 11 shows the pressure distribution around the rigid plate for the two angle of attack cases. The first case with the positive angle of attack typically represents an inverted rigid plate experiencing a static divergence, whereas the second case represents a conventional rigid plate in a static divergence. For the first case, figure 11(a) shows that the lower surface is at a higher pressure than the upper surface. Hence, the aerodynamic lift acts in the direction of rotation and the rigid plate can undergo a static-divergence instability when the lift force crosses the critical restoring bending force (Fung 1969). On the contrary, from figure 11(b) for the negative angle of attack, it can be observed that the nett aerodynamic lift acts in the direction opposite to the angle of rotation and has a stabilizing effect. Therefore, a cantilever beam clamped at the leading edge is less prone to a static-divergence instability.   Figure 12(a) shows the fixed-point stable foil profile in its undeformed state for K B = 0.56. For K B ∈ [0.5425, 0.2], the flexible foil experiences the onset of flapping instability which develops into large-amplitude periodic oscillations, and this regime is defined as the inverted limit-cycle oscillations (LCO). Figure 12(b-d) shows typical full-body profiles for K B = 0.5425, 0.4 and 0.2 over one oscillation cycle. During the initial stages of the inverted-LCO regime for K B ∈ [0.50, 5425] the flexible foil does not perform flapping about the mean position, as shown in figure 12(b). However, for K B ∈ [0.5, 0.2] the flexible foil performs flapping about its mean position. Moreover, from these figures, it can be observed that the inverted-LCO mode is predominantly the first mode and does not exhibit the typical necking phenomenon that is observed in a conventional foil (Huang 1995;Liu et al. 2014). By decreasing the K B value, the flexible foil completely deforms to one direction and performs relatively small peak-to-peak amplitude oscillations as shown in figure 12(e). Due to the small peak-to-peak amplitudes, this regime has been described as a deformed mode by Kim et al. (2013). Since the peak-to-peak amplitudes for this regime are of the same order as the amplitudes for a conventional foil flapping, we describe this regime as the deformed-flapping regime. By further decreasing the K B value, due to the low flexibility, part of the foil curves around the fixed edge to form a semicircular arc and the remaining part of the foil becomes parallel to the free stream velocity, as shown in figure 12(i). Figure 12(f -i) shows the full-body profiles for the inverted flexible foil for a range of K B < 0.1. The full-body profile shown in figure 12(i) closely resembles the full-body profile of the conventional foil shown in figure 13(a) for K B = 3.27 × 10 −4 , Re = 1000 and m * = 0.1. Additionally, the necking phenomenon that is generally observed in the case of a conventional foil flapping has been observed in this regime. The maximum Reynolds numbers with respect to the semicircular arc diameter, Re δ , is approximately 380 for K B = 0.05, which signifies that our current two-dimensional simulations might not be able to capture the three-dimensional flow structures.

Effects of bending rigidity
To further compare the flapping phenomenon in the flipped-flapping regime to the conventional foil, we plot the phase relation between the streamwise and the cross-stream tip displacements, as shown in figure 14. The phase plot for the flipped-flapping regime at K B = 10 −3 forms an eight-shaped figure (figure 14a), which is similar to the phase relationship shown in figure 13(b) for the conventional foil with K B = 3.27 × 10 −4 , Re = 1000 and m * = 0.1. The phase plot for the inverted LCO presents a unique parabolic profile, i.e. the leading edge follows the same path during the upward and downward strokes, which can be seen in figure 14(c) for K B = 0.2. Moreover, no distinct phase plot has been observed for K B corresponding to the flipped-flapping regime, as shown in figure 14(a,b) for K B = 0.01 and 0.001 respectively. Figure 15 shows the time history of streamwise (---) and cross-stream (--) displacements along with the cross-stream amplitude-frequency spectrum. In this figure, q , represents the time instants for which vorticity contours will be presented in § 5.6. The leading edge tip displacements for the inverted-LCO regime exhibit a sudden rise when K B is decreased from 0.545 to 0.4, as shown in figure 15(a,b). Figure 15(c) shows the leading edge tip-displacement history for the deformedflapping regime at K B = 0.1. In this regime, the leading edge remains steady for a significant time period, which is followed by non-periodic oscillations with relatively small peak-to-peak amplitudes ranging between 0.2 and 0.4. Figure 15(d-g) shows the tip-displacement response history of the flipped-flapping regime. The flapping amplitudes for this regime are significantly smaller than those observed in the inverted-LCO regime and are approximately equal to those observed in the conventional foil. However, the flapping frequencies are smaller than those observed in the conventional foil flapping. The transition from the deformed-flapping to the flipped-flapping regime is predominantly characterized by two frequency oscillations, as shown in figure 15(d). The first low harmonic corresponds to the flapping frequency and the second frequency may be attributed to an additional higher-harmonic frequency of shear layer roll-up over the semicircular configuration. According to Strouhal's law, it may be estimated as f H = St(U 0 /D foil ), where D foil is the effective diameter of the semicircular body and St is the Strouhal number. The higher harmonics 0.6838 and 1.025 appearing in figure 15(d,e) approximately coincide with the Strouhal frequencies corresponding to the vortex shedding of the semicircular body. In addition, we also observe a complex vortex-foil interaction and a reversed flow inside the deformed cavity-like region. A detailed investigation of these phenomena is beyond the scope of this study. The wake topology will be further illustrated in § 5.6 with the aid of the vorticity evolution in a flapping period. The high-harmonic frequency gets further damped as K B is lowered from 0.01 to 0.001, and the flapping response for K B = 0.001 exhibits a single distinct frequency (figure 15h). Figure 16 summarizes the effect of K B on the maximum leading edge cross-stream tip displacement. For large K B values, the flexible foil remains stable. A sharp rise in the maximum tip displacement can be observed for K B ∈ [0.3, 0.5]. The maximum flapping amplitudes stabilize to an approximate value of 0.85, which is similar to the experimental observation made by Kim et al. (2013).

Effects of mass ratio
To analyse the effects of mass ratio on the inverted flexible foil, we perform a series of numerical simulations for m * ∈ [0.1, 10]. We first examine the effects of m * on the inverted-LCO regime for a fixed K B = 0.4 and Re = 1000. Figure 17(a) shows the leading edge cross-stream tip-displacement history for m * = 0.5 and 8. From this figure, it can be seen that the flapping frequency is strongly influenced by the increase in m * , while the cross-stream flapping amplitudes are weakly affected. It may be noted that with increase in m * , the unsteady flapping instability takes a greater time to develop due to the inertial effects. Figure 17(b,c) summarizes the effect of m * on the flapping amplitudes and frequencies.  Figure 18 shows the leading edge cross-stream displacement histories for m * = 0.1, 0.5, 1 and 2. From figure 18(a,b), we observe regular periodic oscillations for m * = 0.1 and 0.5, which transform into non-periodic oscillations by increasing m * to 1. These non-periodic oscillations are characterized by variable amplitudes and frequencies. Figure 18(c,d) shows the non-periodic oscillations for m * = 1 and 2. A similar observation of transition from regular periodic to nonperiodic oscillations with increase in m * was made by Connell & Yue (2007) for conventional foils. The mass ratio m * can also affect the transition from the inverted-LCO to the deformed-flapping regime. Figure 19 shows the leading edge cross-stream displacements for m * = 0.1 and 0.5 for a fixed K B = 0.1. With an increase in m * from 0.1 to 0.5, a delay in the transition from the inverted LCO to the deformed flapping can be seen in figure 19. where C d denotes the drag coefficient and 1.189ρ f L represents the quasistatic fluid mass coefficient for a two-dimensional plate per unit width (Lighthill 1960;Newman 1977). The quasistatic equilibrium foil profile can be predicted by solving ( C d = 2.12. This demonstrates the existence of a critical K B value of 0.125 above which the flexible foil can no longer be supported by the fluid force. The quasistatic equilibrium between the restoring structural and fluid forces plays a role in the transition from the inverted-LCO to the deformed-flapping regime.

Formation of a leading edge vortex
The static divergence of an inverted flexible foil explains why the inverted foil is more prone to the flapping instability. However, this does not fully explain the cause of the large-amplitude oscillations observed in the case of the inverted limit-cycle flapping regime presented in § 5.2. The peak-to-peak amplitudes observed for this regime as a function of K B , for fixed Re = 1000 and m * = 0.1, are approximately six times greater than the peak-to-peak amplitudes observed for the conventional foil counterpart. Since vortex structures play a key role in determining the periodic loading on a flexible body, we analyse the wake topology for K B = 0.4, as shown in figure 22. A leading edge vortex (LEV) slowly develops behind the flexible foil during the upstroke. This leading edge vortex remains attached throughout the stroke and grows continuously until the leading edge reaches the maximum displacement prior to the stroke reversal. At the stroke reversal, the leading edge vortex pairs up with a counter-rotating vortex from the lower surface of the foil at the trailing edge and is shed into the wake. Meanwhile, a rotational starting vortex (RSV) is being formed at the leading edge, as shown in figure 22(c,d). The fact that the leading edge vortex produces a large low-pressure region behind the flexible foil results in large-amplitude oscillations. Similar leading edge vortices have been observed in the cases of insect flight (Wang 2000) and delta wings (Ol & Gharib 2003). These leading edge vortices are generally attributed to the effective aerodynamic behaviour observed in the case of insect flight. foil generates a regular von Kármán vortex street with two counter-rotating vortices. However, the flapping inverted foil exhibits a wide range of vortex patterns depending on the non-dimensional bending rigidity K B . Earlier, Kim et al. (2013) observed that a pair of leading edge vortices are shed per one oscillation. In our numerical study, the number of leading edge vortex pairs can be greater than one during the flapping oscillation, as illustrated in figures 23 and 24. For the fixed-point stable regime i.e. for K B 0.545, the flow field typically represents the boundary layer flow over a flat plate. As the flow travels, a symmetric boundary layer flow can be observed on either side of the inverted foil, off the fixed trailing edge, to form a narrow and steady wake. This wake is very similar to the one observed for the case of the fixed-point regime of conventional flexible foil. By decreasing K B < 0.545, the inverted flexible foil experiences large-amplitude oscillations with very low frequency. Figure 23 displays a sequence of plots showing the vorticity contours in the fluid for the case of K B = 0.4 over one downstroke cycle corresponding to time tU 0 /L ∈ [24.7, 29.7]. In this figure, the solid and dashed lines represent the positively and negatively signed vortices. Figure 23(a) displays the instantaneous vorticity contour at tU 0 /L = 24.7, as the foil is about to cross the line Y/L = 0. In this figure, the leading edge vortex (B) that is shed during the previous upstroke can be seen forming a vortex pair with the oppositely signed vortex (C) shed from the lower surface at the trailing edge. The rotation starting vortex (D) formed during the upstroke reversal is still attached to the foil as it crosses the centreline. In figure 23(b), the vortex pair (B)+(C) is shed into the wake; meanwhile, the rotation   Kim et al. (2013), in this case we observe that a pair of leading edge vortices (D) and (H) are shed per half-cycle.
As explained earlier, in the case of the flipped-flapping regime the leading edge of the flexible foil curves around the trailing edge and performs a flapping motion behind the fixed trailing edge. In order to avoid any confusion, we refer to the trailing and leading edges as the fixed and flipped edges respectively. With regard to the wake patterns, the flipped-flapping regime is very similar to the limit-cycle regime of a conventional flexible foil. Figure 25 shows the time history of the steady von Kármán wake of the inverted flexible foil for the flipped-flapping regime with a single distinct frequency at K B = 10 −3 . This figure shows the instantaneous vorticity at six different points over one periodic oscillation. Figure 25 Figure 26 illustrates the time development of the wake topology consisting of a 2S + 2S vortex pattern for the case of the flipped-flapping regime with two distinct frequencies at K B = 0.05. In contrast to the flipped-flapping regime with a single frequency, the wake topology does not represent the von Kármán vortex street. However, a periodic vortex shedding repeats after two complete cycles. The vortex wake in figure 26(f ) corresponding to tU 0 /L = 23.38 is very similar to the wake at tU 0 /L = 18.98, and four vortices (B)-(E) are shed over two cycles. In the flipped-flapping regime, we observe complicated flow features, namely a rolled-up shear layer over the semicircular configuration, vortex-foil interaction and a reversed flow inside the deformed cavity-like region, which are not shown here in detail.

Net energy transfer
In order to analyse the ability of an inverted flexible foil to extract energy from the surrounding fluid flow, we provide a comparison between the energy harvesting estimates for the conventional and inverted flexible foil configurations for low m * = 0.1.  Figure 27 shows the evolution of the non-dimensional bending strain energy (E s ) and the net work done (W s ) on a conventional foil for K B = 3.27 × 10 −4 , Re = 1000 and m * = 0.1. Here, E s and W s are defined as is the curvature of the deformed flexible foil and f (x) is a piecewise polynomial function of sixth order that has been constructed to define the deformed foil profile  figure 27(c-f ) correspond to the points (i)-(iv) indicated in figure 27(a). The strain energy E s of the flexible foil attains a local minimum for the maximum deformation of the trailing edge. In contrast, the flexible foil reaches the maximum E s as the trailing edge is about to cross the centreline. The strain energy E s increases from a local minimum to a local maximum due to the work done by the fluid force on the flexible foil. The flexible foil eventually loses this attained strain energy E s by performing work on the surrounding fluid. Interestingly, for the conventional foil flapping, we observe that the work done by the fluid in developing the structural strain energy is greater than the work done by the structure on the fluid. Therefore, the work done by the fluid per oscillation is greater than zero i.e. W s > 0, and the nett work done by the fluid continuously increases with time at a constant rate, as shown in figure 27(b). This observation reveals that a conventional foil configuration behaves like a non-conservative system. The reason for W s > 0 per oscillation cycle can be attributed to the non-conservative (dissipative) nature of the fluid force and the path followed by the conventional flexible foil. From figures 13 and 27(c-f ), it can be realized that the conventional flexible foil does not follow an identical path during its upstroke and downstroke. Since the work done by a non-conservative force is path-dependent, the work done by the fluid is different from the work performed by the structure. Figure 28 summarizes the maximum and r.m.s. non-dimensional strain energies as a function of K B ∈ [10 −5 , 10 −3 ] for Re = 1000 and m * = 0.1. For sufficiently large K B corresponding to the steady stable mode, the maximum and r.m.s. bending strain energies are close to zero. However, for K B corresponding to the unsteady flapping regime i.e. K B < 7.0 × 10 −4 , the maximum and r.m.s. strain energy values are of O(10 −4 ). The sudden drop in the strain energy for K B < 3.27 × 10 −4 can be attributed to the transition from the fundamental mode k x = 2π to the higher 3π flapping mode. By comparing local maximum strain energies for K B ∈ [7.5 × 10 −5 , 3 × 10 −4 ] and [3 × 10 −4 , 7 × 10 −4 ], we can conclude that lower flapping modes can extract greater bending strain energy compared with the higher bending mode.
In order to determine the effectiveness of the flapping phenomenon in extracting energy, we define the ratio of the maximum strain energy to the total available fluid kinetic energy as where ( E s ) max is the maximum change in the strain energy per oscillation, δ max y represents the maximum tip displacement and T is the time taken to complete one half-cycle. As a typical case we consider K B = 3.27 × 10 −4 , which exhibits maximum non-dimensional bending strain energy, to examine the effectiveness of conventional flexible foils at low m * . On substituting δ max y = 0.12, E max s = 2.3 × 10 −4 , ρ f = 1000 and T = 0.4 observed from the DNS results for K B = 3.27 × 10 −4 , we obtain R = O(10 −7 ), which is negligibly small. Therefore, we show that even though conventional foil can be used to harvest the energy from surrounding fluid flows, the effectiveness of a conventional flexible foil for low m * < 1 in extracting energy from the surrounding fluid flow is very low. A similar observation has been made by Michelin & Doare (2013) for m * < 1. Figure 29(a,b) shows the evolution of the non-dimensional bending strain energy (E s ) and the work done by the fluid force (W s ) on the inverted flexible foil for tU 0 /L ∈ [0, 40] for K B = 0.4, Re = 1000 and m * = 0.1 and 8. Interestingly, we observe that W s closely follows E s in both of the cases. However, the contribution of the structural kinetic energy to the bending strain energy increases with m * . The magnitude of the strain energy developed in the inverted foil is typically O(10 3 ) times the value observed for the conventional foil. In contrast to the E s of the conventional foil, the E s of the inverted foil attains a local maximum for the maximum deformation of the leading edge and a local minimum as the leading edge crosses the centreline. Since the inverted-LCO regime of the inverted foil flapping follows an identical path for both the upstroke and the downstroke (figure 14c), the work performed by the fluid will be the same as the work done by the structure. Therefore, we can observe from figure 29(a,b) that W s ≈ 0 per cycle for tU 0 /L 30 when m * = 0.1 and for tU 0 /L 48.6 when m * = 8. The inverted-LCO flapping regime represents a conservative system, i.e. W s ≈ 0 per cycle, and is more effective in extracting the energy from the surrounding fluid flow in comparison to the conventional foil. Figure 30 summarizes the maximum strain energy developed in a flexible foil as a function of K B ∈ [10 −3 , 1] for Re = 1000 and m * = 0.1. The experimental data (Kim et al. 2013) corresponding to a higher mass ratio of O(1) and Reynolds numbers of O(10 4 ) have also been shown for a qualitative comparison. It can be seen that both the experimental data and the current numerical simulation show a constant maximum strain energy for the inverted-LCO regime. However, the qualitative comparison shows that the current two-dimensional numerical simulations overestimate the maximum strain energy values. We attribute this difference to three-dimensional effects which have been neglected in the current numerical study. The aspect ratios of the flexible foils used in Kim et al. (2013) were in the range of 1-1.3. In addition, we attribute the discrepancy in the maximum strain energy to the large difference in the mass ratio and the Reynolds numbers of the experimental study and the current numerical investigation. Another observation from figure 30 is that the maximum strain energy developed in the flipped-flapping regime is relatively smaller when compared with the inverted-LCO regime. However, the maximum strain energy developed in the flipped-flapping regime is still approximately O(10 3 ) times the maximum value developed in the conventional foil. Figure 31 shows the ratio R of the maximum bending strain energy to the total kinetic energy of the incoming fluid flow, as given in (5.12). It should be noted that (5.12) considers the maximum change in the E s per oscillation, i.e. ( E s ) max , instead of the total stored E s . Use of the total E s can include a large non-oscillatory part; however, that would not be useful for energy recovery. The effects of this consideration can be explained for the case of K B = 0.1, corresponding to the deformed-flapping regime, where there is a large stored strain energy yet the total energy available for energy harvesting is much less than that of K B = 0.2 and 0.3. A comparison with the experimental data of Kim et al. (2013) has also been presented in the figure. The energy efficiency ratio R shows good agreement with the experimental data for K B ∈ [0.1, 0.2]. We also observe a significant difference in the values for K B ∈ [0.2, 0.3]. This difference in the value of R is mainly due to the difference in the maximum bending strain energy. Based on E max s and R values, we confirm that the inverted-LCO regime is the most relevant one for energy harvesting applications among the three flapping regimes of the inverted foil.

Concluding remarks
Through a series of direct numerical simulations, we have studied the flapping instability of the inverted flexible foil configuration in an unbounded axial flow. In this configuration, the trailing edge of the foil is clamped and the leading edge is free to oscillate. We investigated the evolution of the flapping instability as a function of the critical non-dimensional parameters, the non-dimensional bending rigidity K B , the Reynolds number Re and the mass ratio m * , and identified three distinct stability regimes, namely (i) fixed-point stable, (ii) deformed steady and (iii) unsteady flapping state. We analysed the underlying physical mechanism of the onset of flapping instability from the fixed point to the deformed state and the transition from the deformed to the unsteady periodic state. With the aid of DNS results and an analytical model, we have shown that unlike a conventional foil an inverted flexible foil experiences a static-divergence instability much before the onset of flapping. Under the influence of the static-divergence instability, the inverted flexible foil eventually undergoes unsteady flapping stability due to the flow separation at the leading edge. We further examined the effects of Re and m * on the stability regimes of the inverted foil. It was observed that the mass ratio has a small influence on the stability response regimes for a range of m * ∈ [0.1, 8].
Three different unsteady flapping response regimes were qualitatively observed as a function of decreasing K B : (i) inverted limit-cycle oscillation, (ii) deformed-flapping and (iii) flipped-flapping regimes. In contrast to a conventional foil, the inverted flexible foil loses its stability through static divergence for K B 0.55 for Re = 1000 and m * = 0.1. We observed that the static deformation gradually increases with decreasing K B = 0.55 to K B = 0.545. Through a series of numerical simulations for m * = 0.1, we proposed a relation between K B and Re to predict the flapping instability in the inverted flexible foil. With the help of a theoretical model, we showed that the transition from the inverted-LCO to the deformed-flapping regime is due to the existence of a quasistatic equilibrium between the fluid force and the structural restoring force. For mass ratio m * ∈ [0.1, 10], the flapping frequency of the inverted LCO is strongly influenced by the increase in m * , while the cross-stream flapping amplitudes are weakly affected. For the flipped-flapping regime, an increase in m * resulted in a transition from periodic to non-periodic flapping with variable amplitudes and frequencies. The three response regimes can have a wide variety of wakes consisting of up to 14 vortices per oscillation cycle. The inverted limit-cycle flapping regime generates novel 6P + 2S (14 vortices) and 4P + 2S (10 vortices) wakes. The flipped-flapping response regime, on the other hand, is characterized by either a regular von Kármán vortex street or a periodic wake consisting of four vortices.
We have shown that the flapping of inverted flexible foil can generate O(10 3 ) times more strain energy in comparison to conventional flexible foil flapping, which has a profound impact on energy harvesting devices. Although the current two-dimensional simulations were performed at significantly lower Reynolds number and mass ratio, the bending strain energy developed in the inverted foil showed good qualitative agreement with the experimental data presented in Kim et al. (2013). The difference in the magnitude of the strain energy can be attributed to the three-dimensional effects and also to differences in the values simulated. Out of the three response regimes observed, in the context of exploring an alternative energy harvesting mechanism we confirmed the inverted-LCO regime as the most useful regime.