_53-65_ Kusay H. Jabir 53 Flutter Speed Limits of Cantilever Rectangular and Tapered Plates Kusay H. Jabir Foundation Technica l Education Technical College-Baghdad (Received 3 April 2007; accepted 9 September 2007) Abstract: The aerodynamic and elastic forces may cause an oscillation of the structure such as the high frequency of the airfoil surfaces and the dynamic instability occurring in an aircraft in flight and failure may occur at a speed called flutter speed. In this work, analytical and numerical investigations of flutter limits of thin plates have been carried out. The flutter speed of rectangular plates were obtained and compared with some published results. Different design parameters were investigated such as aspect ratio, thickness and their effects on flutter velocity. It was found that the structural mode shape plays an important role in the determination of the flutter speed and the coupling between the bending and torsional mode is the main cause of flutter. Key words:Flutter, Speed, Plates. Introduction: In mechanics, there are three types of forces aerodynamic, elastic and inertia. The interaction between aeroelasticity and elastic forces is called static aerodynamic such as the divergence problem and the interaction between elastic, aerodynamic and the inertiaforces called the dynamic aero elasticity and the flutter is an example for this phenomena in which the amplitude of the oscillations may diverge causing failure. Therefore, this research is concerned with the determination of flutter speed for thin plates which are found in many aeronoutical structures. R.S Srinivasan and B.J. Baba [2] studied the flutter analysis of cantilevered quadrilateral plates and the problem is solved by using a numerical method involving an integral equation based upon finding the strain and kinetic energies of the plate. E. Nissim and I. Lottat [3] suggested an optimization method for the determination of the important flutter modes. The method is based on the minimization of the quadratic values derived from the equation of motion. Al- Khwarizmi Engineering Journal Vol.3 , No.3, pp 53-65 (2007) Al- Khwarizmi Engineering Journal Kusay H. Jabir/Al-khwarizmi Engineering Journal, Vol.3, No.3, pp 53 - 65 (2007) 54 R.S. Srinivasan and B.J. Baba [4] described the free vibration and flutter of laminated quadrilateral plates. They derived the differential equation required to obtain the flutter speed of simple rectangular plates. M.W. Kehoe [5] reviewed the test techniques developed over the last several decades for flight flutter testing of aircraft practical experiences and example test programs were presented to test the effectiveness of the various approaches used. Brain Danowsky [6] studied the development of an aircraft structural variation model, accounting for variations in structural mode shape as well as structural frequency, which has many advantages for the design of aircraft and aircraft flight control systems. J.Hoas, etal [7] described the flutter of circulation control wings and explained the flutter stability of high aspect ratio circulation control wing using lumped approach on conjunction with a modified unsteady aerodynamic strip analysis method. J.E. Sedaghat et al [8] developed to predict the speed and frequency at which flutter occur based upon the use of symbolic programming. The approach performs the computation in a single step and does not require the repeated calculations at various speeds required when using the classical V.G. method. T.S Talib [9] determined analytically using the damping method and calculated the natural frequencies and mode shapes using the finite element technique for subsonic wing structure. He showed that the flutter speed change with changing the skin thickness, material properties and the altitude. His calculations indicate that structural mode shape variation plays a significant role in the determination of wing flutter limits. In this work, the V.G. method together with developed finite element package will be adopted and the flutter speed will be defined for thin plates with different aspect ratios and thickness for both rectangular and tapered plates. Theory: To determine the equation of motion, consider the section shown in Fig.1 with two degrees of freedom α and h. [9] The downward displacement of any other point on the airfoil is: [8] Z=h±xα The strain energy is composed of two components linear and rotational parts 22 . 2 1 2 1 hKKU h+= αα ………………..……(1) And the kinetic energy is: ∫ ∫ ∫++= ]2[2 1 222 dxxxdxhdxhT ραραρ ... (2) Using: ∫ == 22 amrdxxI ρα ∫ == amrxdxS ρα Therefore, 22 . 2 1 2 1 αα αIhmxmhT a ++= Using Lagrange equation in the form: qQq UT q UT dt d = ∂ −∂ − ∂ −∂ = )( ] )( [ & ………… (3) Where q=q (h,α) Inserting equations (1) and (2) into (3) yield the following characteristic equation,         =               +               2 22 2 0 01 mb M mb F b h rW W b h rx x h a a αα ααα && && .. (4) The solution may be assumed as follows: iwtehh 0= And iwte0αα = ……………..…………………(5) equation (4) is reduced to Kusay H. Jabir/Al-khwarizmi Engineering Journal, Vol.3, No.3, pp 53 - 65 (2007) 55         =               +               2 22 2 2 2 0 01 mb M mb F b h rw wb h xx x w h αα αααα α …. (6) The aerodynamic force and moment are given by [9]: )]) 2 1 (([23 nn LaLLb h wbF +−+= ααπρ Where K ciLh 1 21 −= 2 221 2 1 K C K C iL − + −=α }]) 2 1 ())( 2 1 ({ ) 2 1 ([ 2 24 hh hh L a MLaM b h LaMwbM + +++−+       +−= αα πρ k iMM h 1 8 3 , 2 1 −== α Then the equation of motion becomes:               +               −= α α αα αα α b h rw w b h rx x wL h h 22 2 2 2 0 0 1 Simplifying and using 2 2 2 α η w w = and 2 2 2 αw w R h= with 2b m πρ µ = gives:                   ++++−+− +−Ω = αµ αα α b h LaMLaMaM LaLL L hhh hh h 2 2 ) 2 1 ())( 2 1 () 2 1 ( ) 2 1 ( …(7) which can be expressed as a flutter equation in the form:         +Ω=         αα b h MAb h K ijijij ][][ 2 Fig. 1 Airfoil restraind from bending and torsional motion in an airstreem by spring Kh and Kα, acting a distance ba aft of midehord. Also shown are lift L and pithing moment My about the axis twist.[9] Where [Kij] is the stiffness matrix [Mij] is the mass matrix and [Aij] is the aerodynamic matrix Finite Element Method The plate is analyzed by using a suitable element. The element have 8 degrees of freedom with 5 degrees of freedom of each node (u, v, w, θx and θy in the x, y and z directions) and the rotations about x and y axes respectively. The generic displacements at any point of the middle surface are [10]: ∑ ∑ ∑ = = = = −=−= == = 9 1 9 1 9 1 2 2 },,{}{ i ii i xiix i yiiy wNw NZv NZuand wvu θθ θθ δ ……….……….…….(8) ……………..……….(9) …………………….(10) Kusay H. Jabir/Al-khwarizmi Engineering Journal, Vol.3, No.3, pp 53 - 65 (2007) 56 { }∈ { }                 + + +=                 ∈ ∈ =∈ zx yz xy y x zx yz xy y x uw wv vu v u ,, ,, ,, , , γ γ γ [ ]{ }δB )9,...2,1( 0 0 22 0 0 2 0 2 00 =                       − − i Na Nb b h a h b h a h B ii ii ii i i i ζζ ζ ζ 2 h ζ ...... 2 BiAii B h BB ζ+=                 − −− − =                 − − = i ii i Bi ii ii Ai N ba b B Na Nb B 00 000 0 00 000 , 0 0 000 000 000 ηξ ηξ ,22*,12* ,12*,11* NNJb NNJa i i += +=           = ξ ηξ ηξ ,00 0,, 0,, ][ Z YX YX D ∫ ∫ ∫= ζηξ dddJBDBK Te ]][[][][ ηξddJBDB h BDBK B T BA T A e ∫ ∫ += ]][[][6]][[][2][ 2       = shear bending D D D 0 0 ][       + =             −− 10 01 )1)(2.1(2 ; 2 1 00 01 01 )1(12 2 3 v Eh D v v v v Eh D shear bending ηξρ ddJNN h NNM B T BA T A e ) 6 2(][ 2 += ∫ ∫ Strain–Displacement Relationship : The strain matrix is given by : ...........................(11) = And the strain displacement matrix [B] may be written as: ………….…(12) We can isolate terms in bi that multiply by , then And J* is Jacobian inverse matrix. ………..(13) Formulation of Stiffness Matrix : The element stiffness matrix [K]e for the element, using B matrix in eq. (12) as follows : …...(14) Taking integration through thickness; eq . (14) becomes Where {D} is the stress-strain matrix for isotropic material ………………………..(15) The first part of eq (15) is due to transverse shearing deformation where as the second part is associated with flexural deformation. Formulation of Mass Matrix : The formulation of consistent mass matrix for element QPB9 becomes (17) Kusay H. Jabir/Al-khwarizmi Engineering Journal, Vol.3, No.3, pp 53 - 65 (2007) 57 NiNNiN BA           −=           = 001 010 100 , 001 000 000 ∫ ∫ ∫ ∂∂= ζηξξ dddJNNA Te ]][[][][ ∫ ∫ ∫ ∂∂+= ξξ /)()(2]([][ A T BA e NN h NA ∫ ∫ ∂∂= ]/[][2(][ ξATAe NNA AjiAij ≠ { } 0])[][][ =−+− QAKMK λ λ D b V ro 3 2 2 1 ρλ = )/(422 Dtbrwk ssρ= 2w Where The first part of eq. (17) consists of translation inertias and the second part gives rotational (or rotary) inertias. Formulation of Aerodynamic Matrix The Formulation of aerodynamic Matrix for element QPB9 becomes. ζηξηξ dddJN h B ])/[2 ∂∂+ ……….(18) For simplicity it is assumed that the air flowing above the panles is parallel to the X-axis as shown in fig. 2 and the effect of any air entrapped below may be neglected . By integrating Eq . (18) through the thickness produces. ηξη ddJNN h B T B ]/[][2 2 ∂∂+ ………...(19) Where Fig.2 panal under airfolow. Flutter Equation and Solution Method : The flutter equation can be derived as: …………….(20) Where K, M and A are matrices of size (m*m) degrees of freedom. Here Non dimentional dynamic parameter and it is equal to ……………………(21) Where pa is air density, V is the free stream velocity, br is the half length of wing at root, D is the flexural rigidity of the plate and, k2 is the Eigen value and it is equal to [8] ………….…(22) Where is the natural frequency, ρs is the material density and ts is the plate thickness. Equation (20) is a standard eigen value problem. The parameters and k are non-dimensional quantities. Note that for zero flow velocity and k then represent the square of the non-dimensional natural frequency, so that in this case the eigenvalues are real. As the flow velocity increases from zero, two eigenvalues will usually approach each other and coalesce (become complex conjugates) to kcr at a value of 2 which is the critical value of dynamic pressure. For a large system of equation, it is a relatively difficult task to find the eigen values of equation (12), since in general the matrices involved are complex and asymmetric. Therefore, the problem is separated into two parts. First, set zero flow veiocity (w=0) in equation (12) and determine a finite number, say n,of normal modes of the structure ,where n