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 ∙∙∙ selfexcited 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 ∂qi
∂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)
Sα
(static unbalance)
Iα
(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 quasisteady 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)
Quasisteady aero: only circulatory terms due to the bound vorticity. Used for characteristic freq. below 2Hz (e.g., conventional dynamic stability analysis)
ii)
Quasiunsteady 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: “quasiunsteady”+”apparent mass terms” (noncirculatory 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 nontrivial 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 massspring 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 qr ∑ hs + ( x − x0 ) α s qs dxdy ∫∫ 2 r 1 =s 1 =
1 N N = ∑∑ mrs qr qs 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 ∂qr
∂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 ]{qr } + 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
freefree 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] Freefree normal mode
vs
from entire structures
for individual components
M r qr + M rω q = Qr 2 r r
Uncoupled modes
then, combine together by RayleighRitz 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 2D, 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 noncirculatory 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 , qs , qs ) 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 welldeveloped → 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) kmethod (Vg 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) kmethod (Vg 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) kmethod (Vg 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) kmethod (Vg 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) pmethod – 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)
quasisteady aero
−σ
ii) induced lift function
“all points are
iii) flow Eigen solution
valid ones!!”
1 k
iωt [Note] I) kmethod (Vg method), qr = qr e
only valid for single harmonic motion – k ~ ω pt II) pmethod, 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) pk method ∙ The solution is assumed arbitrary (as in pmethod) 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] pk method usually requires handful of iteration to converge. It is more expensive than kmethod. ∙ Alternative: pk 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 ) qs + A1 ( b U ) qs + A0 qs + A3 ys 2
y s + (U b ) β s ys = qs
and governing equation,
= Mq + Cq + Kq
1 2
ρU 2 A2 ( b U ) qs + A1 ( b U ) qs + 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, minimumstate (1991) Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Types of Flutter I) “Coalescence” or “Merging frequency” flutter ∙ coupledmode, bendingtorsion 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 “outofphase” (damping) force are not qualitatively important → may neglect structural damping entirely and use a quasisteady or even a quasistatic aerodynamic assumption → simplified analysis
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
OutofPhase Force (BAH p.528)  2D 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 outofphase aerodynamic damping component vanish.  Rotating complex vector diagram
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
OutofPhase 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 outofphase 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 ∙ outofphase 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 ∙ outofphase forces unimportant ∙ analysis reliable
ωα ω
σ ωα
U bωα
U bωα
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Parameter Effects on Wing Flutter ∙ When one nondimensionalizes 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) 2D 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 (nondimensional) between AC and CG (B.A.H. 922) 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 nondimensional 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 midchord 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 2D airfoil theory → U F → ∞ for some small but finite µ (solid line)
UF bωα
2
1U λ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,
→ secondorder, dampedparameter
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 → outofphase 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 ω α : bendingtorsion 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 95 from the 2D 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. 95 (A)
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Flutter of a simple system 2 D.O.F  Fig. 95 (B)
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Flutter of a simple system 2 D.O.F  Fig. 95 (C)
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Flutter of a simple system 2 D.O.F  Fig. 95 (D)
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Flutter of a simple system 2 D.O.F Fig. 95(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 bendingtorsion flutter Fig. 95(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:  Selfexcited oscillation of the external skin of a flight vehicle when exposed to airflow on that side (supersonic flow)
simplify • For simplicity, consider a 2D 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:
Antisymmetric
• 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) : nondimensional coordinates introduced
Active Aeroelasticity and Rotorcraft Lab., Seoul National University
Panel Flutter • Additional nondimensional 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 standingwave type for A=0 •Standingwave type at low A •Travelingwave 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 reformulated 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.