Aeroelasticity

2. Hz. Active Aeroelasticityand Rotorcraft Lab., Seoul National University. Dynamic aeroelasticity. Then, aeroelasticsystems of equations becomes - Fo...

0 downloads 210 Views 2MB Size
Aeroelasticity

2014 Prof. SangJoon Shin

Active Aeroelasticity and Rotorcraft Lab.

Index 0. Introduction 1. Static Aeroelasticity 2. Unsteady Aerodynamics 3. Dynamic Aeroelasticity 4. Turbomachinery Aeroelasticity 5. Helicopter Aeroelasticity

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic Aeroelasticity

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity  Two principal phenomena -

Dynamic instability (flutter)

-

Responses to dynamic load, or modified by aeroelastic effects

 Flutter ∙∙∙ self-excited vibration of a structure arising from the interaction of aerodynamic elastic and internal loads “response” ∙∙∙ forced vibration “Energy source” ∙∙∙ flight vehicle speed

 Typical aircraft problems -

Flutter of wing

-

Flutter of control surface

-

Flutter of panel Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity  Stability concept If solution of dynamic system may be written or

= y ( x, t )

N

(σ k + iωk ) t ⋅ y x e ( ) ∑ k k =1

a) σ k < 0, ωk ≠ 0 ⇒ Convergent solution : “stable” b) σ= 0, ωk ≠ 0 ⇒ Simple harmonic oscillation : ”stability boundary” k c) σ k > 0, ωk ≠ 0 ⇒ Divergence oscillation : “unstable”

0 ⇒ Continuous convergence : “stable” d) σ k < 0, ωk = e) σ k= 0, ωk= 0 ⇒ Time independent solution : ”stability boundary” f) σ k > 0, ωk = 0 ⇒ Continuous divergence : “unstable”

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity  Flutter of a wing Typical section with 2 D.O.F

Kα , K h : torsional, bending

stiffness

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity  First step in flutter analysis -

Formulate eqns of motion

-

The vertical displacement at any point along the mean aerodynamic chord from the equilibrium z = 0 will be taken as za ( x, t )

za ( x, t ) =−h − ( x − xea )α -

The eqns of motion can be derived using Lagrange’s eqn

d  ∂L  dt  ∂qi

 ∂L = Qi −  ∂qi L= T − U

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity - The total kinetic energy(T)

1  ∂z  T = ∫ ρ  a  dx 2 − b  ∂t  2

b

b

2 1     ρ h + ( x − xea )α  dx = 2 −∫b  b

b

b

1 2 1 h ∫ ρ dx + hα ∫ ρ ( x − xea )dx + α 2 ∫ ( x − xea ) 2 dx = 2 −b 2 −b −b

m

(airfoil mass)



(static unbalance)



(mass moment of inertia about c.g.)

*Note) if xea = xcg , then Sα = 0 by the definition of c.g. Therefore, T=

1 2 1 2 mh + Iα + Sα hα 2 2

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity - The total potential energy (strain energy) U =

-

1 1 kh h 2 + kα α 2 2 2

Using Lagrange’s eqns with

L= T − U

= q1 h= α 1 , q2  mh + Sα α + kh h = Qh ⇒ Qα  Sα h + Iα α + kα α = Where Qh , Qα are generalized forces associated with two d.o.f’s h, α respectively.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity Qh = − L = − L(α , h, α , h, α, h, ⋅⋅⋅) Q M= M (α , h, α , h, α, h, ⋅⋅⋅) = α

ea

ea

Governing eqn. m ⇒  Sα

Sα   h   K h  + Iα  α   0

0   h   −L   =   Kα  α   M ea 

- For approximation, let’s use quasi-steady aerodynamics h L qSCLα (α + ) = U∞ M ac = qSc Cmα α M ea = ( xea − xac )  L + M ac

h ) + qSc Cmα α = eqSCLα (α + U∞

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity *Note) Three basic classifications of unsteadiness (linearized potential flow) i)

Quasi-steady aero: only circulatory terms due to the bound vorticity. Used for characteristic freq. below 2Hz (e.g., conventional dynamic stability analysis)

ii)

Quasi-unsteady aero: includes circulatory terms from both bound and wake vorticities. Satisfactory results for 2 Hz < ωα , ωh < 10 Hz . Theodorsen is one that falls into here. (without apparent mass terms)

iii)

Unsteady aero: “quasi-unsteady”+”apparent mass terms” (non-circulatory terms, inertial reactions:

α, h )

For ω > 10Hz , for conventional aircraft at subsonic speed. Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity

Then, aeroelastic systems of equations becomes  qSCLα  0    qSCLα   h  0   h   Kh  m Sα   h   U ∞   =  +   +    S   I α 0 qSeC K qSeC −   α α α     α Lα Lα  α    0  −qSc Cmα     − U∞   - For stability, we can obtain characteristic eqn. of the system and analyze the roots. neglect damping matrix for first, Sα   h   K h  +  Iα  α   0

  h  0   =   Kα − qSeCLα  α  0  Much insight can be obtained by looking at the undamped system m S  α

qSCLα

(Dowell, pp. 83)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity = α Set

α= e pt , h he pt

(mp 2 + K h ) ( Sα p 2 + qSCLα )   h  pt 0  ⇒   e =   2 2 + − ( ) S p I p K qSeC α   α α α Lα  0   For non-trivial solution, Characteristic eqn.,

det(∆) =0

(mIα − Sα ) p 4 + [ K h Iα + ( Kα − qSeCLα )m − qSCLα Sα ] p 2 + K h ( Kα − qSeCLα ) = 0 A

B

C

− B ± B 2 − 4 AC ∴p = 2A 2

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity The signs of A, B, C determine the nature of the solution.

A > 0, C > 0 (if q < qD) B Either (+) or (-) B = mKα + K h Iα − [me + Sα ]qSCLα •

If [me + Sα ] < 0, B > 0 for all q



Otherwise B < 0

when

Kα K h Iα  Sα  + − 1 + qSeCLα < 0  e me  me 

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity - Two possibilities for B ( B>0 and B<0) i)

B>0: ① B − 4 AC > 0, p are real, negative, so p is pure imaginary → 2

2

neutrally stable ② B 2 − 4 AC < 0, p 2 is complex, at least one value should have a positive real part → unstable ③ B 2 − 4 AC = 0 → stability boundary • Calculation of

qF

DqF2 + EqF + F = 0 ←(from B 2 − 4 AC = 0, stability boundary)

− E ± E 2 − 4 DF qF = 2D Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity where,

{

D ≡ [ me + Sα ] SCLα

{

}

2

}

E ≡ −2 [ me + Sα ][ mKα + K h Iα ] + 4  mIα − Sα2  eK h SCLα F ≡ [ mKα + K h Iα ] − 4  mIα − Sα2  K h Kα 2

① At least, one of the

qF

must be real and positive in order for flutter

to occur. ② If both are, the smaller is the more critical. ③ If neither are, flutter does not occur. ④ If Sα ≤ 0 (c.g. is ahead of e.a), no flutter occurs(mass balanced)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity ii) B<0: B will become (-) only for large q

0 will occur before B=0 since A>0,C>0 B 2 − 4 AC =

∴ To determine Examine

qF , only B>0 need to

be calculated.

p as q increases ±iω1 , ±iω2 ( B 2 − 4 AC > 0) Low q → p =

q→ p= ±iω1 , ±iω2 ( B 2 − 4 AC = 0) → stability boundary −σ 1 ± iω1 , σ 2 ± iω2 ( B 2 − 4 AC < 0) → More higher q → p = Higher

Even more higher



q → p = 0, ±iω1 (C = 0) →

dynamic instability stability boundary

Flutter condition: B 2 − 4 AC = 0

Torsional divergence: C = 0

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity Graphically,

= ωα2

Kα 2 K h = , ωh Iα m *Undamped system

U1 = ωα

U 2 = ωh

σ

Im(p)

q qF

0

Re(p)

σ > 0(unstable) q σ > 0( stable)

- Effect of static unbalance In Dowell’s book, after Pines[1958]

2 ω q Sα ≤ 0 → avoid flutter, if Sα = 0, F = 1 − h2 qD ωα Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity If If

qD < 0(e < 0) ωh < 1.0 ⇒ qF < 0 no flutter qD > 0 and

ωα ωh > 1.0 ⇒ ωα

no flutter

- Inclusion of damping→ “can be a negative damping” for better accuracy,

0, where mq + cq + Kq =

 qSCLα  U∞  c=  qSCL α −  U ∞

 0    −qScCmα  

The characteristic equation is now in the form of

A4 p 4 + A3 p 3 + A2 p 2 + A1 p + A0 = 0 Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity A4 p 4 + A3 p 3 + A2 p 2 + A1 p + A0 = 0 ⋅⋅⋅ * • Routh criteria for stability ; At critical position, the system real part becomes zero, damping becomes zero. Substitute

p = iω into (*), we get,

 A4ω 4 − A2ω 2 + A0 = 0  3 i − A ω + A1ω ) = ( 0 3 

A

2 1 From the second eqn, ω = A , substitute into first equation, then,

2

A  A  A4  1  − A2  1  + A0 = 0 or A A  3  3

3

A4 A12 − A1 A2 A3 + A0 A32 = 0

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Dynamic aeroelasticity p as q increases,

And, we can examine Low

q→ p= −σ 1 ± iω1 , −σ 2 ± iω2 → damped natural freq.

Higher

q→ p= −σ 1 ± iω1 , ±iω2

More higher q → p = −σ 1 ± iω1 , ±σ 2 ± iω2 → dynamic instability. ω ω2

ω1

q

qF q

σ q

- Static instability ∙∙∙ | K |= 0 - Dynamic instability a) frequency coalescence (unsymmetric K ) b) Negative damping (Cij < 0) c) Unsymmetric damping (gyroscopic)

σ Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing Consider disturbance from equilibrium

Using modal method, the displacement

( wea )

and rotation (θ ea ) at

elastic axis can be expressed as N   wea = ∑ hr ( y )qr (t )  r =1 where  N θ = α ( y )q (t ) r r  ea ∑ r =1

qr (t ) : generalized (modal) coordinate hr ( y ), α r ( y ) : mode shape N : total number of modes Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing For N = 4,

a )= h1 1,= α1 0 : rigid translation mode (= ω1 0 ) b= ) h2 x= 0 : rigid pitch mode about c.g.= ( ω2 0 ) 0 ,α 2 c) h3 ( y ), α 3 ( y ) : 1st bending of wing (ω3 ≠ 0 ) d ) h4 ( y ), α 4 ( y ) : 1st torsion of wing (ω4 ≠ 0 )

h3

1st bending mode

α3

(bending dominated) y

y

Modes can be assumed, or calculated from mass-spring representation. The displacements and rotations at any point

w ( x, y, t ) = wea + ( x − x0 ) θ ea =

θ ( x, y,= t ) θ= ea

N

∑ h + ( x − x ) α r =1

r

0

r

 qr (t )

N

∑ α q (t ) r =1

r

r

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing The kinetic energy (T ) is

T=

1 2  m w ( ) dxdy 1 ∫∫ aircraft 2 2

N N 1 = m∑  hr + ( x − x0 ) α r  qr ∑  hs + ( x − x0 ) α s  qs dxdy ∫∫ 2 r 1 =s 1 =

1 N N = ∑∑ mrs qr qs 2 =r 1 =s 1 where, mrs=



l

0

 Mhr hs + Iα α rα s + Sα ( hrα s + hsα r )  dy

TE

M = ∫ mdx : mass/unit span LE

= Sα = Iα

∫ ( x − x ) mdx : static unblance/unit span ∫ ( x − x ) mdx : moment of inertia about E.A./unit span TE

LE

0

TE

LE

2

0

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing The potential energy (U ) is 2

U

2

 ∂θ ea  1 l  ∂ 2 wea  1 l EI dy GJ +    dy 2  ∫ ∫ 0 0 2 2  ∂y   ∂y  N N N N 1 l 1 l EI ∑ hr′′qr ∑ hs′′qs dy + ∫ GJ ∑ α r′qr ∑ α s′qs dy ∫ 0 0 2 =r 1 =s 1 2 = r 1 =s 1 1 N N = ∑∑ K rs qr qs 2 =r 1 =s 1

= K rs where,



l

0

l

EIhr′′hs′′dy + ∫ GJ α r′α ′dy 0

′′ h= ′′ 0 and α= ′ α= ′ 0 [Note] K rs = 0 for rigid modes 1,2, since h= 1 2 1 2

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing Finally, the work done by airloads, l

l

N

0

0

r =1

δW = − ∫ Leaδ wea dy + ∫ M eaδθ ea dy − LHT δ wHT + M HT δθ HT = ∑ Qrδ qr subscript where, [Note]

HT : horizontal tail contribution (rigid fuselage assumption)

Qr =∫ ( −hr Lea + α r M ea ) dy − hr ( HT ) LHT + α r ( HT ) M HT l

0

1 1 → Q1 = r= − ∫ Lea dy − LHT = − LTotal 0 2 1 r =2 → Q2 = M Total (C.G ) 2 l

place T , U , and Qr into the Lagrange’s equation

d  ∂T  dt  ∂qr

 ∂T ∂U + = Qr −  ∂qr ∂qr

yield the equation of motion Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing Equation of motion in matrix form

0 0 [ mrs ]{qr } + 0  0

0 0 0 0 0 K 33 0 K 43

0   q1  0  q2  {Qr }  = K 34   q3   K 44  q4 

zeros are associated with rigid body modes N

[Note] If we used normal modes, w ( x, y, t ) = ∑ φr ( x, y ) qr (t ) r =1

free-free normal mode The equation of motion would be uncoupled

  [ mrs ] →  mrr  ,  

  [ K rs ] →  mrrωr2    Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing

[Note] Free-free normal mode

vs

from entire structures

for individual components

M r qr + M rω q = Qr 2 r r

Uncoupled modes

then, combine together by Rayleigh-Ritz method,

∑m

q + ∑ krs qs = 0

rs s

(more accurate)

(more versatile)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Straight Aircraft Wing Now, let’s introduce the aerodynamic load by considering 2-D, incompressible, strip theory

ea + Uθea − baθea  + 2πρUbC (k )  w ea + Uθ ea − b Lea πρ b 2  w =  ea + U 12 − a θea − b 18 + a 2 θea  = M ea πρ b3  aw   +2πρUb 2 12 + a C (k )  w ea + Uθ ea − b 12 − a θ ea   

(

(

ωb U

=

(

)

lift deficiency fn.

∗ k=

)

3 4

) (

(

1 2

)

− a θ ea  

)

c airspeed (downwash)

ωc 2U

x ∗ a =ea b Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Unsteady Aeroelasticity - Unsteady Aeroelasticity in Incompressible Flow (B.A.H p.272 and B.A. p.119)

( M << 1)

- For incompressible flow

a separation can be made between circulatory and non-circulatory airloads - When the airfoil performs chordwise rigid motion. the circulatory lift depends only on the downwash at the

w3 c =  w ea + Uθ ea − b  4

(

1 2

)

− a θ ea  : downwash at 

3 4

1 4

c ”

c station

c

ea + Uθea − baθea  + 2πρUbC (k )  w ea + Uθ ea − b = Lea πρ b 2  w  “always acts at

3 4

(

1 2

)

− a θ ea  

lift deficiency fn. Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Unsteady Aeroelasticity  wea = ∑ hs qs  s However,  θ ea = ∑ α s qs s  and placing these into Lea , M ea yields

Qr =∫ ( −hr Lea + α r M ea ) dy + H .O.T =Qr ( qs , qs , qs ) l

0

coupled set of homogeneous differential equations. pt For stability analysis, assume qr (t ) = qr e

where

p= σ + iω , and for a) + σ , ω ≠ 0

"flutter" t

b) + σ , ω = 0

"divergence" t Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M - Solutions of the Aeroelastic Equations of Motion (Dowell pp.100~106) - Two Groups: a) Time domain and b) Frequency domain a) Time domain: fundamentally, a step by step solution for the time history ∙ Direct integration method ① equilibrium satisfied at discrete time t ② assumed variation of variables

( q, q , q)

within the time interval ∆t ∙ Examples of methods ① central difference ② Newmark ③ Houbolt Ref. Bathe, “Finite Element Procedures”, Chap. 9 Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M ∙ When selecting a method, three main issues to be aware ① efficient scheme ② numerical stability conditionally stable – dependent on ∆t unconditionally stable ③ numerical accuracy amplitude decay period elongation ∙ Advantage and disadvantage of time domain analysis ① Advantage: straight forward method ② Disadvantage: aerodynamic loads may be a problem → theories are not well-developed → intensive numerical calculation for small number of frequency ( k ) Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M b) Frequency domain: most popular approach ∙ Main issue: aerodynamic loads are well developed for simple harmonic motion iωt ∙ consider simple harmonic motion qr (t ) = qr e

and corresponding lift and moment, (Ref. Drela, last page)

Lea = Lea eiωt M ea = M ea eiωt w   = Lea πρ b3ω 2 lh ( k , M ∞ ) ea + lα ( k , M ∞ ) θ  b   where, w   M ea πρ b 4ω 2  mh ( k , M ∞ ) ea + mα ( k , M ∞ ) θ  = b   lh , lα , mh , mα are dimensionless complex fn. of ( k , M ∞ ) (Refs. Dowell, p.116 and B.A. pp. 103~114) Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M ∙ Then, the governing equation becomes

0 −ω 2 [ M ]{q } + [ K ]{q } + ω 2  A ( k , M ∞ )  {q } = aerodynamic operator (aero. mass matrix) It is presumed that the following parameters are known.

(

M , Sα , Iα , ωh , ωα , b = 12 c inertia

)

stiffness

The unknown quantities are

 ωb  q ,ω, ρ , M ∞ , k  =   U  determined by p

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M I) k-method (V-g method) ∙ consider a system with just the right amount of structural damping, so the motion is simple harmonic

−ω 2 [ M ]{q } + (1 + ig ) [ K ]{q } + ω 2 [ A]{q } = 0 structural damping coefficient [Note] structural damping – restoring force in phase with velocity,

F0 = − g ( q q ) q phase

but proportional to displacement

g required > g available : unstable

displacement

* viscous damping - Fc = −cq

g required = g available : neutral g required < g available : stable

∙ Rewrite equation

(1 + ig ) [ M − A]{q } =

ω 2 [ K ]{q }

Re[Λ ] 1 ω 2 ,= Im[Λ ] g ω 2 Λ,= Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M ∙ Solution process ① Given M , Sα , Iα , ωh ωα , b ② Assume ρ (fix altitude), M ∞ = U a∞ ③ For a set of k values, solve for eigenvalues for Λ

ω

ωα

ωh g

g1

0.03

1 kF

g2

U 1 = ωb k g available ≈ 0.03 > 0 U 1 = ωb k

“the only valid point”

④ for g1 =0 → ω =ωF ( k F =bωF U F ) ⑤ matching problem

UF → M F = M∞ Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M I) k-method (V-g method) (Dowell, p.106) 2 2 ∙ Structural damping is introduced by multiplying at ωh , ωα

× (1 + ig ) , g : structural damping coefficient pure sinusoidal motion is assumed → ω ≡ ωR , ωI ≡ 0 for a given U , the g required to sustain pure sinusoidal motion is determined ∙ Advantage – the aero. force need to be determined for real frequencies

= U U= ω= 0) ∙ Disadvantage – loss of physical sight, only at F (ω R , ωI the mathematical solution will be meaningful ∙ Following parameters are prescribed

M , Sα , Iα , ωh ωα , k , m 2 ρ∞ bS then, the characteristic equation becomes a complex polynomial in unknowns

(ωα

ω )(1 + ig )

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M I) k-method (V-g method) (Dowell, p.106) ∙ A complex roots are obtained for ωα ω and g From ωα ω and the previously selected k ≡ ωb U ,

ωα b U∞

=

ωα k ω

Then, plot g vs U ∞ bωα (typical plot for two d.o.f system below)

g : value of structural damping required to sustain neutral stability → If the actual damping is g available, then flutter occurs when g = g available

g

I

g available

ω ωα

U bωα

UF bωα

II If g < g available , U < U F → no flutter will occur

II

I U bωα

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M I) k-method (V-g method) (Dowell, p.106) ∙ Uncertainty about g available in a real physical system, flutter speed is defined as minimum value of U F bωα for any g > 0

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M pt II) p-method – time dependent solution q= qe , p= σ + iω

∙ The equation,

 A ( p, M )  {q } p 2 [ M ]{q } + [ K ]{q } =

Now the aero becomes more approximate i)

quasi-steady aero

−σ

ii) induced lift function

“all points are

iii) flow Eigen solution

valid ones!!”

1 k

iωt [Note] I) k-method (V-g method), qr = qr e

only valid for single harmonic motion – k ~ ω pt II) p-method, q= qe , p= σ + iω

 A ( p, M )  [ M ]{q} + [ K ]{q } =

- “true damping” (H. Hassig)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M III) p-k method ∙ The solution is assumed arbitrary (as in p-method) However, the aero. is assumed to be

A ( p, M ) ≅ A ( k , M ) Then, the eqn. becomes:

0 { p [ M ] + [ K ] −  A ( k , M )}{q } = 2

∙ Solution process ① specify ki , M i

σ 0 + iω0 ② solve for p= 0 k0

③ check for double matching

k 0 = ki M F = Mi Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Solutions of the Aeroelastic E.O.M [Note] p-k method usually requires handful of iteration to converge. It is more expensive than k-method. ∙ Alternative: p-k method (Dowell)

h, α ~ e pt is assumed, p= σ + iω in aero. terms, only a k ≡ ωb U is assumed The eigenvalues p are computed → new ω → new k → new aero. terms – iteration continues until the process converges For small σ , i.e., σ << ω , σ ~ true damping solution

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Padé Approximation Method The generalized forces Qr are computed for harmonic motion

Qr = 12 ρU 2Qrs qs eiωt

(πρω

= Qrs where

( Qrs )real + i ( Qrs )img :

2

Ars qs eiωt

)

complex function of reduced frequency

( Qrs )img

( Qrs )real

k

k one can fit above by Padé Approximation in Laplace transform domain p of from

 b U) p  ( 1 2 2 2 Qr ρU  A2 ( b U ) p + A1 ( b U ) p + A0 + A3 =  qs 2 ( b U ) p + β1   mass

damping

stiffness

lag

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Padé Approximation Method For harmonic motion p = iω

 β1k   1 k2   2  = Qr ρU  − A2 + A0 + A3 2 q  + i  A1k − A3 2 3  s 2 + + k β k β 1  1   

( Qrs )real

( Qrs )img

and then evaluate coefficients A2 , A1 , A0 , A3 , β1 to fit Qrs over certain range of k , 0 ≤ k ≤ 2 ( k ≡ ωb U ) [Note] for better fit, use more lag terms,

Qr

N   b U) p ( 1 2 2 2 ρU  A2 ( b U ) p + A1 ( b U ) p + A0 + ∑ Am  qs 2 b U p β + ( ) m =3  m−2  

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Padé Approximation Method Next, introduce new augmented state variables ys , defined as

(b U ) p q p qs = s p + (U b ) β s (b U ) p + βs

ys =

pys + (U b ) β s ys = pqs Return to time domain,

= Qr

1 2

ρU 2  A2 ( b U ) qs + A1 ( b U ) qs + A0 qs + A3 ys  2



y s + (U b ) β s ys = qs



and governing equation,

= Mq + Cq + Kq

1 2

ρU 2  A2 ( b U ) qs + A1 ( b U ) qs + A0 qs + A3 ys  2

    y =  y s +  Uβ b  s qs   Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Padé Approximation Method or

M *   0   0

0 M* 0

0   q   0    0   q +  K * I   y   0

−M * C* −I

0  q    G   q  = 0 H   y 

*  M= M − 12 ρ b 2 A2  1 * = − C C ρ bA1  2 *  K= K − 12 ρU 2 A0  G = 12 ρU 2 A3 where       H =  U β b      

 q    = and then,  q  y   

q  [ A]  q   y  



 = AX X

Ref.: Karpel, minimum-state (1991) Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Types of Flutter I) “Coalescence” or “Merging frequency” flutter ∙ coupled-mode, bending-torsion flutter (2 d.o.f flutter) ∙ for U > U F , one of ωI → (+ ) and large (stable pole) the other ωI → (−) and large (unstable pole)

ωR remain nearly the same torsion mode being unstable  ∙ although   the airfoil is  bending mode being stable  undergoing on oscillation composed of both

ωα ω σ ωα

→ torsional mode usually

U bωα U bωα

goes unstable → flutter mode contains significant contributions of both bending and torsion

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Types of Flutter (Dowell. P.103) I) “Coalescence” or “Merging frequency” flutter ∙ the “out-of-phase” (damping) force are not qualitatively important → may neglect structural damping entirely and use a quasi-steady or even a quasi-static aerodynamic assumption → simplified analysis

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Out-of-Phase Force (BAH p.528) - 2-D rigid airfoil with a torsional spring (1 d.o.f system)

Iα α + Kα α = My by assuming

α = α o eiωt Iα πρ b 4

where

my =

  ωα  2  0 1 −   + my =    ω  

My

πρ b ω α o e 4

2

iωt

, function of only k = ωb

U

Substituting into (1), flutter occurs when the out-of-phase aerodynamic damping component vanish. - Rotating complex vector diagram

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Out-of-Phase Force (BAH p.528) - Rotating complex vector diagram

α

− Kα α My

− Iα α

U < UF

Im {M y } < 0

U > UF

Im {M y } > 0

Im {M y } = 0

Im {M y }

M y which lags that motion ( Im {M } < 0 ), removes energy from the y oscillation, providing damping.

{ }

This out-of-phase component, Im M y

, is the only source of

damping or instability from the system.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Types of Flutter II) Single d.o.f. flutter ∙ frequency of mode almost independent of reduced velocity ∙ results from negative damping ∙ out-of-phase part of aerodynamic operator is very important ∙ typical of systems with large mass ratio at large reduced velocity (e.g. turbo machinery, bridge, …)

ω ωα

ωh ωα

σ ωα

U bωα U bωα

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Types of Flutter (Dowell. P.103) II) Single d.o.f. flutter ∙ frequencies, ωR , independent of the airspeed

(U

bωα ) variation

∙ true damping, ωI , also moderate change with airspeed ∙ one of the mode (usually torsion) becomes slightly (−) at U F → very sensitive to structural and aerodynamic damping forces → since those forces are less precisely described, analysis gives less reliable results ∙ Since the flutter mode is virtually the same as that of the system at zero airspeed, the flutter mode and frequency are predicted rather accurately (mass ratio < 10) ∙ Airfoil blades in turbo machinery and bridges in a wind.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Types of Flutter III) Divergence ∙ flutter at zero frequency ∙ very special type of single d.o.f. flutter ∙ out-of-phase forces unimportant ∙ analysis reliable

ωα ω

σ ωα

U bωα

U bωα

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Parameter Effects on Wing Flutter ∙ When one non-dimensionalizes the flutter determinant (2D), 5 parameters will appear:

m mass ratio = 2 πρ b Sα distance CG aft of EA xα = = mb b Iα radius of gyration about EA γα = = mb 2 b e distance EA aft of midchord a= = b b = µ

ωh = uncoupled bending to torsion frequency ratio ωα

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Parameter Effects on Wing Flutter [additional]

ωα t = nondimensional time M = Mach. No. (compressibility effect) ωα b = Kα = reduced frequency U 1 = reduced velocity   ω UF = f  µ , xα , γ α , a, h , M  bωα ωα  

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends The trends are: a)

, (CG. Ahead of EA) - frequently no flutter

b)

1

0.3

M

dip can be quite severe and approach to zero

- Structural damping can remove dip completely

c) 2-D airfoil theory

5 - For large

,

10

constant, for small ,

helicopter rotor blade

constant (dashed line) Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends

d)

Various 1

altitude flutter

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Approximate Formula An approximate formula was obtained by Theodorsen and Garrick for small large .

Distance (non-dimensional) between AC and CG (B.A.H. 9-22) Recall divergence:

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends . non dimensionalize the typical section equation of motion

S I h m e ωh U M = F1 (ωα t : α , α 2 , , , , , ) 2 b mb mb ρ (2b) b ωα bωa

α = F2 (ωα t :

Sα Iα m e ωh U M , 2, , , , , ) 2 mb mb ρ (2b) b ωα bωa

- Choice of non-dimensional parameters: . not unique, but a matter of convenience i) non dimensional dynamic pressure, or ‘aeroelastic stiffness No.’

U 1 4 ρU 2 instead of a non dimensional velocity, λ≡ = µ Kα 2 mωα 2 bωα 2

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends ii) ωα t

nondimensional time

Sα mb I γα 2 ≡ α 2 mb m µ≡ ρ (2b) 2 e a≡ b

γα ≡

ωh ωα

radius of gyration (squared) mass ratio location of e.a measured from a.c or mid-chord frequency ratio Mach number

M ka =

static unbalance

ωα b U

Reduced frequency

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends -

For some combinations of parameters, the airfoil will be dynamically unstable (‘flutter’)

-

Alternative parametric representation Assume harmonic motion = h Eigenvalues

= ω ωR + iωI

ω ωR U ) = GR ( xα , rα , µ , a, h , M , bωα ωα ωa -

iωt he = , α α eiωt

,

ω ωI U = GI ( xα , rα , µ , a, h , M , ) bωα ωα ωa

For some combinations, ωI < 0 , the system flutters.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends ― the coalescence flutter , conventional flow condition (no shock oscillation and no stall) I) Static unbalance, xα … if xα < 0 , frequently no flutter occurs

ωh ω II) Frequency ration … U F bωα is a minimum when h ≈ 1 ωα ωα III) Mach No. M … aero pressure on an airfoil is normally greatest near M = 1 → flutter speed tends to be a minimum U2 For M >> 1, from aero piston theory, p ≈ ρ M 12 For M ≥ 1 and constant µ , U F ≈ M

― for flight at constant altitude,  U 

α



ρ (hence µ ) and α



(speed of sound)

 = M  ∞  → compatibility relation fixed. U = Mα ∞ →   bωα   bωα  U F bωα compatibility relation

≅ M 1 2 for large M 1.0

and for a flat plate

M Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Parameters Trends ― by repeating flutter calculation for various altitudes (various

ρ ,α ∞ ,

various µ and α ∞ bωα ) Altitude flutter

No flutter

IV) Mass ratio µ … For large For small

M

µ → constant flutter dynamic pressure µ → constant flutter velocity (dashed line)

for M ≡ 0 and 2-D airfoil theory → U F → ∞ for some small but finite µ (solid line)

UF bωα

2

1U  λF =  F  ≅ µ  bωα 

constant for large

µ

µ Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter Prevention - Flutter Prevention - add mass or redistribute the mass (“mass balance”)

xα < 0

- increase - move

away from 1

- add damping, mainly for single D.O.F flutter - use composite materials • couple bending and torsion • shift ωα away from ωh - limit flight envelope by “fly slower”

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Physical Explanation of Flutter (BA p. 258) • Purely rotational motion of the typical section

Iα α + Kα α = My

∂M y ∂M y    π 3 1 2  α +  Kα − α= 0 ― Approximate form:  Iα + ρ 0b S  + a   α −  ∂α ∂α  2 8    ∂M y ∂M y if ∂α , ∂α are known,

→ second-order, damped-parameter

system with 1DOF 2 ― Laplace transform variable p, characteristic polynomial a0 p + a1 p + a2 two possible ways of instability I)

α coeff. (+) → (−), a2 ≤ 0 in Routh’s criterion → “torsional

divergence” …negative “aerodynamic spring” about E.A. overpowers Kα

II)

∂M y ∂α

(−) → (+), a1 ≤ 0 in Routh’s criterion → dynamic instability

entirely aerodynamic “negative” damping I m {M y } = 0

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Physical Explanation of Flutter • Qualitative explanation of negative damping ― principal part Lo … due to the incremental a.o.a, in phase with α, e.a. at ¼ chord → adding to the torsional spring Kα when a < − 12 ― bound circulation Γo … changing with time. Since the total circulation is const., countervortices strength are induced shed from the trailing edge → wake vortex sheet → out-of-phase loading is induced (upwash) at low k ― upwash… produces additional lift L2 ⇒ when e.a. lies ahead of ¼ chord, the moment due to L2 is in the same sense of α → net positive work per cycle of the wing “negative damping” ― at higher k , damping becomes (+) more cycles of wake effects upwash, bound circulation lags behind α , center of pressure of lift oscillates Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Physical Explanation of Flutter (BA p.258) i) Pure rotational system (1 D.O.F) Iα α + Kα α = My

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

2 D. O. F. system ∂G  h  2  m h + ωh h + Sα α = − qS α +  ∂α  U ∂G  h  2   Sα h + Iα α + ωα α = qS e α +  ∂α  U

[

]

[

]

 1  + b a   2   e= b  a + ( γ + 1) M Aw  Piston theory 2    4 2 b  

―Dimensionless frequency and damping I) U = 0 ≈ 1 2 critical U bωα … mode shape remains the same as for free vibration, involving pure rotation about an axis II) rotation axis moves forward, as indicated by falling amplitude of bending III) gradual suppression of h… caused by lift variation due to torsion, lift, in phase with α, drives the bending freedom at ω greater ωh → response to it has a maximum downward amplitude at the instant of maximum upward force

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

2 D. O. F. system ―Dimensionless frequency and damping IV) simultaneously ω drops… lift constitutes a negative “aerodynamic spring” on the torsional freedom with “spring constant” ~ dynamic pressure V) small advances in

ϕ h… due to lift,

due to

h

VI) flutter occurrence … bending amplitude =0, only pure rotational oscillation about E.A., no damping acts

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F (BAH p. 532)

― flutter from coupling between the bending and torsional motions the most dangerous but not the most frequently encountered ― Equations of motions

2  + S α   mh α  + mω h h = Qh = − L   2   α ω S h I I + +  α α α α = Qα = M y  α

― Simple harmonic motion

b

h = h0eiωt , α = α 0ei (ωt +ϕ ) = α 0eiωt 2 2 2   − ω mh − ω Sα α + ωh mh = − L ⇒ 2 2 2 S h I Iα α = M y ω ω α ω − − +  α α 

b

b

ba bxα Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F ― Aerodynamic operator

 h  1   L = −πρb 2ω 2  Lh +  Lα − Lh ( + a )α  2    b   1 1 1 h    M y = −πρb 2ω 2  M h − Lh ( + a ) +  M α − ( Lα + M h )( + a ) + Lh ( + a ) 2 α  2 2 2 b    

function of Lh , Lα , M α (incompressible) K , M α 1 2

Plugging the aerodynamic operator, and set the coefficient determinant to zero • characteristic eqn. for ωα ω … implicitly dependent on the 5 dimensionless system parameters

a : axis location

ωh ω α : bending-torsion frequency ratio xα = Sα mb : dimensionless static unbalance

rα = Iα mb 2 : radius of gyration

m πρb 2 : density ratio • parametric trends of U F in terms of 5 parameters Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F ― Divergence speed U p

Up

2 Kα = ρ eccl α

Kα 1  2πρ b 2  + a  2  1 b( + a ) 2

UD bωα

from a.c. to e.a.

Kα b 2 Iα 1 m U= p bωα Iα mb 2 πρ b 2 [1 + 2a ]

rα2 m πρ b 2 [1 + 2a ]

both UD above and the flutter speeds in Fig 9-5 from the 2-D aerodynamic strip theory → the predicted U F will not exceed U D

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F - Fig. 9-5 (A)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F - Fig. 9-5 (B)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F - Fig. 9-5 (C)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F - Fig. 9-5 (D)

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Flutter of a simple system 2 D.O.F Fig. 9-5(A),(C)… dip near ωh ωα ≅ 1 → can bring up with small amounts of structural friction (B)… density ratio increase → raise flutter speed (flutter speed vs. altitude) “mass balancing”… flutter speed is more sensitive to a change of → Not much balancing is needed to assure safety form bending-torsion flutter Fig. 9-5(D)… flutter is governed by (a + xα ) chordwise c.g. Garrick and Theodorsen (1940):

rα2 UF m ≈ bωα πρ b 2 1 + 2(a + xα ) 2   From a.c. to c.g. Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter - Panel Flutter: - Self-excited oscillation of the external skin of a flight vehicle when exposed to airflow on that side (supersonic flow)

simplify • For simplicity, consider a 2-D simply supported panel in supersonic flow; for a linear panel flutter analysis, the equation of motion is: , where

Eh3 D= 12 (1 −ν 2 )

(isotropic, plate stiffness)

mass/unit,

: thickness

= aerodynamic pressure Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter

• Putting all together, the governing equation becomes:

• It is subject to :

• They are the simply supported B.C • Using Galerkin Method • satisfies all the B.C.’s Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter •

By setting:



we get:

Anti-symmetric

• where

: speed of sound,

: critical speed param.

: lowest natural frequency

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter - A typical result :

[Note] If

constant

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter Theoretical considerations of panel flutter at high supersonic mach numbers (AIAA J, 1966) •Basic Panel Flutter Eqn. and its Sol. •A rectangular panel simply supported on all 4 edges and subject to a supersonic , elastic foundation flow over one side, midplane compressive force structural damping

•Aerodynamic pressure for high supersonic Mach No.

(1)+(2) : non-dimensional coordinates introduced

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter • Additional non-dimensional parameters

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter • Simply supported B.C’s • Solution procedure

O.D.E

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter General solution of O.D.E This along with the B.C, the determinant must be: - Equal to zero for nontrivial solutions

- For low values of the determinant, the eigenvalues are real. - Above a certain value of A, they become imaginary.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Panel Flutter Complete panel behavior - Plotting - The Frequency coalescence: Instability occurs at

- Flutter Frequency:

- Deflection shapes •Simple sine shape standing-wave type for A=0 •Standing-wave type at low A •Traveling-wave type at high values of A

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Computational Aeroelasticity

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Computational Aeroelasticity Structure

Aerodynamics

• With the abundance of computational resources and algorithms, there has been a great development in two areas: • CFD: Computational Fluid Dynamics • CSD: Computational Structural Dynamics • CAE: Computational Aeroelasticity

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Computational Aeroelasticity •

Difficulties arise from the nature of the two methods. •

CFD: Finite difference discretization procedure based on Eulerian (spatial) description

• CSD: finite element method based on Lagrangian (material) description. • Define the nature of the coupling when combining the two numerical schemes.

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Computational Aeroelasticity i) Tightly (or closely) coupled analysis: - Most popular - Interaction between CFD and CSD codes occurs at every time step - Guarantee of convergence and stability ii) Loosely coupled analysis: - CFD and CSD are solved alternatively - Occasional interaction only => Difficulties in convergence iii) Intimately coupled (unified) analysis: - The governing equations are re-formulated and solved together Active Aeroelasticity and Rotorcraft Lab., Seoul National University

Computational Aeroelasticity i) – Tightly (or closely) coupled analysis: CFD

Initial

Solve Governin g Equation

t ← t+∆t

Compute Pressure Distribution

NO

CSD A

Integrate Surface Loads

Governing Equation

Update CFD Grids

Surface Deflection

Equilibrium ?

YES A

Active Aeroelasticity and Rotorcraft Lab., Seoul National University

End of Chapter III

Active Aeroelasticity and Rotorcraft Lab.