On limit cycle oscillations of typical aeroelastic section with different preset angles of incidence at low airspeeds

Journal of Fluids and Structures 74 (2017) 19–34

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

On limit cycle oscillations of typical aeroelastic section with different preset angles of incidence at low airspeeds Carlos R. dos Santos a, *,1 , Daniel A. Pereira b,2 , Flávio D. Marques a,3 a b

São Carlos School of Engineering, University of São Paulo (EESC/USP), São Carlos, SP, Brazil Aeronautics Institute of Technology (ITA), São José dos Campos, SP, Brazil

highlights • • • •

Bifurcation regions with the dominance of structural and aerodynamics nonlinearities were distinguished. Jumps from damped oscillations to LCOs were detected. The computational model have successfully characterized the dynamic stall loads during LCOs. The flutter onset velocity increases for preset angles from zero up to ten degrees and then reduces quickly.

article

info

Article history: Received 20 April 2017 Received in revised form 26 June 2017 Accepted 7 July 2017

Keywords: Aeroelasticity Dynamic stall Flutter LCO Stall-induced oscillations Hopf bifurcation Semi-empirical model

*

a b s t r a c t Helicopter blades and wind turbines are examples of aeroelastic systems that can reach high angles of attack and can vibrate due to the effects of the dynamic stall, thereby leading to fatigue problems or performance loss. Structural and aerodynamic nonlinearities influence the aforementioned behavior and their modeling is crucial for phenomena characterization. Such system modeling requires proper knowledge of the physical events during the stall, which can be better achieved by validating the model with experimental data. This work investigates the nonlinear dynamics of a NACA 0012 airfoil under the influence of structural and aerodynamic nonlinearities due to dynamic stall effects at high angles of attack. Experimental and numerical analyses are carried out. Moreover, different preset incidence angles for the typical aeroelastic section are also considered. The aeroelastic signals are used for estimating the Hopf bifurcation onset and to build the bifurcation diagrams. By using a typical section model with two degrees of freedom coupled to the Beddoes–Leishman aerodynamic model, numerical results have been able to capture with good precision experimental features. The onset of the Hopf bifurcations allows the determination of the flutter critical airspeed. Results for zero preset angle show that limit cycle oscillations from small to moderate displacements are mostly driven by the hardening nonlinearity. After reaching larger angles of incidence the dynamic stall nonlinearities supplant those from structural sources. For higher preset angles, the dynamic stall effects tend to increase the energy associated with pitching motion and to reduce amplitudes in plunge motion. Another effect related to the aerodynamic nonlinearities relies on the increase of the flutter velocity by around 10% for preset angles ranging from zero up to

Corresponding author. E-mail addresses: [email protected] (C.R. dos Santos), [email protected] (D.A. Pereira), [email protected] (F.D. Marques).

1 M.Sc. degree candidate. 2 Ph.D. candidate. 3 Associate Professor. http://dx.doi.org/10.1016/j.jfluidstructs.2017.07.008 0889-9746/© 2017 Elsevier Ltd. All rights reserved.

20

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

ten degrees. For higher preset angles an abrupt reduction in the flutter onset velocity is observed. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Nonlinear aeroelastic phenomena can arise from structural or aerodynamic sources. Limit cycle oscillations (LCO) initiating at the Hopf bifurcation point have been often observed in aeroelastic systems (Dowell et al., 2003). Nonlinear structural behavior is usually related to concentrated or distributed influences resulting from large geometric displacements or material properties. Concentrated structural nonlinearities have been widely used in aeroelasticity research in the form of hardening and softening stiffness effects (Pereira et al., 2016). Either viscous or compressibility effects influence nonlinear aerodynamic loading. For instance, flow separation during pre- and post-stall is driven by viscous behavior effects, while shock excursion in unsteady transonic regime leads to typical nonlinear behavior due to compressibility effects (Dowell and Tang, 2002). The design of rotary wings admitting combined effects of structural and aerodynamic nonlinearities represents a formidable challenge, including the prediction of regimes of operation subject to dynamic stall effects. From the authors’ best knowledge, it is still important to investigate simplified concepts, such as typical airfoil sections and wings. Aeroelastic analysis of combined structural and aerodynamic nonlinearities for finite wings has been explored by Tang and Dowell (2002), where a highly flexible, long aspect ratio wing with a tip concentrated mass was tested. Supercritical and subcritical LCOs were found over a range that involved the stall regime. Stanford and Beran (2013) have used the same slender wing to explore a direct formulation to the Hopf bifurcation point and to construct the bifurcation diagram without the need for extensive time integration of aeroelastic equations of motion. Zakaria et al. (2015) have observed the influence of preset angles on the energy harvesting from an oscillating composite beam. When exposed to a flow field, the beam reaches self-excited motion with large geometric deformations associated with a nonlinear aerodynamic loading. Regarding nonlinear aeroelastic investigations, Camilo et al. (2013) presented a numerical study on the behavior of a typical section by considering structural nonlinearities and transonic flow regime. At low airspeeds, an important source of aerodynamic nonlinearity is associated with the flow separation at high angles of incidence. Dimitriadis and Li (2009), Razak et al. (2011), and Vasconcellos et al. (2016) have shown experimental results where LCOs driven by dynamic stall effects were observed. Such nonlinear phenomenon leads to responses that are related to complex trailing and leading edge vortex interactions (McCroskey, 1982; Carr, 1988; Mulleners and Raffel, 2013; Choudhry et al., 2014). This paper presents an experimental and numerical investigation on structural and aerodynamic nonlinear effects on a two-degree-of-freedom typical aeroelastic section for different preset angles of incidence. The preset angle of incidence must be understood as a fixed incidence to the airfoil from which the aeroelastic oscillations occur. The aerodynamic nonlinear behavior is expected when higher angles of attack are attained, while the structural nonlinearity comes from the hardening effect present in the pitch spring stiffness. The experimental procedure is carried out with a rigid wing in a suspension prepared to ensure pitch and plunge motions. Adjustments to the airfoil incidence guarantee a range of preset angles. Aeroelastic signals are acquired and used for validating a numerical model based on the typical section equations of motion and unsteady aerodynamics computed through a semi-empirical approach, that is, the Beddoes–Leishman dynamic stall model. Limit cycle oscillations are observed by exposing the airfoil to the flow field, and aeroelastic signals are collected from this motion for subsequent analysis. Bifurcation diagrams are achieved from aeroelastic time histories, regarding the airspeed evolution for different preset angles of incidence. The results demonstrate that the effect of stall-induced oscillations can surpass the effect of inherent structural nonlinearities when the angle of attack is higher than the stall angle of attack for the NACA 0012, by changing considerably the behavior of the system and the flutter onset velocity. Simulations using the numerical model allow the assessment of bifurcation onset points and inferences on the physical aspects that occur in the aeroelastic system. 2. Experimental apparatus The experiments were performed in a blower-type wind tunnel (open jet) with an output section of 0.5 × 0.5 m and able to reach approximately 23 m/s maximum airspeed. The experimental apparatus comprises a rigid wing with a NACA 0012 airfoil profile with a span length of 0.6 m supported by an elastic suspension. The device is supposed to be correlated with the typical aeroelastic section considering pitch and plunge motions. Fig. 1 sketches the wing and its elastic suspension. While the plunge motion is restrained by four elastic steel beams (two on each side at the top and at the bottom), the pitching motion is restrained by a piano wire going through the elastic axis. Nonlinear structural behaviors coming from pitch and plunge elastic features were characterized, and are presented further in this text. Structural damping effects estimations were also included in the dynamical model. A preset angle of incidence is defined as the structural equilibrium position of the airfoil in wind-off condition, including non-zero incidences. Therefore, aeroelastic oscillations will occur about this angle. Fig. 2 depicts the plate used for adjusting the preset angle, thus enabling a fixed initial condition in incidence to be set, which is responsible for enhancing the effect of aerodynamic nonlinearities.

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

21

Fig. 1. Sketch of the experimental apparatus.

Fig. 2. Details of the preset adjustment device.

Linear and angular encoders (US Digital optical encoders) were used to acquire aeroelastic information from the wind tunnel testing for the plunge and pitching motions, respectively. Data acquisition has been filtered by a low pass filter, cutting R out frequencies higher than 50 Hz, and processed with dSPACE⃝ (DS1006 Processor and Controller Board) connected to a PC computer to control the data assessment. Wind tunnel control was done manually by the user and true airspeed was calculated from static pressure and temperature readings during the tests.

3. Aeroelastic equations of motion The airfoil aeroelastic model is illustrated in Fig. 3, where pitch and plunge concentrated stiffnesses represented by springs are fixed to the airfoil of chord length c at the so-called elastic axis. The pitch spring is supposed to behave in a nonlinear manner. The plunge motion is assumed positive downwards and a preset angle θp is added to the structural pitch angle θ¯ (t) to yield the airfoil incidence angle θ (t) = θ¯ (t) + θp . Unsteady aerodynamic force and moment are applied at the elastic axis.

22

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

Fig. 3. Typical aeroelastic section.

The aeroelastic equations of motion are given by:

[ µe

xθ rθ2

xθ

[ 2 ωh

+

0

] {¨ } [ ] {˙ } ¯ ¯ d1,1 d1,2 h(t) h(t) ¨θ¯ (t) + d2,1 d2,2 θ˙¯ (t) ]{ { } } ¯ 4U 2 0 −Cl (t) h(t) = , θ¯ (t) rθ 2 ωθ 2 F (θ¯ ) µπ c 2 2Cmea (t)

(1)

¯ = where h(t) is the dimensionless plunge displacement, µe is the ratio between the total apparatus mass, mT , and the c airfoil mass, ma , µ = 4ma /(πρ c 2 ) is the airfoil mass ratio, xθ is the dimensionless distance from the elastic axis to the airfoil center of gravity (CG), rθ is the radius of gyration about elastic axis, ωh is the uncoupled natural frequency in plunge, ωθ is the uncoupled natural frequency in pitch, U is the airspeed, F (θ¯ ) is a function to represent the nonlinearity applied to the pitch stiffness value, di,j are damping constants determined by the viscous damping Rayleigh model, Cl (t) = (2L(t))/(ρ U 2 c) is the aerodynamic lift force coefficient and Cmea (t) = (2Mα (t))/(ρ U 2 c 2 ) is the aerodynamic moment coefficient measured at the elastic axis, and ρ is the air density. The unsteady aerodynamic loading is computed through the Beddoes–Leishman method (Leishman and Beddoes, 1989). The Beddoes–Leishman method includes corrections to account for the effects of dynamic stall of pitching airfoils at high angles of attack. It is important to distinguish the incidence angle θ (t) (which is related to structural displacements) from the airfoil aerodynamic angle of attack α (t) (which is related to the aerodynamic loads) as the latter is given by: 2h(t)

α (t) = θ (t) + arctan

[˙

h(t)

]

U

= θ¯ (t) + θp + arctan

[˙

h(t) U

]

.

(2)

A dimensionless pitch rate definition is convenient to the aerodynamic model. If this rate is evaluated at the same point where the elastic axis is located, it can be expressed as: q(t) = θ˙ (t)

c U

.

(3)

The Beddoes–Leishman method can be conveniently written in a state space representation. This formulation adopts eight states to represent the unsteady potential (linear) aerodynamics of an airfoil with the fully attached flow consideration. The methodology generalizes the airfoil indicial response as function of the Mach number (Leishman and Nguyen, 1990) and p p the potential normal force coefficient, Cn , and aerodynamic pitching moment at quarter of the chord, Cm , are given in state space as,

{˙x} = [A] {x} + [B]

{

Cnp

p Cm

{ } α q

} = [C ] {x} + [D]

,

{ } α q

(4)

,

(5)

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

23

where {x} is the state vector of 8 aerodynamic states and [A] is the state matrix given by:

[A] = diag[a11 a22 a33 a44 a55 a66 a77 a88 ],

(6)

with coefficients:

( ) ⎧ 2U 2 ⎪ ⎪ a11 = −b1 β ⎪ ⎪ c ⎪ ⎪ ( ) ⎪ ⎪ ⎪ 2U ⎪ 2 ⎪ ⎨a22 = −b2 β c

⎧ 1 ⎪ ⎪ a55 = − ⎪ ⎪ b 3 Kα M T I ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎪ ⎨a66 = − b K T 4 αM I ( ) 2U ⎪ 2 ⎪ a77 = −b5 β ⎪ ⎪ ⎪ c ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎩a88 = −

,

⎪ 1 ⎪ ⎪ a33 = − ⎪ ⎪ Kα T I ⎪ ⎪ ⎪ ⎪ ⎪ 1 ⎪ ⎩a44 = −

,

(7)

KqM TI

Kq T I

√

where β = 1 − M 2 , M is the Mach number, Kα = 0.75[(1 − M) + πβ 2 M 2 (A1 b1 + A2 b2 )]−1 , Kq = 0.75[(1 − M) + 2πβ 2 M 2 (A1 b1 + A2 b2 )]−1 , Kα M = (A3 b4 + A4 b3 )[b3 b4 (1 − M)]−1 , KqM = 7[15(1 − M) + 3πβ M 2 b5 ]−1 , TI = c /Vs , with Vs denoting the sound speed, and the remaining parameters are: A1 = 0.3, A2 = 0.7, A3 = 1.5, A4 = −0.5, b1 = 0.14, b2 = 0.53, b3 = 0.25, b4 = 0.1, and b5 = 0.5. The input matrix [B] (cf. Eq. (4)) is:

[

1 0.5

[B] =

1 0.5

1 0

0 1

1 0

1 0

0 1

]T

0 1

,

(8)

and the output matrix [C ] (cf. Eq. (5)) by:

[

c [C ] = 11 c21

c12 c22

c13 0

c14 0

0 c25

0 c26

0 c27

]

0 , c28

(9)

with

( ) ⎧ 2U 2 ⎪ ⎪ A1 b1 c11 = Cnα β ⎪ ⎪ c ⎪ ⎪ ( ) ⎪ ⎪ ⎪ 2U ⎪ ⎪ c12 = Cnα β 2 A2 b2 ⎪ ⎪ c ⎪ ⎨ ) ( 4 −1 c = ⎪ 13 ⎪ ⎪ M Kα T I ⎪ ⎪ ) ( ⎪ ⎪ ⎪ −1 1 ⎪ ⎪ c14 = ⎪ ⎪ M Kq T I ⎪ ⎪ ⎩ c21 = c11 (0.25 − xac )

,

⎧ c22 = c12 (0.25 − xac ) ⎪ ⎪ ) ( ⎪ ⎪ −A3 −1 ⎪ ⎪ ⎪ c = 25 ⎪ ⎪ M b 3 Kα M T I ⎪ ⎪ ) ( ⎪ ⎪ ⎪ − 1 −A4 ⎨c = 26 M b 4 Kα M T I ⎪ ( ) ⎪ ⎪ C 2U ⎪ n α 2 ⎪ ⎪c27 = − b5 β ⎪ ⎪ 16 c ⎪ ⎪ ( ) ⎪ ⎪ ⎪ −7 −1 ⎪ ⎩c28 = 12M

,

(10)

KqM TI

where Cnα is the normal force coefficient curve slope with angle of attack in the linear range and xac is the position of the aerodynamic center from the leading edge of the airfoil as fraction of the chord length. The matrix [D] in Eq. (5) is given by:

⎡ 4 ⎢ [D] = ⎣ M 1 −

1

−

M

⎤

M ⎥ ⎦. 7

(11)

12M

An effective angle of attack due to circulatory terms, αE , can be defined as:

αE (t) = β 2

(

2U c

)

(A1 b1 x1 + A2 b2 x2 ) .

(12)

To impose nonlinear characteristics due to dynamic stall effects, more dynamical states must be included in the model. These new states depend on dimensionless time constants. However, they can be converted to a dimensional time approach by multiplying each one by the factor c /(2U). By doing so, the state x9 imposes a lag in the previously calculated normal force by considering a time constant Tp , that is:

( x˙ 9 =

2U c

)

p

Cn (t) − x9 TP

,

(13)

24

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

leading to the delayed normal force, Cn′ : Cn′ (t) = x9 .

(14)

Two new states, x10 and x11 , are related to the flow separation from the airfoil trailing edge. Leishman and Beddoes (1989) defined an empirical approach based on the Kirchhoff flow model to represent the portion f of the airfoil with respect to the chord length, where the flow is attached under static conditions. Then, f as function of an arbitrary angle αˆ is given by:

⎧ ⎨

f (αˆ ) =

1 − 0.3e

|αˆ |−α1

if ⏐αˆ ⏐ ≤ α1

⏐ ⏐

s1

⎩0.04 − 0.66e α1s−2|αˆ |

(15)

if ⏐αˆ ⏐ > α1 ,

⏐ ⏐

where α1 is defined when f (α1 ) = 0.7, and s1 and s2 are constants determined from experimental data. Under dynamic conditions, the delay observed in the leading edge separation is modeled as dependent on a time constant Tf0 . Then, the state x10 , which is related to an equivalent angle of attack associated with Cn′ , is

( x˙ 10 =

2U

) f ( Cn′ (t) ) − x Cn α

c

10

,

Tf

(16)

fa (t) = x10 ,

(17)

where the time constant Tf varies ( )from Tf0 during the vortex shedding phase and fa (t) is the value of the delayed separation point considering the angle

Cn′ (t) Cnα

.

For the state x11 , which is related to the separation point by considering the instantaneous angle of attack, it follows:

( x˙ 11 =

2U

)

c

2(f (α ) − x11 )

,

Tf0

(18)

fb (t) = x11 ,

(19)

with fb (t) being similar to fa (t). Taking fˆ (t) as the maximum value between fa (t) and fb (t), the aerodynamic coefficients contributions considering the trailing edge flow separation are given by:

⎧ √ ( )2 1 + fa (t) ⎪ f ⎪ ⎪ C (t) = C αE (t) nα ⎨ n 2 { } f ⎪ Cm (t) = K0 + K1 (1 − fˆ (t)) + K2 sin(π fˆ (t)m ) Cnα αE (t) + Cm0 ⎪ ⎪ √ ⎩ f Ct (t) = η Cnα fa (t) αE2 (t),

(20)

where typically m = 2, η represents the percentage √ of tangential force measured experimentally under fully attached conditions in relation to the value expected from Cnα fa (t)αE2 (t), K0 , K1 , and K2 are empirical constants that match the pitching moment coefficient curve under static conditions, and Cm0 is the pitching moment coefficient when the normal force is zero. The vortex shedding phase influences the airfoil unsteady aerodynamic behavior. Some methodologies have been proposed to represent the dynamic stall onset. For low wind speeds Sheng et al. (2006) and (2007) have used a lagged value of the angle of attack to trigger the beginning of the dynamic stall onset. Such flow event occurs when the value of the delayed angle of attack reaches a critical value. For the present work, the classical theory is based on the delayed value of the normal force coefficient. When |Cn′ | ≥ Cn1 , with Cn1 being the value of Cn at the stall imminence under static conditions, the vortex shedding phase is triggered. A dimensionless counter τv is established and marches according to the dimensionless time 2Ut /c. During this condition, the parcel of the normal force coefficient associated with the vortex shedding Cnv is given by the state x12 as:

( x˙ 12 =

2U

) ˙ Cv (t) − x12

c

Tv

,

(21)

Cnv (t) = x12 ,

(22)

where:

{ Cv =

Cnα [1 − 0.25(1 + 0

√

fa (t))2 ]αE (t)

, if τv ≤ 2Tv l , , if τv > 2Tv l ,

(23)

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

25

with Tv l being an empirical time constant defined as the dimensionless time window in which the shedding vortex travels from the leading edge to the trailing edge. The vortex travel over the airfoil also changes its pressure center position and a term of pitching moment coefficient due to this motion can be expressed as:

[

v (t) = −0.25 1 − cos Cm

(

πτv

)]

Tv l

Cnv .

(24)

Moreover, during the vortex shedding phase, the time constants Tf and Tv , and the angle α1 (cf. Eq. (15)) are modified. Chantharasenawong (2007) models these variations as follow:

Tf =

⎧ Tf0 ⎪ ⎪ ⎪ ⎪ ⎨ 1 Tf0 3 ⎪1

⎪ Tf ⎪ ⎪ ⎩2 0 4Tf0

⎧ T ⎪ ⎨ v0 0.25Tv0 Tv = ⎪0.5Tv0 ⎩ 0.9Tv0 α1 =

{

if 0 ≤ τv ≤ Tv l and α α˙ ≥ 0 if Tv l < τv ≤ 2Tv l and α α˙ ≥ 0 if 0 ≤ τv ≤ 2Tv l and α α˙ < 0

(25)

if 2Tv l < τv , if 0 ≤ τv ≤ Tv l and α α˙ ≥ 0 if Tv l < τv ≤ 2Tv l and α α˙ ≥ 0 if 0 ≤ τv ≤ 2Tv l and α α˙ < 0 if 2Tv l < τv ,

α10 α10 − [1 − fa (t)]0.25 δα1

if α α˙ ≥ 0 if α α˙ < 0,

(26)

(27)

with the subscript 0 denoting the value of each parameter before the vortex shedding and δα1 the maximum change in α1 during this phase. After the vortex shedding phase, the reattachment phase is triggered when |Cn′ (t)| < Cn1 and |α (t)| is decreasing. Moreover, it is also necessary to take α1 = α10 , Tv = Tv0 , and a better representation to Tf variation is given by:

{ Tf =

Tf0 2Tf0

if fa (t) ≥ 0.7 if fa (t) < 0.7.

(28)

The total aerodynamic coefficients can be summarized as:

⎧ ⎨Cn (t) = Cnp (t) + Cnf (t) + Cnv (t) − Cnα αE (t) p f v Cm (t) = Cm (t) + Cm (t) + Cm (t) ⎩ f Ct (t) = Ct (t).

(29)

By decomposing the normal and tangential forces coefficients, the lift coefficient is: Cl (t) = Cn (t) cos[α (t)] − Ct (t) sin[α (t)],

(30)

and the pitching moment coefficient at the elastic axis results: Cmea (t) = Cn (t)(xea − xca ) + Cm (t),

(31)

where xea is the elastic axis position from the airfoil leading edge in fraction of the chord length. 4. Results and discussions The characterization of the wing model suspension under wind-off condition is carried out to assess the structuraldynamics parameters. Fig. 4 depicts the plunge and pitch springs restoring force and moment with the respective displacements. Fig. 4(a) shows the measured spring force curve for the plunging motion and linear extrapolation for those measured values was assumed. For the restoring pitching moment, Fig. 4(b) can be used to observe the hardening effect associated with the spring stiffness. Polynomial fitting can be used to express a cubic pitching restoring moment, that is: 30.7 θ¯ 3 + 6.96 θ¯ (Nm). Linear extrapolation was also included in Fig. 4(b) to compare with the cubic extrapolation that represents the hardening effect at higher angles of incidence, as well as to estimate the uncoupled pitch natural frequency. After certifying that the experimental model presents linear plunge spring and hardening effect in the pitch spring, the remaining parameters required by the numerical model can be seen in Table 1. Aeroelastic experimental tests were performed by establishing a wind tunnel airspeed and disturbing the airfoil with an input in either plunge or pitch displacement. Experimental tests were performed for five preset angle of incidence: 0◦ , 11.5◦ , 12.5◦ , 17.5◦ , and 19.5◦ . These preset values were taken based on observing the aeroelastic responses for which stall-induced oscillations were clearly present.

26

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

(a) Plunging motion.

(b) Pitching motion.

Fig. 4. Plunge and pitch spring stiffness values determination from the respective restoring force and moment curves.

Table 1 Experimental model constants. Variables

Description

Values

c

xθ xea rθ

Chord length (m) Air density (kg/m3 ) Airfoil mass (kg) Total mass (kg) Decoupled plunge natural frequency (rad/s) Decoupled pitch natural frequency (rad/s) Dimensionless distance between elastic axis and CG of airfoil Dimensionless elastic axis position from the leading edge Dimensionless rotational inertia term about elastic axis

0.25 1.079 1.5 3.6733 21.77 24.85 0.66 0.25 0.7303 [

di ,j

Damping matrix elements using Rayleigh approach

ρ ma mT

ωh ωθ

5.49 10.03

(n)

Linear flutter velocity (numeric) (m/s)

13.71

(e)

Linear flutter velocity (experimental) (m/s)

14.2

Uf

Uf

10.03 32.48

]

Fig. 5 shows the aeroelastic time series of experimental results for zero preset angle of incidence. Figs. 5(a) and 5(b) present the responses for the wind-off condition, and the characteristic coupled damped signals in plunge and pitch. As already expected, at certain wind speed the flutter occurs. Figs. 5(c) and 5(d) show the aeroelastic responses for flutter velocity, which has been measured at 14.2 m/s. The plot reveals that small flutter amplitudes have been observed for zero preset angle. Consequently, linear models can be considered fairly reliable with the classical pitch–plunge modal coupling. Fig. 6 depicts the power spectra evolution of the experimental signals acquired at different freestream velocities. It is worth observing that the frequencies measured for pitch and plunge in wind-off condition are not equal to their uncoupled natural frequencies, but they reflect the response of a coupled aeroelastic system. From this result, it is clear that the linear flutter mechanism with the plunge and pitch modes coalescence. When the flow velocity is increased, after 14.2 m/s the pitch/plunge uncoupled natural frequencies become the same equal to 3.38 Hz which is called the fundamental frequency in the flutter regime. After crossing the critical flutter boundary, which leads to increasing aeroelastic displacements, nonlinear behavior arises. Then, at moderately large displacements in pitch, the hardening and mild nonlinear aerodynamic effects drive the responses to LCOs. Fig. 7 shows the power spectra of the plunge and pitch motions when the onset flutter point has been reached during experiments. One can note from the spectra that the typical LCO frequency content compound of a fundamental harmonic (3.38 Hz) and its even and odd super-harmonics. Experimental data have shown the transitional dynamics from mostly linear behavior to LCO at flutter regime. If the airspeed increases from the bifurcation point the LCO amplitudes should also increase. Further tests have been made for incremented values of airspeed. Fig. 8 presents the bifurcation diagrams of measured plunge and pitch LCO amplitudes from the flutter observation at 14.2 m/s at zero preset angle condition. Here, the numerical model was also used and the results of the simulated bifurcation diagram were included in the respective pitch (cf. Fig. 8(a)) and plunge (cf. Fig. 8(b)) experimental results. The computational model is simulated using the parameters from Tables 1 and 2. The parameters required for applying the Beddoes–Leishman model were taken from Leishman and Beddoes (1986). Adequate agreement between the measured amplitudes and those predicted by the complete computational model is attained regardless of the

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

(a) Plunge response at wind-off.

(b) Pitch response at wind-off.

(c) Plunge response at 14.2 m/s.

(d) Pitch response at 14.2 m/s.

27

Fig. 5. Time series for the plunging and pitching motion when the preset angle is zero degrees.

use of a classical formulation for the aerodynamic model. A closer observation of the simulated bifurcation diagram reveals that the flutter onset happens at smaller airspeed when compared to the experimental observation, that is, at 13.71 m/s. The simulated LCOs grow smoothly with airspeed increments. Moreover, the computational model indicates that near 16 m/s the LCO amplitudes change their pattern as airspeeds increase. After that point, pitch LCOs reduce in amplitude growth rate (cf. Fig. 8(a)), while for the plunge LCOs a slight reduction in amplitudes is observed (cf. Fig. 8(b)). It is worth noting that when the LCO amplitudes reach angles of attack near and larger than the static stall values, dynamic stall effects are more pronounced. Therefore, one may infer that the changes in the LCOs evolution with airspeed observed from simulated results are related to the stall-induced oscillations effects. When the same sort of simulations are performed with a linear aerodynamics (using only the first eight states of the Beddoes–Leishman method), but keeping nonlinear pitch stiffness, the LCOs amplitudes for velocities below 16 m/s are kept like those simulated with the full Beddoes–Leishman aerodynamic model. Discrepancies between the results from the numerical models with different aerodynamics approaches are observed at higher incidence angles near the stall. To ensure that linear aerodynamics fail to predict stall-induced oscillations, the same comparison was performed for the preset angle of 11.5◦ and the results are illustrated in Fig. 9. It is observed that a model based on linear aerodynamics does not accurately predict the aeroelastic behavior. Although LCOs are present only due to nonlinear structural effects, linear aerodynamics does not provide coherent bifurcation onset and amplitudes when oscillating at higher angles of attack near and over the stall condition. Fig. 10 illustrates how LCO frequencies vary with the airspeed for the case of zero preset angle. Experimental and computational results are compared considering the complete nonlinear model. It is possible to observe that the computational model tends to predict oscillations with lower frequencies than the frequencies measured experimentally and with higher magnitudes. It can be associated with inaccuracies coming from the characterization of structural features, possibly related to uncertainties in the damping matrix.

28

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

7 6 5 4 3 2 14

1 12 10

0 10

8 8

6 6

4

4

2

2

0

0

Fig. 6. Study of the experimental relation between the pitch–plunge coupling frequencies responsible for leading to the instability of the system.

102

10-6 100

10-8 10-2

10-10 10-4

10-12 10-6 0

2

4

6

8

10

12

14

16

18

20

22

24

0

2

(a) Plunge power spectra.

4

6

8

10

12

14

16

18

20

22

24

(b) Pitching power spectra.

Fig. 7. LCO frequency content of the plunge and pitch motions at the observed flutter onset (U = 14.2 m/s) at zero preset angle. Table 2 Empirical parameters for the Beddoes–Leishman method (Leishman and Beddoes, 1986). Cnα (rad−1 )

α10

s1 (deg)

s2 (deg)

K0

K1

K2

η

Cn1

Tp

Tf 0

T v0

Tv l

δα1

(deg)

6.188

15.25

3.00

2.30

0.0025

−0.135

0.04

0.965

1.45

1.7

3.0

6.0

7.0

2.10

(deg)

Preset angles of incidence can be used to take the system to particular incidence motion regions, thereby allowing the examination of aerodynamic nonlinearities related to dynamic stall effects. Figs. 11 and 12 show the bifurcation diagrams for pitch and plunge experimental results, respectively, considering all preset angles cases. In all these cases, LCO amplitudes reach high angles of attack above the static stall angle of the NACA 0012 airfoil. Experimental results were also compared to the numerical model simulations, which are depicted in Fig. 11 and 12 by continuous lines. A clear change in the dynamic behavior is revealed from each preset angle case. The closest points to experimentally acquired flutter onset are

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

(a) Pitch.

(b) Plunge.

Fig. 8. Bifurcation diagrams for the zero degrees of preset angle.

(a) Pitch.

(b) Plunge. Fig. 9. Bifurcation diagram for the 11.5◦ of preset angle.

Fig. 10. Frequency content of LCOs for the zero of preset angle.

29

30

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

Fig. 11. Bifurcation diagrams of pitching motion for all preset angles cases.

Fig. 12. Bifurcation diagrams of plunging motion for all preset angles cases.

presented in the bifurcation diagrams. It can be observed that the numerical model anticipates the flutter onsets for all preset cases different from the zero one. A reason for the prediction errors can be related to other structural nonlinearities and aerodynamic effects (asymmetries and other airfoil geometrical imperfections). Considering the good match between the LCO amplitudes for all preset angle cases, the predictive capabilities of the numerical code has been considered to be adequate. For the particular case of 19.5◦ preset angle, no flutter point has been depicted in Figs. 11 and 12, because experimental signals at airspeeds below the LCO onset have presented very irregular oscillations with small amplitudes. For preset angles of 11. 5◦ , 17.5◦ and 19.5◦ the numerical model predicts a sudden growth of LCO amplitudes, that is, a fast jump from damped oscillations to LCO with considerably large amplitudes. This fact suggests how difficult it is to capture the flutter onset velocity experimentally. In the case of the preset angle of 12.5◦ a smooth bifurcation shape is also observed. Here the peak-to-peak LCOs amplitudes have reached the largest values among the pitching angles. For the preset angles of 17.5◦ and 19.5◦ the numerical simulations also indicate that even at low airspeeds small amplitudes of oscillation are observed. This phenomenon can be related to unsteady transient loads due to flow separation at high angles of attack, above the static stall angle. The bifurcations diagrams for pitch and plunge motions have shown different behavior with respect to the increase of preset angles. While the pitch LCOs keep a typical range of peak-to-peak amplitude values, the plunge oscillations have much smaller amplitudes and its mean value of oscillation reduces as the freestream velocity increases. In fact, except for the zero degrees preset angle, all the plunge motion cases have a significant reduction in LCO amplitudes. Moreover, even in the zero degrees preset case for airspeeds above 16 m/s a reduction in plunge amplitudes is expected. It can be concluded that the dynamic-stall-related nonlinearities on the LCOs drives the flutter mechanism in such a way that the energy spent in pitch

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

31

Fig. 13. Analysis for two different preset conditions at U = 15.5 m/s.

motion is larger than that corresponding in the plunge motion. In order to investigate this hypothesis, the aerodynamic loading coefficients were observed. For the sake of comparison, the case with zero degrees preset angle and that with 11.5◦ at the same airspeed of 15.5 m/s were considered. Fig. 13 presents the aerodynamic lift and pitch moment – at the elastic axis – coefficients obtained computationally and the phase portraits for pitch and plunge motions in both situations. It is worth noting that while the airfoil with zero degrees preset angle presents characteristics of linear aerodynamics, the airfoil with a preset angle of 11.5◦ presents aerodynamical behavior of an airfoil already dominated by the dynamic stall with hysteretical lift and a distinguishable pitch down moment. The aerodynamic moment tends to elevate the rate of variation of the pitch angle transmitting energy to this motion. As the freestream velocity is the same in the two preset cases, a reduction in the energy of the plunge motion is expected. Such phenomenon is verified from both experimentally and numerically results by the reduction of plunge amplitudes and their rate of variation. Moreover, there are also consequences on the lift coefficient hysteresis that reaches higher values at higher angles of incidence when compared with the case with zero degrees of preset. Nonetheless, as the angle of incidence decreases the lift coefficient becomes lower than in the linear case. For the most severe preset angle cases, the motion features and loading behavior have also been verified. Fig. 14 presents how the aerodynamic loading coefficients vary for preset angles cases 12.5◦ , 17.5◦ , and 19.5◦ at the airspeed of 12 m/s. In these cases, the pitching down moment is stronger as the preset angle increases and the lift coefficient has an hysteresis progressively larger indicating strong deep stall behavior, especially for presets of 17.5◦ and 19.5◦ . As exposed during the formulation of the Beddoes–Leishman aerodynamic model, the methodology is influenced by physical events that occur during the dynamic stall. Fig. 15 indicates how two of these events vary with the airspeed and the preset angle during aeroelastic motions, that is: the angle where the vortex shedding phase starts (i.e., when |Cn′ (t)| is growing and becomes larger or equal to Cn1 ), and the angle where the flow reattachment phase starts (i.e., when |Cn′ (t)| is decreasing and becomes lower or equal to Cn1 ). If one of those conditions is not attained, the angles were not plotted for the respective airspeed in Fig. 15. Comparing the results in Fig. 15 to those in Fig. 11, it is possible to conclude that the incidence angle where the vortex shedding phase starts is almost the maximum incidence angle attained by the airfoil. This fact reinforces the importance of the large pitch down moment, which is consequence of the vortex shedding phase and quickly reduces the airfoil angle of attack during aeroelastic motions. Aerodynamic nonlinearities are not only related to the reduction of plunge amplitudes of motion but also explain the modification of the flutter onset velocity as a function of the preset incidence angle. This variation can be observed in Fig. 16

32

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

Fig. 14. Analysis for three different preset conditions at U = 12 m/s.

Fig. 15. Variation of the angle where the vortex shedding phase starts (continuous lines) and the reattachment phase starts (dashed lines).

for numerical solutions for the critical flutter speed. Here the simulations have been carried out assuming steps of 0.1 m/s to yield a flutter speed Uc . A normalization with respect to the zero degrees preset angle case critical flutter speed Uc0 was admitted. The computational model indicates that the flutter onset velocity increases in a rate of 10% as a function of the preset angle until approximately 10◦ . A similar result was obtained by Bichiou et al. (2014) when simulating aeroelastic airfoils with an unsteady lumped vortex aerodynamic model. Given a preset angle and before the flutter velocity, they observed that the equilibrium position of the airfoil increases as the freestream velocity increases. The same increasing effect

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

33

Fig. 16. Effect of the preset angles on the flutter onset.

on the equilibrium position of the airfoil is observed by increasing the preset angle at the same determined airspeed. Then, the increase in the critical flutter velocity is a consequence of the increase in the equilibrium position of the airfoil due to the preset angle. An analogous result was observed by Patil et al. (2000) where structural nonlinearities lead to the increase in the flutter velocity even if linear aerodynamics is considered. However, for preset angles higher than 10◦ the critical flutter speed presents severe reduction. This feature can be related to the fact that the airfoil is near to its stall condition as the preset angles become larger, thereby speeding up the process of stall induction leading to oscillations. 5. Concluding remarks The present paper has presented an investigation on airfoil pitch and plunge LCOs at low airspeeds and considering preset angles of incidence. Experiments were conducted in a blower-type wind tunnel and different preset angles were assumed to acquire LCO signals for a range of freestream velocities. Experimental results were compared with a numerical typical aeroelastic section model where the unsteady aerodynamic loading is provided by the Beddoes–Leishman method for dynamic stall corrections. Experimental apparatus characterization and numerical model validation were presented, where good agreement between experimental and computational results was observed. The computational model has indicated that for some preset angle conditions the amplitudes of motion do not have a smooth bifurcation onset, but instead a fast jump from damped oscillations to LCOs with considerably large amplitudes. It reflects on some difficulties on the experimental assessment of the flutter onset velocity for these cases. When studying the case with a preset angle of zero degrees the existence of two regions was noticed. One region presents dominance of structural nonlinearities on the LCOs and another indicates stall related aerodynamic nonlinearities with more influence on the motion. For higher angles of preset the aerodynamics are dominated by the dynamic stall and the deep stall. These situations have led to pitching motion absorbing more energy than the plunge motion, therefore with small amplitudes. Aerodynamic nonlinearities also influence the flutter onset velocity. For preset angles below 10 degrees the flutter speed increases up to 10% in relation to that calculated with a preset angle of zero degrees. However, after this preset angle the flutter velocity quickly reduces as an effect of the proximity with the static stall angle. Acknowledgments The authors acknowledge the sponsorship of the Brazilian agencies: grant #2016/24522-4, São Paulo Research Foundation (FAPESP) and grants #307658/2016-3 and #145439/2015-1, National Council for Scientific and Technological Development (CNPq). References Bichiou, Y., Nuhait, A.O., Abdelkefi, A., Hajj, M.R., 2014. Unsteady aeroelastic behaviors of rigid airfoils with preset angles of attack. J. Vib. Control 22 (4), 1010–1022. http://dx.doi.org/10.1177/1077546314537106. Camilo, E., Marques, F.D., Azevedo, J.L.F., 2013. Hopf bifurcation analysis of typical sections with structural nonlinearities in transonic flow. Aerosp. Sci. Technol. 30 (1), 163–174. http://dx.doi.org/10.1016/j.ast.2013.07.013.

34

C.R. dos Santos et al. / Journal of Fluids and Structures 74 (2017) 19–34

Carr, L.W., 1988. Progress in analysis and prediction of dynamic stall. J. Aircr. 25 (1), 6–17. http://dx.doi.org/10.2514/3.45534. Chantharasenawong, C., 2007. Nonlinear Aeroelastic Behaviour of Aerofoils Under Dynamic Stall (Ph.D. thesis), Imperial College London (University of London). Choudhry, A., Leknys, R., Arjomandi, M., Kelso, R., 2014. An insight into the dynamic stall lift characteristics. Exp. Therm Fluid Sci. (ISSN: 0894-1777) 58, 188–208. http://dx.doi.org/10.1016/j.expthermflusci.2014.07.006. Dimitriadis, G., Li, J., 2009. Bifurcation behavior of airfoil undergoing stall flutter oscillations in low-speed wind tunnel. AIAA J. 47 (11), 2577–2596. http://dx.doi.org/10.2514/1.39571. Dowell, E.H., Edwards, J., Strganac, T.W., 2003. Nonlinear aeroelasticity. J. Aircr. 40 (5), 857–874. http://dx.doi.org/10.2514/2.6876. Dowell, E.H., Tang, D., 2002. Nonlinear aeroelasticity and unsteady aerodynamics. AIAA J. 40 (9), 1187–1196. http://dx.doi.org/10.2514/2.1853. Leishman, J.G., Beddoes, T.S., 1986. A generalised model for airfoil unsteady aerodynamic behaviour and dynamic stall using the indicial method. In: Proceedings of the 42nd Annual Forum of the American Helicopter Society, Washington DC, pp. 243–265. Leishman, J.G., Beddoes, T.S., 1989. A semi-empirical model for dynamic stall. J. Amer. Helicopter Soc. 34, 3–17. http://dx.doi.org/10.4050/JAHS.34.3. Leishman, J.G., Nguyen, K.Q., 1990. State-space representation of unsteady airfoil behavior. AIAA J. 28 (5), 836–844. http://dx.doi.org/10.2514/3.25127. McCroskey, W.J., 1982. Unsteady airfoils. Annu. Rev. Fluid Mech. 14, 285–311. http://dx.doi.org/10.1146/annurev.fl.14.010182.001441. Mulleners, K., Raffel, M., 2013. Dynamic stall development. Exp. Fluids (ISSN: 0723-4864) 54 (2). http://dx.doi.org/10.1007/s00348-013-1469-7. Patil, M.J., Hodges, D.H., Cesnik, C.E.S., 2000. Nonlinear aeroelastic analysis of complete aircraft in subsonic flow. J. Aircr. 37 (5), 753–760. http://dx.doi.org/ 10.2514/2.2685. Pereira, D.A., Vasconcellos, R.M.G., Hajj, M.R., Marques, F.D., 2016. Effects of combined hardening and free-play nonlinearities on the response of a typical aeroelastic section. Aerosp. Sci. Technol. 50, 44–54. http://dx.doi.org/10.1016/j.ast.2015.12.022. Razak, N.A., Andrianne, T., Dimitriadis, G., 2011. Flutter and stall flutter of a rectangular wing in a wind tunnel. AIAA J. 49 (10), 2258–2271. http: //dx.doi.org/10.2514/1.J051041. Sheng, W., McD Galbraith, R.A., Coton, F.N., 2006. A new stall-onset criterion for low speed dynamic-stall. J. Solar Energy Eng. 128 (4), 461–471. http: //dx.doi.org/10.1115/1.2346703. Sheng, W., McD Galbraith, R.A., Coton, F.N., 2007. Improved dynamic-stall-onset criterion at low Mach numbers. J. Aircr. 44 (3), 1049–1052. http://dx.doi. org/10.2514/1.29163. Stanford, B., Beran, P., 2013. Direct flutter and limit cycle computations of highly flexible wings for efficient analysis and optimization. J. Fluids Struct. 36, 111–123. http://dx.doi.org/10.1016/j.jfluidstructs.2012.08.008. Tang, D., Dowell, E.H., 2002. Limit-cycle hysteresis response for a high-aspect-ratio wing model. J. Aircr. 39 (5), 885–888. http://dx.doi.org/10.2514/2.3009. Vasconcellos, R.M.G., Pereira, D.A., Marques, F.D., 2016. Characterization of nonlinear behavior of an airfoil under stall-induced pitching oscillations. J. Sound Vib. 372, 283–298. http://dx.doi.org/10.1016/j.jsv.2016.02.046. Zakaria, M.Y., Al-Haik, M.Y., Hajj, M.R., 2015. Experimental analysis of energy harvesting from self-induced flutter of a composite beam. Appl. Phys. Lett. 107 (2), 023901. http://dx.doi.org/10.1063/1.4926876.