Preparation of Papers in a Two Column Model Paper Format _________________________________________ *Corresponding author e-mail: esalimipour@qiet.ac.ir ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 41 STALL FLUTTER CONTROL OF A WING SECTION BY LEADING EDGE MODIFICATIONS E. Salimipour 1* , M. Saei Moghaddam 2 , Sh. Yazdani 3 1,3 Department of Mechanical Engineering, Quchan University of Advanced Technology, Quchan, Iran 2 Department of Chemical Engineering, Quchan University of Advanced Technology, Quchan, Iran ABSTRACT A solution procedure is described for determining the two-dimensional and two- degrees of freedom flutter characteristics for wings at large angles of attack. This procedure requires a simultaneous integration in time of the solid and fluid equations of motion. The fluid relations of motion are the unsteady, compressible Navier-Stokes equations, solved implicitly by second-order Roe’s approximation scheme in a moving coordinate system. The solid equations of motion were integrated in time by use of fourth-order Runge-kutta method. In this paper, the stall flutter of a rectangular wing with section of NACA 0012 is studied. Therefore, the aeroelastic responses for the system were calculated by applying mode shapes for vibrating wing. Then the obtained responses resulted from several changes in leading edge shape of wing are compared. Results showed that these different leading edge shapes cause the changes on oscillating parameters of the system. In these changes, applying a camber with 25 o angle had the best result in this study. KEYWORDS: Stall flutter; Navier-Stokes equations; Leading edge shape; Aeroelastic responses; Mode shapes. 1.0 INTRODUCTION Aeroelasticity is defined as the interaction of aerodynamics, elasticity and dynamics. Classical theories of aeroelasticity assume that the aerodynamic and structural forces are linear. For many decades, the classical approach has been successful in providing approximate estimates of aircraft response to gusts, turbulence and external excitation. The flutter boundaries are often quite accurately predicted when compared to flight test results. On the other hand, these classical methods are unable to capture phenomena arising from structural and aerodynamic nonlinearities. Aerodynamic nonlinearities are often encountered at transonic speeds or high angles of attack where flow separation occurs (Ghadiri & Razi, 2007). Flutter phenomenon is a kind of dynamic instability that results from the interaction of inertia, mailto:esalimipour@qiet.ac.ir Journal of Mechanical Engineering and Technology 42 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 elastic and aerodynamic forces and causes the vibration of the wing to diverge. At transonic and supersonic regime, the shock wave incidence affects the oscillating parameters of the system, because this event can alter generalized forces, which operate on the system and at high angles of attack the flutter encounters a phenomenon which is called dynamic stall. Dynamic stall is a phenomenon caused by vortex shedding on the surface of oscillating airfoils at high angles of attack. This causes a huge decrease in lift and increase in drag force and pitching moment too. If, during part or all of the time of oscillation, the flow was separated, then the flutter phenomenon exhibits some different characteristics and is called stall flutter (Fung, 2007). In past decades existing analytical tools for the prediction of stall flutter treated this problem as a fluid-structural interaction problem. Since the unsteady aerodynamics associated with this phenomenon is complex, researchers have in the past relied on experimental data for airfoil static and dynamic characteristics. These data are usually synthesized from a series of analytical expressions for different parts of the dynamic stall loop and used in a flutter analysis (Jiunn-Chi, Kaza & Sankar, 1987). The fluid-structural equations can be integrated in time simultaneously. The aerodynamic loads determined from the integration of the unsteady, compressible Navier-Stokes equations, drive the structural dynamics equations. The results of structural deformations alter the aerodynamic loads. Under certain conditions, the coupling between the aerodynamic loads and the structural motion by imposing small disturbances on the wing angle of attack cause the aeroelastic responses grow rapidly in a divergent oscillatory fashion, and finally, flutter occurs. There are many studies of stall flutter at high Reynolds numbers starting from about 1950. Maybe, the most detailed of the former investigations are those of Halfman, Johnson & Haley (1951) and Rainey (1958). In these studies, the wing was placed at a mean angle of attack and then forced to vibrate, either in pitch or heave. One of the first stall flutter analyses using numerical solution of navier-stokes equations carried out by Jiunn-Chi et al. (1987). They studied NACA 0012 airfoil flutter by simultaneous integration in time of the solid and fluid equations of motion at high angles of attack. Price & Keleris (1955) investigated NACA 0012 airfoil flutter by applying nonlinear effects of aerodynamic loads at stall angle of attack by using semi-experimental methods. More recent studies such as Razak, Andrianne & Dimitriadis (2011) have performed PIV 1 measurements for an elastically mounted wing undergoing stall flutter. The wing was free to move in pitch and Heave, but they observed that in stall flutter the pitching mode was predominant. An analysis based on modified Leishman-Beddoes model at low mach number was carried out by Song, Qinghua, Chenglin & Xianping (2011). The main modifications for L-B model included a new dynamic stall criterion and revisions of normal force and pitching moment coefficient. Bhat & Govardhan (2013) experimentally studied the stall flutter boundaries of a NACA 0012 airfoil at low Reynolds numbers by measuring the forces and flow fields around the airfoil when it is forced to oscillate. These measurements indicated that for large mean angles of attack of the airfoil, there is positive energy transfer to the airfoil over a range of reduced frequencies, indicating that there is a possibility of airfoil excitation or stall flutter. Sun, Haghighat, Liu & Bai (2015) developed a nonlinear time-domain aeroservoelastic model to study stall flutter and design flutter suppression control systems. A review 1 Particle Image Velocimetry Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 43 research for Nonlinear flutter wind tunnel test and numerical analysis of folding fins with freeplay nonlinearities was performed by Ning, Nan, Xin & Wei (2016). The scope of the present work is to find a way for delaying of stall flutter effects by limiting the vibration amplitudes of the wing. One of the ways is the eliminating of stall by moving forward the leading edge separation point. By modifying the leading edge shape, we can alter the position of flow separation on the airfoil surface. In this study, Raynolds Average Navier-Stokes (RANS) equations with two-DOF structural dynamics model of wing are solved simultaneously in a moving computational grid by writing a computer code. The RANS equations are applied on a full unsteady turbulent compressible flow. For this purpose, the turbulent shear stresses were modeled by using an eddy viscosity concept. In this study, the Spalart-Allmaras one-equation model is used because of its adequate accuracy for flows with adverse pressure gradient (Blazek, 2001). Generally, the solution procedure is that the aerodynamic loads obtained from the solution of the RANS equations were placed into the solid equations of motion in each time step. Then, the aeroelastic responses and the new velocities of moving grid will be achieved. These values will be used in the fluid flow equations for new time step. Since there is a need to carefully validate any aerodynamics model prior to its application to flutter calculations, a number of code validation studies are first presented, including some cases where deep dynamic stall and stall flutter occurs. Following the code validation, the studies on leading edge modifications are presented. 2.0 MATEMATICAL AND NUMERICAL FORMULATION 2.1 Flow-field Equations In this section, according to Jiunn-Chi et al. (1987), the numerical procedure used to compute the unsteady viscous flow is briefly described. In the (x,y,t) coordinate system, the two-dimensional, unsteady Navier-Stokes equations may be written as v v F GU F G t x y x y              (1) Where U is conservative variables vector, F and G are convective flux vectors and Fv and Gv denote viscous term vectors, which may be written as     2 2 0 0 , , , , xx xy v v yx yy xx xy yx yy u v u p uvu U F G F G uv v pv u v u vu E p v E pE                                                                          (2) E is the total energy per unit volume and is defined as Journal of Mechanical Engineering and Technology 44 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 2 2 1 2 p u v E            (3) The shear stress components are expressed as   4 2 3 3 xx t u v x y              (4)   4 2 3 3 yy t v u y x              (5)  xy yx t u v y x                (6) where μ and μt represent the dynamic viscosity and turbulent addy viscosity, respectively. All the calculations were performed in a transformed coordinate system (ξ,η,τ), which is linked to the physical plane (x,y,t) according to the following relationship:    , , ; , , ;x y t x y t t       (7) The Jacobian of transformation J is given by 1 J x y y x       (8) and the metrics of transformation are given by    ; ; ; ; ;x x t tz z t t t ty x y x x y x y x y x y J J J J J J                         (9) ξt and ηt include the grid points velocities in the coordinate system. In the (ξ,η,τ) coordinate system at non-dimensional form, the two-dimensional, unsteady Navier-Stokes equations may be written as * ** * * Re v v F GMU F G                       (10) where M∞ and Re denote the free-stream mach number and Reynolds number, respectively. the quantities F * , G * , Fν * and Gν * are given by Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 45         * * * * * / ; / ; / ; / ; / x y t x y t v x v y v v x v y v U U J F F G U J G F G U J F F G J G F G J                      (11) Above equations are discretized using an implicit, finite-volume scheme based on Roe’s second-order approximate Riemann solver (Blazek, 2001). 2.2 Dynamic equations of wing Resisting elastic forces are developed that are proportional to the airfoil torsional and translational displacement. In Figure 1, a sketch of the two-DOF rectangular cantilever wing is shown. Note that the pitching axis may offset from the center of mass of the wing, leading to a coupling between the pitching and heaving degrees of freedom. The governing equations of the two-DOF system motion are (Fung, 1945): 2 2 4 2 2 4h h h h m S C EI L t t t y               (12) 2 2 2 .2 2 2 E A h I S C GJ M t t t y                  (13) where y is the coordinate axis along the wing, h is the wing vertical displacement, α the angular displacement (angle of attack), Iα the wing section moment of inertia (per unit span) about the pitching axis, m and Sα=mxα the mass and static moment per unit span, and Ch and Cα the structural damping for the pitching and plunging motion coefficients, respectively. The EI and GJ are the torsional and bending stiffnesses, respectively. The Sα depends on the offset between the wing pitching axis and the wing center of gravity. The lift force (L) and pitching moment (ME.A) are calculated by following relationships:     dspL yyyxyx   (14)    dLxxM AEAE .. (15) where ds is the surface element. p, τxy and τyy are calculated by solving the unsteady, compressible Navier-Stokes Equation (10). Figure 1. Schematic diagram of the two-DOF system Journal of Mechanical Engineering and Technology 46 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 The partial differential Equations (12) and (13) are solved based on Galerkin’s method. For this purpose, the h and α can be decomposed to product of two independent variables as follows:   1( )h h t  (16)   1( )t    (17) More details of above method are found in (Hodges & Pierce, 2002). Φ and Ψ are the torsional and bending mode shapes for vibrating wing and can be written as (Marzocca, Librescu & Silva, 2002); 1 1 1 1 1 1 1 1 1 sinh sin (cos( ) cosh( ) sinh( ) sin( ) cosh cos K                        (18) 2 2 sin( )K    (19) where K1 and K2 are constant coefficients and, β1=0.5969π, β2=0.5π and η=y/l. By replacing the Equations (16) and (17) into the Equations (12) and (13) the non- dimensional forms of wing equations of motion can be written as below: 2 3 6 3 1 5 1 1 3 12 3 2 h h A C A L A H A x H A H b b                           (20) 2 24 7 . 4 1 5 1 1 4 14 4 2 E A A C A M A r A x H A r b b                       (21) where: 1 1 2 1 4 2 2 2 3 4 ; ; ; ; ; h Ih A EI A GJ m H t r b A ml A I l b mb                  (22) The symbol apostrophe used in Equations (20) and (21) represents derivative with respect to τ and the coefficient A1 to A7 are constant values that obtained from vibrating mode shapes where are determined in appendix. In order to reach aeroelastic responses, first the lift force and pitching moment will be obtained at each time step by solving the Navier-Stokes equations, then by placing them into the Equations (20) and (21) the new values of h1 and α1 will be calculated. Consequently, the flow conditions will be determined for next time step. We can see that the variables fo Equations (12) and (13) will be updated at each time step. Thus, these equations are fully unsteady. Figure 2 shows a flowchart of this solution procedure. In the present work, the Equations (20) and (21) era integrated in time using a fourth-order Runge-kutta method and the aeroelastic responses are calculated at the end of the wing where y=l (or η=1). Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 47 Figure 2. Solution procedure flowchart 3.0 RESULTS AND DISCUSSIONS 3.1 Code Validation Study Before the application of this computer code in flutter applications, the flow solver was extensively calibrated for a number of test cases. In this work, only a few of these studies are reported. A 575×81 C-type grid with excellent orthogonality was used in this study which is shown in Figure 3. Figure 3. View of C-type grid used in flow computations In Figure 4, the surface pressure distribution for a NACA 0012 airfoil at 11.0264 deg angle of attack is shown. The free-stream Mach number and Reynolds number were Solve RANS equations and obtain L and ME.A Put L and ME.A in structural equations Solve structural equations and obtain h, α, Add α to wing angle of attack Obtain computational grid velocities using Go to new time step Journal of Mechanical Engineering and Technology 48 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 0.301 and 3×10 6 , respectively. For comparison, experimental data obtained by Charles, Acquilla & William (1987) are also shown. The present results and the experiments are seen to be in good agreement everywhere. Figure 4. Comparisons of theory versus experiment for the surface distribution over a NACA 0012 airfoil at 11.0264 deg angle of attack (M∞ = 0.3, Re = 3×10 6 ) In order to illustrate the capability of the Navier-Stokes solver to obtain time-accurate results in highly separated flow, the lift, drag, and moment hysteresis loops are shown and compared with experiments (McAlister, Carr & McCorosky, 1982) in Figure 5 for a NACA 0012 airfoil oscillating in pitch. The mean angle and amplitude of oscillation were 14.91 and 9.88 deg, respectively. The reduced frequency of oscillation, normalized with respect to the half-chord, was k=0.151. The free-stream Mach number and Reynolds number were 0.283 and 3.45×10 6 , respectively. It is seen that the theory correctly predicts the near linear increase in lift during the upstroke, the dynamic stall causing rapid variations in lift, drag, and moment alike, and the poststall recovery phase of the flow during the downstroke. The fact that the flow solver is able to capture much of the dynamic stall flow features increases the confidence in the capability of this code to handle stall flutter problems. It must be mentioned that use of wall function approach for turbulence modeling could not be accurately reliable in this study. Because the law of wall doesn’t always hold for flow near solid boundaries, most notably for separated flows (Wilcox, 1994). Experiment (Charles et al., 1987) Present study Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 49 Figure 5. Comparisons of theory and experiments for the unsteady airloads on a NACA 0012 airfoil experiencing dynamic stall (M∞ = 0.283, Re = 3.45×10 6 , k = 0.151). As a final test of the above solver's ability to handle aeroelastic responses manner, the stall flutter calculations by use of the Navier-Stokes/structural dynamics solver explained above were considered. In this test case, the starting point was the steady viscous flow over a NACA 0012 airfoil at 0.3 Mach number and 9×10 6 Reynolds number at 15 o angle of attack. At this angle of attack, the airfoil is on the verge of stall. The airfoil was given a small amplitude perturbation in its angle of attack and the subsequent motion was obtained. The pitch and heave responses are plotted as a function of non-dimensional time and are compared with obtained results performed by Jiunn-Ch et al. (1987) as shown in Figure 6. Very good agreement is observed between the two solvers. It was found that the airfoil returned to steady state following a period of damped oscillations. The non-dimensional aeroelastic parameters are assigned as below: 0 0 0.2, 0.25, 100, 51.5, 0.5, 15 , 0.25 h x r h b              Figure 6. Time response of a two-degree of freedom solid-fluid system experiencing stall flutter (NACA 0012 airfoil, M∞ = 0.3, Re = 9×10 6 , initial angle of attack 15 o ) 3.2 Flutter Calculations We use the current unsteady solver in the mentioned fluid-structure interaction method for the two-dimensional NACA 0012 wing section. This model simulates the bending and torsional motion of the wing. It consists of two degrees of freedom, heaving and pitching, for a NACA 0012 wing section. Figure 7 shows the variation of flutter speed Jiunn-Chi et al. (1987) Present Experiment (McAliste et al., 1987) Present study Journal of Mechanical Engineering and Technology 50 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 index versus the freestream Mach number (M∞) for inviscid flow. Flutter speed index includes the parameters which are affected on flutter incidence and is defined as (Kirshman & Liu, 2006)     f f f U V b      (23) where (U∞)f and (ωα)f are the values leading to happening flutter. The non-dimensional aeroelastic parameters are assigned as 0 0 0.5; 0.25; 26.35; 0.5; 1 ; 0.1 h x r h b           It is seen that by increasing the freestream Mach number, the flutter speed is reduced. Also, at transonic speeds the curve slop is increased and clearly, the system stability will be lower than subsonic one. Figure 7. Variation of flutter speed index versus Mach number 3.3 Leading edge modifications In this section, we attempt to find a way to limit the stall flutter phenomenon. The most effective factor in this kind of flutter is dynamic stall incidence, which alters the aerodynamic loads. Dynamic stall is a phenomenon caused by vortex shedding on the surface of oscillating wings at high angles of attack. This causes a huge decrease in lift and increase in drag force and pitching moment. Therefore, by eliminating of dynamic stall, we can suppress the stall flutter. One of the most effective portions of the wings on flow separation is the leading edge. By modifying the leading edge shape, we can alter the position of flow separation on the wing surface. We applied several modifications on the leading edge shape of the NACA 0012 wing section and compared the aeroelastic Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 51 responses, lift and pitching moment coefficients with the standard case. The following values are selected for the non-dimensional aeroelastic parameters: 6 0 0 0.5; 0.25; 100; 51.5; 0.5 0.4; Re 2.5 10 ; 15 ; 0.25 h x r M H                     In the first case, we increased the sharpness of the leading edge and then solved the aeroelastic equations. In Figure 8, the results are compared to the standard leading edge shape. The vibration amplitude has been reduced in heave motion; although the pitch amplitude is rising in compare to the standard case responses. Therefore, this change could not be suitable. Figure 8. Comparisons of sharped versus standard leading edge for aeroelastic responses The second case involves a comparison of blunted versus standard NACA 0012 leading edge for the aeroelastic responses with same parameters as shown in Figure 9. We can see that by increasing the leading edge thickness, the vibrations amplitudes have been descended. The reason is that the leading of the relocation of separation onset toward Journal of Mechanical Engineering and Technology 52 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 the trailing edge as shown in Figure 10. By moving forward the separation onset point, the vortex growth opportunity will be reduced. Also, the separation in blunted leading edge occurs at the higher angle of attack than the standard one; consequently, the stall will exist in smaller regions. Figure 9. Comparisons of blunted versus standard leading edge for aeroelastic responses 0.07 c 0.14 c α=21.5 ᵒ α=24.5 ᵒ Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 53 Figure 10. Comparisons of blunted versus standard leading edge for separation onset Figure 11. Comparisons of small cambered versus standard leading edge for aeroelastic responses For the next case, a small camber in leading edge is applied (Carr, McAlister & McCorosky, 1977). In Figure 11, the results with same parameters are compared with the standard leading edge shape. It is seen that the amplitudes are decreased. Because, in this case, the stall depth has been reduced. At the last case, we applied a camber with 25 o angle on the leading edge where in Figure 12, the results with same parameters are compared with the standard leading edge shape. By observing the results, the reduction of amplitudes is determined. Figures 13 and 14 show the flow field for some positions of the wings in a period of vibrations in two cases of standard and 25 o cambered leading edge. It can be seen flow separation condition changes by cambering of leading edge Angle of attack as positive and negative, changes in each period on these vibrations. Since in two last cases downward cambering of leading edge tend to delay in flow separation on positive angles of attack and vice versa. Since the initial angle of attack is 15 o , there is not any separation at the beginning unless angle of attack was negative, flow separation had occurred slightly; however, when it is positive again, separation disappeared. It can be observed that there is not any deep stall in the flow field until the end of vibrations. Journal of Mechanical Engineering and Technology 54 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 Figure 12. Comparisons of 25 o cambered versus standard leading edge for aeroelastic responses. Figure 13. Flow field for Standard NACA 0012. Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 55 Figure 14. Flow field for 25 o cambered NACA 0012 leading edge CONCLUSIONS The turbulent compressible flows were solved in the time domain. With existence of dynamic stall, the amplitude and frequency of responses has been affected. It was observed that by producing a camber on the leading edge the flutter onset had been delayed. Because the stall depth was reduced at lift force and pitching moment. Specially applying a camber with 25 o angle had the best result in this study. It should be noted that the changes of leading edge could affect in flutter decreasing greatly, but in some cases (for example: cambering over of leading edge) these changes could eliminate aerodynamic properties of wing. APPENDIX A By replacing the Equations (16) and (17) into the Equations (12) and (13), based on Galerkin's method, the space dependent terms will appear to integral relations and respect to the wing geometry, will have constant values where form coefficients A1 to A7. These values have been obtained as   2 2 1 2 1 12 0 22.94429 d A d K d              (A.1)   2 21 2 2 2 0 8 d A d K d               (A.2)    21 2 3 1 0 1.85598A d K    (A.3)    2 21 2 4 0 2 K A d    (A.4)     1 5 1 2 0 0.92348A d K K       (A.5)   1 6 1 0 1.06667A d K     (A.6)   1 7 2 0 0.63662A d K    (A.7) Journal of Mechanical Engineering and Technology 56 ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 NOMENCLATURE b Semichord length rα Radius of gyration about elastic axis c Airfoil chord length xα Static unbalance in the structural model Ch Damping coefficients in heave Vf Flutter speed index Cα Damping coefficients in pitch u, v Cartesian velocity components Cd Drag coefficient p Pressure Cd Drag coefficient t Time Cl Lift coefficient ρ Density Cm Moment coefficient τ Non-dimensional time Cp Pressure coefficient E Total energy per unit volume DOF Degree of freedom k Reduced frequency based on semichord length, ωf b/ U∞ h Airfoil vertical displacement Φ, Ψ Torsional and bending mode shapes Iα Moment of inertia α Angle of attack J Jacobian of transformation η Non-dimensional span length of wing l Wing width ξ,η Transformed coordinates L Lift force ωf Frequency of vibration M∞ Freestream Mach number ωh Natural frequency in plunging ME.A Pitching moment around elastic axis ωα Natural frequency in pitching U∞ Freestream velocity μ Airfoil-to-air-mass ratio E.A Elastic axis τxx,yy,xy Shear stresses C.G Center of gravity REFERENCES Ghadiri, B., & Razi, M. (2007). Limit cycle oscillations of rectangular cantilever wings containing cubic nonlinearity in an incompressible flow. Journal of Fluids and Structures, 23, 665-680. Fung, Y. C. (1945). An introduction to the theory of Aeroelasticity. NewYork: Dover Publication Inc. Jiunn-Chi, W., Kaza K. R. V. & Sankar, L. N. (1987). Technique for the prediction of airfoil flutter characteristics in separated flow. Journal of Aircraft, 26(2), 168- 177. Halfman, R. L., Johnson, H. C. & Haley, S. M. (1951). Evaluation of high angle-of- attack aerodynamic derivative data and stall-futter prediction techniques. NACA TN 2533. Rainey, A. G. (1958). Measurement of aerodynamic forces for various mean angles of attack on an airfoil oscillating in pitch and on two finite span wings oscillating in bending with emphasis on damping in the stall. NACA Report 1305. Stall Flutter Control of a Wing Section by Leading Edge Modifications ISSN: 2180-1053 Vol. 8 No.1 January – June 2016 57 Price, S. J., Keleris & J. P. (1955). Non-linear dynamics of an airfoil forced to oscillate in dynamic stall. Journal of Sound and Vibration, 193, 823-864. Razak, N. A., Andrianne, T. & Dimitriadis, G. (2011). Flutter and stall flutter of rectangular wing in a wind tunnel. AIAA Journal, 49(10), 2258–2271. Song, Sh., Qinghua, Z., Chenglin, Zh. & Xianping, N. (2011). Airfoil aeroelastic flutter analysis based on modified leishman-beddoes model at low Mach number. Chinese Journal of Aeronautics, 24, 550-557. Bhat, S. S., Govardhan & R. N. (2013). Stall flutter of NACA 0012 airfoil at low Reynolds numbers. Journal of Fluids and Structures, 41, 166-174. Sun, Zh., Haghighat, S., Liu, H. & Bai, J. (2015). Time-domain modeling and control of a wing-section stall flutter. Journal of Sound and Vibration, 340, 221–238. Ning, Y., Nan, W., Xin, Zh. & Wei, L. (2016). Nonlinear flutter wind tunnel test and numerical analysis of folding fins with freeplay nonlinearities. Chinese Journal of Aeronautics, 29(1), 144-159. Blazek, J. (2001). Computational Fluid Dynamics: Principles and Applications. Amsterdam, London, New York: Elsevier. Hodges, D. H. & Pierce, G. A. (2002). Introduction to Structural Dynamics and Aeroelasticity. Cambridge University Press. Marzocca, P., Librescu, L. & Silva, W. A. (2002). Aeroelastic response and flutter of swept aircraft wings. AIAA Journal, 40, 801–812. Charles, L. L., Acquilla. S. H. & William, G. J. (1987). Pressure distributions from high Reynolds number transonic tests of an NACA 0012 airfoil in the langley 0.3- meter transonic cryogenic tunnel. NASA Langley Research Center 15. McAlister, K., Carr, L. W. & McCorosky, W. J. (1982). Experimental study of dynamic stall on advanced airfoils. NASA TM, 84245-TR, 82-A-8. Wilcox, D. C. (1994). Turbulence Modeling for CFD. California: DCW industries Inc. Kirshman, D. J. & Liu, F. (2006). Flutter prediction by an Euler method on non-moving Cartesian grids with gridless boundary conditions. Computers & Fluids, 35, 571–586. Carr, L. W., McAlister, K. W. & McCorosky, L. J. (1977). Analysis of the development of dynamic stall based on oscillating airfoils. NASA TN, D-8382.