Flutter of a rectangular plate

Flutter of a rectangular plate ⋆ Christophe Eloy ∗ , Claire Souilliez, Lionel Schouveiler IRPHE, CNRS, Universités Aix-Marseille I & II, 49 rue Joliot Curie, 13013 Marseille, France Abstract We address theoretically the linear stability of a variable aspect ratio, rectangular plate in an uniform and incompressible axial flow. The flutter modes are assumed to be two-dimensional but the potential flow is calculated in three dimensions. For different values of aspect ratio, two boundary conditions are studied: a clampedfree plate and a pinned-free plate. We assume that the fluid viscosity and the plate viscoelastic damping are negligible. In this limit, the flutter instability arises from a competition between the destabilising fluid pressure and the stabilising flexural rigidity of the plate. Using a Galerkin method and Fourier transforms, we are able to predict the flutter modes, their frequencies and growth rates. The critical flow velocity is calculated as a function of the mass ratio and the aspect ratio of the plate. A new result is demonstrated: a plate of finite span is more stable than a plate of infinite span. Key words: flow-induced vibration, cantilevered flexible plate, flutter instability, variable aspect-ratio PACS: 46.40.-f, 46.40.Ff, 46.70.De 1 Introduction In his seminal paper on the instability of jets, Lord Rayleigh (1879) suggested that his theoretical approach could be used to prove that an infinite flag is always unstable. Indeed, it can be proved easily that an elastic plate of infinite dimension (both in the span- and streamwise direction) is always unstable ⋆ An earlier version of this paper was presented in the 7th FSI, AE & FIV+N Symposium, held within the 2006 PV & P Conference in Vancouver, BC, Canada. ∗ Corresponding author. Email address: christophe.eloy@irphe.univ-mrs.fr (Christophe Eloy). Preprint submitted to Journal of Fluid and Structure 18 December 2006 when immersed in an axial potential flow. However, this fluid-structure interaction problem becomes far more complex mathematically when the finite dimensions of the flag are explicitly taken into account. This is the subject of the present paper. Using analytical tools of airfoil theory, Kornecki et al. (1976) have shown that a plate of infinite span but finite chord was stable for flow velocities below a critical velocity (in the following, we will use the term plate instead of flag to emphasise the importance of the finite density and finite bending stiffness of the material). Kornecki et al. (1976) assumed an elastic plate and used two different theoretical approaches to model the flow around this plate. They first assumed a potential flow with zero circulation. Then, using a method introduced by Theodorsen (1935), they added a distribution of vorticity in the plate wake to smooth out the trailing edge singularity of the pressure field. This results in an unsteady circulatory flow. More recently these approaches have been used again with better computer accuracy (Huang, 1995; Watanabe et al., 2002a). Another theoretical approach has been used by Guo and Paı̈doussis (2000). They solved the two-dimensional problem (assuming infinite span) in the Fourier space for a potential flow. The present paper can be viewed as an extension of the analysis of Guo and Paı̈doussis (2000) to take into account the finite span of the plate. An important aspect of these theoretical models is the way they deal with the flow at the boundary conditions. Kornecki et al. (1976) in their zero-circulation model had pressure singularities both at the trailing and the leading edge. The use of the Theodorsen (1935) theory in the second model of Kornecki et al. (1976) suppresses the trailing edge singularity by the use of the Kutta condition. In this unsteady problem, even if the mathematical proof has yet to be given (Frederiks et al., 1986), the Kutta condition can indeed be used because the Reynolds number is usually very large and the flutter frequency is of order one (Crighton, 1985). Finally, the theoretical model of Guo and Paı̈doussis (2000) solves the pressure distribution problem in the Fourier space assuming implicitly no singularities. This means that an incoming “wake” has to be added to the flow to suppress the leading edge singularity. This “wake” cannot be justified on physical grounds. Surprisingly, these three different twodimensional models give almost the same results for the critical velocity of the instability (see Watanabe et al., 2002a). This means that the wakes added upstream or downstream do not have a great influence on the stability of this fluid-structure interaction. Shayo (1980) first attempted a three-dimensional stability analysis to understand the dependence of the critical velocity on the plate span. In his study, he made several mathematical assumptions to simplify the calculation that led him to conclude that a flag of infinite span is more stable than a finite one. This latter result is in contradiction with a slender body approach (Lighthill, 2 1960; Datta and Gottenberg, 1975; Lemaitre et al., 2005). This discrepancy was reexamined by Lucey and Carpenter (1993) who addressed theoretically the linear stability of a compliant panel of finite aspect ratio embedded in a rigid wall (this case bears strong similarity with the flutter of a flag in axial flow). They found that a compliant wall of finite span is always more stable than its infinite counterpart (see also the review of Yadykin et al., 2003, on the added mass of a flexible plate of finite aspect ratio). This seems in contradiction with the results of Shayo (1980). More recently, Argentina and Mahadevan (2005) studied the flutter instability of a flag with a simple twodimensional model based on the theory of Kornecki et al. (1976). Using the data of a three-dimensional numerical simulation, their model was then adjusted to take into account the finite the aspect ratio of the flag. They showed qualitatively that the finite span tends to stabilise the system. As we shall see in the present paper, a theoretical analysis can be carried out without any assumptions and it shows that the smaller the plate span, the more stable the system, in agreement with the results of Lucey and Carpenter (1993) and Argentina and Mahadevan (2005). Experiments with plates made of metal, paper or plastic sheets have been carried out by Taneda (1968), Datta and Gottenberg (1975), Kornecki et al. (1976) and more recently by Yamaguchi et al. (2000); Watanabe et al. (2002b); Shelley et al. (2005); Souilliez et al. (2006) (see the book of Paı̈doussis, 2004, for a comprehensive list of references). These experiments showed that the flutter modes observed at threshold are always two-dimensional. They also showed that the instability threshold is always larger than the theoretical predictions. The work of Watanabe et al. (2002a) shows that the critical velocity measured in the experiments is at least twice as large as the analytical and numerical predictions for all experimental parameters. So far, no satisfactory explanation of this apparent discrepancy has been given. Cantilevered plates in axial flow have also been modelled numerically by Watanabe et al. (2002a), Balint and Lucey (2005) and Tang and Paı̈doussis (2006). In these studies, a two-dimensional solver based on the Navier-Stokes equations or on a vortex method has been combined to a linear beam model for the plate. The critical velocities obtained with these numerical simulations are similar to the results of Kornecki et al. (1976) and Guo and Paı̈doussis (2000). In their papers, Guo and Paı̈doussis (2000) and Balint and Lucey (2005) also provide a mechanism of irreversible energy transfer from the fluid to the structure that gives another insight into the the instability mechanism. In this paper, we present an analytical stability analysis of an elastic plate immersed in an axial uniform flow. The analysis assumes a two-dimensional flutter motion (as it has been observed in the experiments so far), but takes into account explicitly the finite span and chord of the plate to calculate the surrounding flow. This paper is divided into four sections. We will first present 3 the principle of the stability analysis involving a Galerkin method. The threedimensional flow will then be calculated in the Fourier space. The resulting eigenmodes and critical velocities and their dependence on the aspect ratio will then be discussed. Finally, we will see how these results may explain the observed critical velocities in the experiments. For the sake of clarity, the most technical parts of our analytical developments have been postponed to Appendices at the end of the paper. 2 Principles of the stability analysis In this section, the different steps of the stability analysis are detailed. Starting from the equation of motion for the plate and expanding its two-dimensional deflection on Galerkin modes, the initial problem is transformed into an eigenvalue problem for the flutter modes and their frequencies. 2.1 Equation of motion As shown on Fig. 1, the plate lies in the vertical OXZ plane and the wind is parallel to OX. Assuming a two-dimensional deformation of the plate, the lateral deflection Y is a function of the streamwise direction X and the time T only. The equation of motion of the elastic plate is given by ρp ∂2Y ∂4Y + D = h∆P iZ , ∂T 2 ∂X 4 (1) where ρp is mass per unit surface of the plate, D its flexural rigidity (given by D = Eh3p /[12(1 − ν 2 )], E being the Young’s modulus, hp the plate thickness and ν the Poisson’s ratio), ∆P is the pressure difference between the two sides of the plate and h.iZ denotes the average along the spanwise direction Z H — 2 Y X L U Y(X,T) 0 H -— 2 Fig. 1. Schematic of the plate subject to a two-dimensional deflection Y (X, T ) in an axial flow of velocity U . L is the plate chord and H its span. 4 for −H/2 < Z < H/2. Because of the antisymmetry of the flow with respect to the plate plane OXZ, perturbation pressures on upper and lower sides are equal and opposite. Therefore ∆P = −2P + , where P + is the perturbation pressure on the side Y > 0. In Eq. (1) above, we have assumed that the viscoelastic damping of the material and the tension due to the viscous boundary layers are negligible. If needed, the form of the present analysis facilitates the introduction of these terms in future studies. The plate chord L and the flow velocity U are used as characteristic length and velocity respectively. Using lowerscript letters for dimensionless quantities, we have Y Z TU X x= , y= , z= , t= . (2) L L L L The above equation of motion (1) is made dimensionless ∂2y 1 ∂4y + = M ∗ hpiz , ∂t2 U ∗ 2 ∂x4 (3) where p is the dimensionless pressure difference, U ∗ is the reduced flow velocity and M ∗ the mass ratio given by 2P + , ρa U 2 r ρp ∗ LU, U = D ρa L M∗ = , ρp p=− (4) (5) (6) with ρa the density of the surrounding fluid (usually air). For an inviscid surrounding fluid and a perfectly elastic material with no damping, the motion of the plate is entirely governed by three dimensionless parameters: U ∗ , M ∗ and the aspect ratio H ∗ which is equated as H∗ = H . L (7) 2.2 Galerkin decomposition Equation of motion (3) becomes linear if p is linear with respect to y (this is the case if the flow is supposed to be potential). In this case, we can assume a Galerkin decomposition of the lateral deflection and switch to the frequency 5 domain such that the deflection y can be written as y(x, t) = N X an yn (x) eiωt +c.c., (8) n=1 where c.c. denotes the complex conjugate, N is the truncation of the Galerkin expansion and yn are the beam eigenfunctions in vacuo satisfying the free boundary-condition at the trailing edge and either a clamped or pinned boundarycondition at the leading edge (see Paı̈doussis, 1998, 2004). The dimensionless angular frequency ω is such that ω = ΩL/U (where Ω is the dimensioned angular frequency). Note that these Galerkin modes satisfy the following relation d4 yn = kn4 yn , 4 dx (9) where kn are wavenumbers sorted in ascending order. The dimensionless wavenumbers kn are related to the dimensioned ones Kn through kn = Kn L. To solve the equation of motion (3) using the Galerkin expansion (8), one needs to determine the pressure pn (x, z) associated with the deflection yn (x) of the plate at angular frequency ω. This will be done in the next section. Upon defining the standard scalar product as f ⊗g = Z 0 1 f (x)g(x) dx, (10) n the Galerkin modes satisfy the orthogonality condition ym ⊗yn = δm , where δij is the Kronecker delta. By multiplying Eq. (3) by the modes ym , the equation of motion becomes a system of N linear equations with N unknowns (the amplitudes an ). The solvability condition imposes   1 det −ω I + ∗ 2 Q − M ∗ P = 0, U 2 (11) where I is the N × N identity matrix, Q is the diagonal matrix with the kn4 in ascending order on its diagonal and Pmn (ω) = ym ⊗ hpn (x, z)iz . (12) Note that the matrix P is a function of ω as it will be shown below. For a given mass ratio M ∗ and dimensionless velocity U ∗ , the solvability condition (11) has 2N solutions ωi corresponding to each flutter mode. For an infinite Galerkin truncation N , there is an infinite, countable number of pairs of conjugate numbers ωi solution of Eq. (11). The frequencies of the flutter modes are simply given by ℜ(ωi ) and their growth rates by σi = −ℑ(ωi ). To determine the stability of the plate flutter, one now needs to calculate the 6 matrix P and search the velocity U ∗ needed to have at least one flutter mode with a positive growth rate. 3 Three-dimensional flow In this section, we calculate the potential flow above the plate in the limit of large aspect ratio H ∗ . The calculation is done in the Fourier space. As we shall see, we will need to solve two singular integral equations to obtain the pressure on the plate, one along the span and the other along the chord. The presence of these singularities stress the intrinsic singular aspect of the flow in this problem. 3.1 Perturbation potential in the Fourier space Assuming an inviscid and irrotational flow, the pressure pn is given by the unsteady Bernoulli equation ∂ϕn + iωϕn pn = 2 ∂x ! , (13) y=0 where ϕn (x, y, z) is the perturbation velocity potential at frequency ω in the half-plane y > 0. The factor 2 comes from an opposite perturbation pressure on the other side of the plate. The potential ϕn satisfies a Laplace equation ∂ 2 ϕn ∂ 2 ϕn ∂ 2 ϕn + + = 0, ∂x2 ∂y 2 ∂z 2 (14) and the following boundary conditions ϕn |x→±∞ = ϕn |y→∞ = ϕn |z→±∞ = 0, ∂ϕn ∂y = iωyn + y=0 ! (15) dyn , for (x, z) ∈ D, dx (16) ∂ϕn = 0, for (x, z) ∈ / D, pn = 2 + iωϕn ∂x y=0 where the domain D corresponds to the plate   H∗ H∗ D = (x, z) | x ∈ [0 1] and z ∈ − 2 2  (17) . (18) The second boundary condition (16) derives from the continuity of the normal velocity above the plate. The third boundary condition (17) ensures the 7 continuity of the pressure outside the plate in the Oxz plane. Note that this boundary condition allows a solution of the form ϕn (x, y = 0) = ϕ0 exp(−iωx). This is equivalent to a vorticity distribution in the Oxz plane written as a progressive wave travelling at the velocity U (coming back to dimensioned units). Downstream of the trailing edge, this wave models the “wake” of the plate; upstream of the leading edge, it has no physical meaning. The difficulty of the mathematical problem arises from the heterogeneity of the boundary conditions. If Eq. (17) was a condition on the y-derivative of ϕn as Eq. (16), we would have a Neumann problem, which would be simpler to solve. Here, the boundary conditions (16) and (17) involve different derivatives of ϕn leading to singularities along the edges of the plate. To solve the set of equations (14–17), we switch to the Fourier space both in x- and z-direction. Using the Laplace equation (14), we have ϕn (x, y, z) = Z ∞ −∞ Z ∞ −∞ ϕen (kx , kz ) e−κy eikx x eikz z dkx dkz , where ϕen (kx , kz ) is the Fourier transform of ϕn (x, 0, z) and κ= q kx2 + kz2 . (19) (20) The rest of the present section involves lengthy and technical calculations; readers not interested in these technical details may jump to the results section. 3.2 Boundary conditions in the Fourier space Using Eqs. (13) and (19) gives the pressure above the plate in the Fourier space Z ∞ 2i (kx + ω) pen (kx , z) eikx x dkx , (21) pn (x, z) = −∞ where pen is pen (kx , z) = Z ∞ −∞ ϕen (kx , kz ) eikz z dkz . (22) The function pen is proportional to the pressure along the span for a given wavenumber kx . For an infinite plate in the z-direction (as it has usually been assumed in the previous theoretical studies), pen is constant along the z-axis. Taking the z-derivative of Eq. (22), executing an inverse Fourier transform and using the boundary condition shown in Eq. (17) gives ϕen (kx , kz ) = 1 Z H ∗ /2 ∂ pen −ikz z e dz. 2πikz −H ∗ /2 ∂z 8 (23) Now taking the x-derivative of Eq. (21), its inverse Fourier transform and using again the boundary condition (17) gives Z 1 1 ∂pn −ikx x e pen (kx , z) = − dx. 4πkx (kx + ω) 0 ∂x (24) These two latter equations allow the expression of the boundary condition (17) in the Fourier space. To solve the problem, we now need to use the boundary condition (16) expressing the continuity of normal velocity above the plate. Using the Fourier transform (19) with the boundary condition (16) gives Z ∞ ven (kx ) eikx x dkx = iωyn + −∞ with ven (kx ) = Z ∞ −∞ dyn , for (x, z) ∈ D, dx −κϕen (kx , kz ) eikz z dkz . (25) (26) The right hand side of Eq. (25) [iωyn + dyn / dx] is a function of x only. Therefore ven (kx ) which should be a function of both kx and z does not depend on z for −H ∗ /2 < z < H ∗ /2. Using Eq. (23) with (26) and inverting the integral signs gives Z H ∗ /2 −H ∗ /2 Z ∞ κ ikz (z−v) ∂ pen e dkz dv = −2π ven (kx ), (kx , v) ∂v −∞ ikz (27) where v is a dummy variable in place of z (therefore, this variable v is not related to ven ). Solving this integral equation for ∂ pen /∂v allows the expression of the dependency of the pressure on the z-direction. 3.3 Pressure along the span Upon rescaling z, v and kx with the characteristic length H ∗ /2 in such a way that 2v 2z |kx |H ∗ v = ∗ , z = ∗ , kx = , (28) H H 2 equation (27) becomes Z 1 −1 p′n (v)F (k x (z − v)) dv = π, for z ∈ [−1 1], (29) where F is given by the integral F (X) = Z 0 ∞ s 1+ Z X K1 (α) 1 sin(αX) dα = − dα, α2 α 0 9 (30) where Kν is the modified Bessel function of the second kind and pn depends on pen in the following way pn (v) = − |kx | pen (kx , v), ven (kx ) (31) where the prime denotes derivation with respect to v. The function pn implicitly depends on k x . The solution of the singular integral equation (29) for p′n is outlined in the Appendix A. The main results are presented on Fig. 2. For large k x (which means that the wavelength of the deflection 2π/kx is small compared to the span H ∗ ), the pressure profile is constant except close to the edges of the plate. It converges towards a top-hat function for infinite k x , or equivalently infinite span (this is what was expected from two-dimensional theory). For small k x , the pressure profile along the span is elliptic q pn (v) ∼ k x 1 − v 2 when k x → 0. (32) We are interested in the average of pn along the span which is a function of k x only. In Fig. 3, this average is plotted: the symbols represent the results of a symbolic-numerical calculation and the two plain curves are asymptotic expansions valid in the limit of either small or large k x . For small k x , the average is obtained from Eq. (32) and is simply hpn iz =  2 πk x when k x → 0. + O kx 4 (33) 1.2 1 100 10 0.8 1 0.6 pn 0.4 0.2 0.1 0 −0.2 −1 −0.5 0 0.5 1 v Fig. 2. Repartition of the rescaled pressure pn as a function of the rescaled spanwise direction v. Each curve corresponds to a different value of k x : 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100 (from bottom to top). 10 This result allows the recovery of the slender body approximation of Lighthill (1960) (see also Lemaitre et al., 2005) as detailed in Appendix C. For large k x another expansion can be used: hpn iz = 1 −   1 + o e−kx when k x → ∞. 2k x (34) We have called this expansion “large span approximation” since it is only valid for sufficiently large aspect ratio H ∗ . In Fig. 3, the dashed line shows an approximation valid through the whole domain of k x hpn iz ≈ 1 − h 1 2k x + exp (π/4 − 2)k x i. (35) This approximation gives the correct average pressure within 3% error (this is the maximum difference between the dashed curve and the crosses in Fig. 3). It also has the advantage to give the correct mathematical expansion given by Eqs. (33) and (34) in the limit of small and large k x . Coming back to the former dimensionless variables kx and z, ven can be expressed from Eqs. (31) and (35) ven (kx ) = −hpen (kx , z)iz |kx |α (|kx |H ∗ ) , (36) where h.iz still denotes the average along the span for −H ∗ /2 < z < H ∗ /2 1 large span approximation 0.8 0.6 �pn �z 0.4 0.2 0 −2 10 slender body approximation −1 0 10 10 1 10 2 10 kx Fig. 3. Average of the pressure along the span hpn iz as a function of k x . Crosses are the results of a symbolic-numerical calculation. The plain curves are the slender body and large span approximations, given by Eqs. (33) and (34) respectively. The dashed line is the approximation given by Eq. (35). 11 and α is the following function α (X) ≈ 1 + 1 . X − 1 + exp [(π/8 − 1)X] (37) 3.4 Pressure along the chord Inserting Eqs. (36) and (24) into Eq. (25) and inverting the integral signs gives 1 Z1 ′ dyn hpn (u, z)iz G(x − u, H ∗ ) du = iωyn + , 4π 0 dx (38) where the prime denotes the derivation with respect to the dummy variable u (in place of x) and ∗ G(X, H ) = Z ∞ −∞ |kx | α (|kx |H ∗ ) eikx X dkx , kx (kx + ω) (39) where α is defined in Eq. (37). Equation (38) is a singular integral equation for p′n similar to Eq. (22) in Guo and Paı̈doussis (2000). The only difference is that G depends explicitly on the aspect ratio H ∗ in the present paper. Equation (38) can be solved by multiplying by (iω + d/ dx), leading to dyn d2 yn 1 Z1 ′ hpn (u, z)iz H(x − u, H ∗ ) du = −ω 2 yn + 2iω , + 4π 0 dx dx2 with H(X, H ∗ ) = −2 Z 0 (40) ∞ α (kx H ∗ ) sin (kx X) dkx . (41) This function can be approximated by ! 2 8−π 8 H(X, H ) ≈ − − sgn(X). − ∗ ∗ X H H + 2|X| ∗ (42) This latter expression is a very good approximation (within 2%) of Eq. (41). Using Eqs. (33) and (34) this expression is shown to be the correct exact expansion in the limit of small H ∗ (slender body approximation) and large H ∗ (large span approximation) as shown in Appendices C and D. The average pressure hpn (x, z)iz can be obtained from the singular integral equation (40) on its derivative. Since this integral equation is linear, the pressure can be decomposed into three terms, corresponding to the three terms on the right hand side of Eq. (40) ) (G) (K) hpn (x, z)iz = −ω 2 p(M n (x) + 2iωpn (x) + pn (x), (43) where the superscripts M , G and K correspond to the added mass, the gyroscopicity and stiffness respectively (see Paı̈doussis, 2004). 12 The solution of the integral equation (40) is found by a Gauss-Chebyshev method which is detailed in Appendix B. For a given aspect ratio H ∗ , each pressure term on the right-hand side of Eq. (43) can be evaluated numerically. In the limit of large aspect ratio, an expansion of the pressure in powers of 1/H ∗ is found. We call this mathematically exact expansion “large span approximation” (its derivation is detailed in Appendix D). In the limit of small aspect ratio, the slender body approximation holds. It has been first given by Lighthill (1960) but it can be derived from the present formulation (see Appendix C for details). 4 Results We have now calculated the fluid load on the plate for any given Galerkin mode yn and any given aspect ratio H ∗ . We are thus able to determine the stability of the system by using the solvability condition (11). In this section, the eigenfrequencies and growth rates of the flutter modes are given as a function of the dimensionless flow velocity U ∗ , the mass ratio M ∗ and the aspect ratio H ∗ . 4.1 Eigenmodes and eigenfrequencies For a given aspect ratio H ∗ and a given mass ratio M ∗ , we can now calculate the complex frequencies of the flutter modes ω as a function of the dimensionless velocity U ∗ . The results of such a calculation are given in Fig. 4 for a clamped-free plate when M ∗ = 1 and H ∗ = ∞. The top curve shows the real part of the frequency ℜ(ω) multiplied by U ∗ , and the bottom curve show the growth rates σ = −ℑ(ω) for the first 4 eigenmodes. We chose to plot the product ℜ(ω)U ∗ instead of ℜ(ω) alone because it is simply equal to kn2 in vacuo. This can be seen by taking the limit of vanishing M ∗ in Eqs. (3) and (9). Therefore, one expects the eigenfrequencies to be close to those values if the flow simply excites the natural modes of vibration in vacuo. The first values of kn2 are 3.5, 22, 62, 121 and we observe that indeed, in the limit of vanishing U ∗ , the eigenfrequencies in Fig. 4 are close to those values but slightly below. This shift is simply due to an added mass effect. For higher values of the flow velocity U ∗ , both the eigenfrequencies and the growth rates change. For U ∗ < 5.12, the growth rates are all negative meaning that all modes are stable. For U ∗ > 5.12, the second mode becomes unstable and for higher flow velocity, the third and fourth modes eventually become unstable too as seen in Fig. 4. Slightly above the destabilisation of the third 13 mode (at U ∗ = 11.0), the second mode becomes stable again for U ∗ > 13.5. The critical flow velocity Uc∗ is defined as the lower value of the flow velocity for which a mode becomes unstable. In this example, we have Uc∗ = 5.12. This critical flow velocity is implicitly a function of the mass ratio M ∗ and the aspect ratio H ∗ . We shall see below how these two dimensionless parameters modify the critical velocity. In Fig. 4, the first mode appears to be always stable, as it is the case for all values of parameters M ∗ and H ∗ we have tested. For 6.32 < U ∗ < 14.2, the branch of this first flutter mode splits into two divergence modes (of zero frequencies) labelled 1a and 1b; however, these divergence modes always remain 120 4 100 80 ℜ(ω)U ∗ 60 3 40 2 20 1 1 1a+1b 0 2 4 6 8 10 U 12 14 16 18 20 18 20 ∗ 2 1.5 3 1 0.5 2 σ 4 0 −0.5 1a 1 −1 1 −1.5 −2 1b 2 4 6 8 10 U 12 14 16 ∗ Fig. 4. Frequencies (top) and growth rates (bottom) of the first 4 modes as a function of the dimensionless flow velocity U ∗ . The mass ratio of the clamped-free plate is M ∗ = 1 and its aspect ratio H ∗ is infinite. The labels number the different flutter modes in order of ascending frequency. 14 stable in the configurations tested (Guo and Paı̈doussis, 2000, have shown that these divergence modes can be unstable for different boundary conditions such as a clamped-clamped plate). 4.2 Critical curves We now explore the variation of the critical velocity Uc∗ as a function of the aspect ratio H ∗ and the mass ratio M ∗ . In Fig. 5, for the same mass ratio M ∗ = 1 as in Fig. 4, we have plotted the critical velocity Uc∗ as a function of H ∗ for a clamped-free plate. Crosses show the results of the exact calculation for different aspect ratios H ∗ and the solid lines are the results of the slender body and large span approximations valid in the asymptotic limit of small and large H ∗ respectively. As seen in Fig. 5, the critical velocity is a monotonically decreasing function of H ∗ . For instance, for H ∗ = 1, the present theory gives Uc∗ = 7.36 whereas it is equal to 5.12 for H ∗ = ∞ (for H ∗ = 1, the large span approximation gives Uc∗ = 8.21). It means that for this value of M ∗ , the aspect ratio has a strong influence on the resulting critical velocity and that it should be taken into account to model the experiments. To emphasise the importance of the aspect ratio in determining the critical velocity, we have plotted on Fig. 6 the critical velocity as a function of the 40 35 30 slender body approximation 25 Uc∗ 20 15 large span approximation 10 5 0 −2 10 −1 0 10 H ∗ 10 1 10 Fig. 5. Critical velocity as a function of the aspect ratio H ∗ for M ∗ = 1 and for a clamped-free plate. Crosses are the results of the exact calculation. The curves on the left and right are the results of the slender body and the large span approximation respectively. 15 mass ratio M ∗ for three different values of H ∗ . The bottom curve shows the classical result obtained for plates of infinite span : H ∗ = ∞. In this case, it is the exact same result as the one obtained by Guo and Paı̈doussis (2000). The two top curves show the same result taking into account a finite aspect ratio H ∗ = 1 and H ∗ = 0.2. As it can be seen from the figure, for small values of the mass ratio (M ∗ . 1), the critical velocity increases as aspect ratio decreases. However, for larger mass ratios (M ∗ & 1.5), the curves tend to gather to give a critical velocity around Uc∗ ≈ 10. This is because, for the aspect ratios considered in this figure, the wavelength of the flutter mode is small compared to the span when M ∗ & 1.5 and thus a two-dimensional theory (H ∗ = ∞) give approximatively the accurate result. On Fig. 6, one can also see two different modes of instability corresponding to the two lobes of the curves. For H ∗ = 1, for M ∗ . 1.5, the first unstable mode is the second mode as it was the case in Fig. 4 for M ∗ = 1. For 1.5 . M ∗ . 6, the third mode is the first unstable one. For larger M ∗ , it is the fourth mode and for even larger M ∗ eventually the 5th, 6th, etc. (not shown on the figure). The transition between modes is slightly shifted when the aspect ratio H ∗ is changed. This explains why the critical curve is higher for H ∗ = ∞ than for H ∗ = 1 in the range 1.32 < M ∗ < 1.42. For particular values of M ∗ , for instance for M ∗ = 1.5 and H ∗ = 1, the linear theory predicts that for increasing U ∗ , the flow becomes unstable for the second mode for 9.2 < U ∗ < 11.2, then stable when 11.2 < U ∗ < 11.4, 30 25 20 H ∗ = 0.2 Uc∗ 15 10 5 0 −1 10 H∗ = 1 H∗ = ∞ 0 10 M 1 10 ∗ Fig. 6. Clamped-free plate. Critical velocity as a function of the mass ratio M ∗ for three different values of the aspect ratio H ∗ as labelled. The symbols are experimental results: + and × for H ∗ = 1 with two different masts (Souilliez et al., 2006); filled circles for 0.6 < H ∗ < 1.5 (Huang, 1995); open diamonds and circles for H ∗ ≈ 0.25, AISI 304 and waxed paper respectively (Yamaguchi et al., 2000). 16 then unstable again for the third mode for U ∗ > 11.4. This peculiar behaviour has already been reported by Lemaitre et al. (2005) in the case of plates of very small aspect ratio. It has never been observed experimentally for two reasons. Firstly, because the instability bifurcation is subcritical meaning that as soon as the instability starts, the amplitude of the mode jumps to another branch not predicted by the present linear stability analysis; secondly, a small material damping added in the theory would smooth out the cusps observed in the curves of Fig. 6 and thus suppress this peculiar property. As Lemaitre et al. (2005) have shown, the addition of a material damping of Kelvin-Voigt type always smooth the critical curves but it may either decrease or increase the value of the critical velocity depending on the values of M ∗ and H ∗ (this is also discussed in Paı̈doussis, 1998). Note that, in the present paper, the material damping could be easily included since it is a linear effect. As shown in Fig. 6, the present theoretical stability curve for aspect ratio H ∗ = 1 compares well with experimental for M ∗ . 1. However, for larger mass ratios M ∗ & 2, the agreement is not as good (for the other materials considered in Yamaguchi et al., 2000, the critical velocity is even higher; typically they find 15 < Uc∗ < 45). We suggest that it is due to three-dimensional effects in the deflection of the plate prior the flutter instability (such as it is the case for a real fabric flag). These 3D effects tend to rigidify the plate and thus delay the instability. Figure 7 shows the critical curves for a different set of boundary conditions: a pinned-free plate. The results bear similarity with those of a clamped-free plate. We recover that the smaller the span, the more stable the system. The same modes are unstable for small M ∗ (first the second mode, then the third). 16 14 12 H∗ = 1 10 Uc∗ 8 6 4 2 0 −1 10 H∗ = ∞ 0 1 10 10 2 10 M∗ Fig. 7. Same as Fig. 6 for a pinned-free plate. As a reference, we plot the results for a clamped-free plate and H ∗ = 1 as a dashed curve. 17 However, there are three main differences: the velocity threshold is slightly higher for the second mode whereas it is lower for the third mode; the transition between the second and third mode is shifted to M ∗ ≈ 0.8; and the critical velocity of the third mode and fourth mode (for M ∗ & 1 and M ∗ & 5) monotonously decreases as M ∗ increases. Again, the peculiar behaviour and the cusps observed for M ∗ & 5 would probably be suppressed if one added a material damping to the present model. Note also that the present theory could also be used to calculate the stability of any set of boundary conditions. 5 Conclusions In this paper, we have developed a new theoretical model that enables the modelling of the flutter of a rectangular cantilevered plate in axial flow. We have assumed that the plate deflection is two-dimensional and immersed in a three-dimensional potential inviscid flow. Starting with the equation of motion of the plate, the stability analysis has been performed using a Galerkin method. The inviscid fluid forces on the plate have been calculated in the Fourier space assuming a finite plate span. The eigenmodes and the corresponding eigenfrequencies and growth rates have been calculated as a function of the three dimensionless parameters of the problem: U ∗ the flow velocity, H ∗ the plate aspect ratio and M ∗ its mass ratio. This theory has also been expanded in the asymptotic limit of large span, giving rise to what we called the “large span approximation”. We have shown that for a given mass ratio, the critical velocity is almost always a decreasing function of the aspect ratio. It means that for a given plate chord, the smaller the span, the more stable the system, as predicted by slender body theory for smaller values of H ∗ . This behaviour could account for the apparent discrepancy between the theoretical models and the experimental results published so far for small M ∗ . Indeed, in the literature, the experiments have always appeared to be more stable than the predictions (Watanabe et al., 2002a). It has been shown in this paper that when aspect ratio is explicitly taken into account, the predictions of the critical velocity are much improved. A detailed experimental study on the flutter of Mylar plates is the subject of a forthcoming paper. The present work could be extended to different configurations. If Eq. (3) that models the plate motion remains linear, many extra effects can be added into the analysis, for instance viscoelastic material damping or tensile effects could be introduced. The boundary conditions on the plate can also be modified by changing the basis of Galerkin modes as we have shown by studying the clamped-free and the pinned-free plate. 18 Acknowledgements We would like to thank Emmanuel de Langre and Olivier Doaré for their encouragements and advice. A Integral equation for the pressure along the span We now want to solve the singular integral equation (29) Z 1 −1 p′n (v)F (k x (z − v)) dv = π, for z ∈ [−1 1], (A.1) for p′n with F (X) still given by Eq. (30). To achieve this purpose, p′n is expanded as a series of Chebyshev polynomials p′n (v) = J X T2j−1 (v) Aj √ , 1 − v2 j=1 (A.2) where Ti are the Chebyshev polynomials of the first kind, Aj the unknown amplitudes and J is the truncation of this expansion [the function F (X) being even, p′n has to be odd, this is why only the odd Ti are needed in the previous expansion]. Introducing the scalar product q 2Z1 f (z)g(z) 1 − z 2 dz, f ⊙g = π −1 (A.3) inserting Eq. (A.2) into Eq. (A.1) and multiplying by the Chebyshev polynomials of the second kind U2q−2 gives the following linear equation BA = Π, (A.4) where A is the unknown column vector with the amplitudes An , Π is the column vector with π on the first row and zeros elsewhere and B is the matrix given by Z 1 T2j−1 (v) √ F (k x (z − v)) dv. (A.5) Bq,j = U2q−2 ⊙ −1 1 − v2 The coefficients Bq,j can be calculated using a symbolic calculation software such as Mathematica. Calculating a 15×15 matrix B for a given k x takes about three hours on a laptop computer. Once these coefficients Bq,j are known, the solution of Eq. (A.4) for A can be evaluated numerically. The Aj solutions can then be inserted into Eq. (A.2), which gives the function pn after integration. This function is plotted in Fig. 2 for different values of k x . 19 B Integral equation for the pressure along the chord We want to solve the singular integral equation (40) for the pressure. To do so, the variables x and u are rescaled such that xb = 2x − 1, ub = 2u − 1. (B.1) The new pressure is Pn (ub) = hpn (u, z)iz , such that Pn′ (ub) = hp′n (u, z)/2iz . The integral equation (40) becomes 1 Z1 ′ c x b−u b , H ∗ ) du b = F (x b), P (ub)H( 4π −1 n (B.2) c x b−u b, H ∗ ) = H(x − u, H ∗ ) and F (x b) is one of the derivatives of yn (x). with H( c The function H can be decomposed into three terms where c X, c H ∗) = H c + H( 1 1 c c (H ∗ ), H2 + H 3 H∗ c (X) c =− 4 , H 1 c X c (X) c = −8 sgn(X), c H 2 c (X, c H ∗ ) = 8 − π sgn(X), c H 3 c ∗ H + |X| (B.3) (B.4) (B.5) (B.6) c /H ∗ dominates and gives rise to the slender For small aspect ratio, the term H 2 body approximation as it is shown in Appendix C. For large aspect ratio, c gives the first order of approximation (independent of H ∗ ) and the term H 1 corresponds to the two-dimensional theory of Guo and Paı̈doussis (2000). To solve the integral equation (B.2) the function Pn′ is first expanded as Pn′ (ub) = M X Ti (ub) Abi √ , 1 − ub2 i=1 (B.7) where Ti are the Chebyshev polynomials of the first kind. Then, injecting the expansion (B.7) into (B.2) and left multiplying by the Chebyshev polynomial of the second kind Uq−1 (using the scalar product ⊙ defined by Eq. (A.3)) gives the linear system M X i=1 Ci,q (H ∗ )Abi = Fq , 20 (B.8) where Fq = Uq−1 ⊙ F and 1 Z 1 Ti (ub) c √ H(xb − ub, H ∗ ) dub. Ci,q (H ) = Uq−1 ⊙ 2 4π −1 1 − ub ∗ (B.9) The matrix C is decomposed into three terms corresponding to the three terms c defined in Eq. (B.3) of H C(H ∗ ) = C (1) + 1 (2) C + C (3) (H ∗ ), ∗ H (B.10) where C (1) is the identity matrix, (2)  √ 4 1 − xb2 Ui−1 (xb) , Uq−1 ⊙ iπ 16 q =− 2 4 , for q, i of same parity, 2 2 π q − 2q (i + 1) + (i2 − 1)2 = 0, for q, i of different parity, Ci,q = (B.11) (B.12) (B.13) and C (3) (H ∗ ) is evaluated numerically for a given aspect ratio. The linear system (B.8) can then be solved to express the amplitudes Abq as a function of the coefficients Fq for a given H ∗ . Using Eq. (B.7) and integrating gives the average pressure hpn (x, z)iz on the plate for a given Galerkin mode yn . The results presented in the present paper were obtained for a Chebyshev truncation M = 20 and Galerkin truncation N = 10. It is important to have a sufficient number of modes in the Galerkin discretisation if one wants to model accurately the instability threshold (especially for large mass ratio M ∗ ) as it has already been shown by Lemaitre et al. (2005). C Slender body approximation For small k x , the average of the pressure along the span is obtained from Eq. (32) and is simply hpn iz = πk x /4. Coming back to the dimensionless variables kx and z, ven can be expressed as ven (kx ) = − 8 hpen (kx , z)iz when H ∗ → 0. ∗ πH (C.1) Following the same procedure as between Eqs. (38) and (41), one finds H(X, H ∗ ) = − 8 sgn(X) when H ∗ → 0. H∗ 21 (C.2) c in the Appendix B. The integral equation This corresponds to the term H 2 (40) is then solved by integrating by parts, leading to ∂yn ∂ 2 yn πH ∗ hpn (x, z)iz = − −ω 2 yn + 2iω + 4 ∂x ∂x2 ! when H ∗ → 0. (C.3) This result is the same as the one obtained by Lighthill (1960) and used by Lemaitre et al. (2005) in their study (they added an additional tensile effect due to tension because their ribbon was hanged vertically). This latter expression can be inserted into Eq. (11) to obtain the eigenmodes and eigenfrequencies for the slender body approximation. This approximation is valid when the wavelength of the modes is large compared to the plate span. D Large span approximation In the limit of large aspect ratio, k x is large and the average of the pressure along the span is obtained from Eq. (34). Following the same procedure as between Eqs. (38) and (41), one finds H(X, H ∗ ) = −  2 1 π − ∗ sgn(X) + O X H H ∗2  when H ∗ → ∞. (D.1) In the new variables defined in Appendix B, we have to solve the singular c is now integral equation (B.2) where the function H c X, c H ∗) = H c + H( 1  1 π c H2 + O ∗ 8H H ∗2  when H ∗ → ∞. (D.2) The linear system (B.8) has now to be solved with the new matrix ∗ C(H ) = C (1)   π (2) 1 . + C +O ∗ 8H H ∗2 (D.3) We have to remember that we are interested in determining the matrix P given by Eqs. (12) and (43) such that P = −ω 2 P (M ) + 2iωP (G) + P (K) . (D.4) As in Eq. (43), the terms on the right hand side correspond to the added mass, gyroscopicity and stiffness pressure respectively. In this large span approximation, these matrices can be obtained easily if the scalar products are all evaluated on the basis of the Chebyshev polynomials of the second kind. Indeed, once the following scalar products have been evaluated numerically 22 (M ) Ym,q = Uq−1 ⊙ ym , dym (G) , Ym,q = Uq−1 ⊙ dx d2 y m (K) , Ym,q = Uq−1 ⊙ dx2 (D.5) (D.6) (D.7) the matrices P (X) (where the superscript X stands for M , G or K) can be obtained easily after some straightforward algebra. They are expressed as powers of 1/H ∗ such that P (X)     π 1 (X) 1 (X) =− P0 + ∗ P1 +O , 4 H H ∗2 (D.8) where M X 1 (X) m,n = (X) m,n =− P0 P1 q=1 q (M ) (X) Ym,q Yn,q , (D.9) M M πX 1 (M ) X (2) (X) Ym,q Ci,q Yn,i , 8 q=1 q i=1 (D.10) The above expressions are very easy to obtain since Eqs. (B.11–B.13) give an explicit expression for C (2) . An analytical expression of P (X) accurate up to order 1/H ∗ 2 could be obtained by the same method as we have shown in a previous paper (Eloy et al., 2006). However, the improvement is not significant and the exact calculus detailed in Appendix B is still necessary for moderate (X) H ∗ . Note that when P1 is set to zero, or equivalently for infinitely large aspect ratio H ∗ , we recover the results of the two-dimensional theory of Guo and Paı̈doussis (2000). References Argentina, M., Mahadevan, L., 2005. Fluid-flow-induced flutter of a flag. Proc. Nat. Acad. Sc. USA 102 (6), 1829–1834. Balint, T. S., Lucey, A. D., 2005. Instability of a cantilevered flexible plate in viscous channel flow. Journal of Fluids and Structures 20, 893–912. Crighton, D. G., 1985. The Kutta condition in unsteady flow. Ann. Rev. Fluid Mech. 17, 411–445. Datta, S. K., Gottenberg, W. G., 1975. Instability of an elastic strip hanging in an airstream. Journal of Applied Mechanics, Transactions of the ASME, 195–198. Eloy, C., Souilliez, C., Schouveiler, L., 2006. Flutter of a rectangular cantilevered plate. Proceeding of the 6th FSI, AE & FIV+N Symposium, ASME PVP 2006 / ICPVT-11 Conference. 23 Frederiks, W., Hilberink, H. C. J., Sparenberg, J. A., 1986. On the Kutta condition for the flow along a semi-infinite elastic plate. J. Eng. Math. 20, 27–50. Guo, C. Q., Paı̈doussis, M. P., 2000. Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow. Journal of Applied Mechanics, Transactions of the ASME 67, 171–176. Huang, L., 1995. Flutter of cantilevered plates in axial flow. Journal of Fluids and Structures 9, 127–147. Kornecki, A., Dowell, E. H., O’Brien, J., 1976. On the aeroelastic instability of two-dimensional panes in uniform incompressible flow. Journal of Sound and Vibration 47 (2), 163–178. Lemaitre, C., Hémon, P., de Langre, E., 2005. Instability of a long ribbon hanging in axial air flow. Journal of Fluids and Structures 20 (7), 913–925. Lighthill, M. J., 1960. Note on the swimming of slender fish. J. Fluid Mech. 9, 305–317. Lucey, A. D., Carpenter, P. W., 1993. The hydroelastic stability of threedimensional disturbances of a finite compliant wall. Journal of Sound and Vibration 165 (3), 527–552. Paı̈doussis, M. P., 1998. Fluid-Structure Interactions: Slender Structures And Axial Flow, Volume 1. Academic Press. Paı̈doussis, M. P., 2004. Fluid-Structure Interactions: Slender Structures And Axial Flow, Volume 2. Elsevier Academic Press. Rayleigh, L., 1879. On the instability of jets. Proc. London Math. Soc. X, 4–13. Shayo, L. W., 1980. The stability of cantilever panels in uniform incompressible flow. Journal of Sound and Vibration 68 (3), 341–350. Shelley, M., Vandenberghe, N., Zhang, J., 2005. Heavy flags undergo spontaneous oscillations in flowing water. Phys. Rev. Lett. 94, 094302. Souilliez, C., Schouveiler, L., Eloy, C., 2006. An experimental study of flag flutter. Proceeding of the 6th FSI, AE & FIV+N Symposium, ASME PVP 2006 / ICPVT-11 Conference. Taneda, S., 1968. Waving motions of flags. Journal of the Physical Society of Japan 24 (2), 392–401. Tang, L., Paı̈doussis, M. P., 2006. A numerical investigation on the dynamics of two-dimensional cantilevered flexible plates in axial flow. Proceeding of the 6th FSI, AE & FIV+N Symposium, ASME PVP 2006 / ICPVT-11 Conference. Theodorsen, T., 1935. General theory of aerodynamic instability and the mechanism of flutter. NACA Report 496. Watanabe, Y., Isogai, K., Suzuki, S., Sugihara, 2002a. A theoretical study of paper flutter. Journal of Fluids and Structures 16 (4), 543–560. Watanabe, Y., Suzuki, S., Sugihara, M., Sueoka, Y., 2002b. An experimental study of paper flutter. Journal of Fluids and Structures 16 (4), 529–542. Yadykin, Y., Tenetov, V., Levin, D., 2003. The added mass of a flexible plate oscillating in a fluid. Journal of Fluids and Structures 17, 115–123. 24 Yamaguchi, N., Sekiguchi, T., Yokota, K., Tsujimoto, Y., 2000. Flutter limits and behavior of a flexible thin sheet in high-speed flow - ii: Experimental results and predicted behaviors for low mass ratios. Transaction of the ASME Journal of Fluids Engineering 122, 74–83. 25