METODA KOMPUTASI DALAM ANALISA STABILITAS STATIK DAN DINAMIK PESAWAT TERBANG NON-CONVENSIONAL Erwin Sulaeman*
Acknowledgement
Metoda Komputasi Dalam Analisa Stabilitas Statik dan Dinamik Pesawat Terbang NonConvensional
Rakesh K. Kapania Liviu Librescu Virginia Polytechnic Institute and State University
Erwin Sulaeman Indonesian Aerospace, IAe
Raphael T. Haftka University of Florida 2
1
*
Metoda Komputasi
Overview
Stabilitas Statik
• Definisi
K
– Stabilitas Statik dan Dinamik – Pesawat Terbang Non-Convensional
• Structural Finite Element Method
Buckling
P Kg
u
0
FEM
– Analytical Formulation for MDO
Stabilitas Dinamik
• Aerodynamic Boundary Element Method – Analytical Integration for Singular Kernel Function
M
u!!
C
u!
K
P Kg
u
Flutter Am
u!!
Ac
u!
Ak
u
• Stabilitas Dinamik Sayap : Flutter • Kesimpulan
FEM
3
Non-Conventional Aircrafts
BEM
4
SBW as an MDO Problem
• 325 Passenger Transports 1995 Conventional
2010 Strut-Braced Wing
2010 Box Wing The strut reduces wing bending moment
• Larger spans • Reduced t/c • Reduced weight Ramp
714,499 lb.
Ramp
504,835 lb.
Ramp
535,963 lb.
Fuel
295,415 lb.
Fuel
169,235 lb.
Fuel
186,821 lb.
Zero Fuel 419,084 lb.
Zero Fuel 335,600 lb.
Zero Fuel 349,142 lb.
Pax/Cargo 68,250 lb.
Pax/Cargo 68,250 lb.
Pax/Cargo 68,250 lb.
Op Equip
Op Equip
Op Equip
Empty
*
17,779 lb. 333,055 lb.
Empty
17,680 lb. 249,670 lb.
Empty
Indonesian Aerospace, IAe PT Dirgantara
17,603 lb. 263,289 lb.
• Less induced drag • Less sweep angle, more laminar flow regions, less friction drag • Reduce engine weight and fuel consumption 5
6
Description of the MDO Process
Optimum Configurations
Baseline Design
• Mission profile follows a typical Boeing 777 200IGW mission
Induced
Initial Design Variables
Drag
Updated Design Variables Geometry Definition
Friction and Form Drag
Structural
Propulsion
Aerodynamics
– 7500 n m i range + 500 n m i reserve range – Mach 0.85 cruise – 325 passengers – 2 GE- 90 class engines
Wave Drag
Optimization Wing bending material weight
SFC
Interference L/D
547,000 lbs.
Drag
Weights Offline CFD
521,000 lbs.
Analysis Range/ Field Performance
Performance
Stability and
Objective Function/
Control
Constraints
* Structural Optimization includes static aeroelasticity
608,000 lbs. 524,000 lbs. 1
Optimizer
8
Proposed Scheme for Computational Aeroelasticity M K M
u
P Kg
u
K
P Kg
0
K
P Kg
P Kg
qc A
u
F
K
P Kg
qd
A
u
0
K
P Kg
u
C
u
q A
u
F
M K
Buckling
K
M
u
C u
Proposed Scheme for Computational Aeroelasticity
u
Modal Analysis
0
q A
u
P Kg
C u
u
K
P Kg
0
q A
u
F
Buckling 0
Modal Analysis
Passive Load Alleviation
K
P Kg
qc A
u
F
Passive Load Alleviation
Divergence
K
P Kg
qd
A
u
0
Divergence
K
P Kg
u
0
M
M
Flutter
u
u
K
C
P Kg
u
u
q A
u
0
Flutter
9
Computational Aeroelastic Research Topics
10
PART 1: Structural Analysis
• Structural Optimization
– Fully-stressed design concept to predict the wing bending weight. Double Plate Model
• Passive Load Alleviation with Geometric Stiffness Effect
– Considering structural flexibility, buckling, and aerodynamic load redistribution (follower force). – B -i linear strut stiffness-displacement diagram
• Flutter Analysis using Detailed FEM Model
– Wing Section Model : Hexagonal Box
• Structural FEM Formulation
– Exact Formulation for non-prismatic beam element – 2 nd Order Accurate Geometric Stiffness Formulation
• Unsteady Aerodynamic Formulation
– Hybrid doublet lattice/doublet point method for non-planar interfering compressible lifting surfaces
• Modified K / PK Flutter Solutions
– Include additional buckling eigenvalue problem
– Computational off-line using Nastran
• The stiffness matrix [ K ] – Formulated analytically for a non-uniform beam element – The beam stiffness varies as an arbitrary polynomial function
• The geometric stiffness [ Kg ] – The initial buckling load is included in the formulation to increase the accuracy
– Fast computational aeroelasticity as part of MDO code 11
12
Hexagonal Box Model
Wing Box Structure Model 0.04
t1
0.03
t
• 1 4 S u b -component Variables – 4 skins : t 1, t 2 , t 4 and t5 – 2 shear webs : t3 and t6 – 4 spar caps: c1 , c 2, c 3 and c4 – 4 stringers: s1 , s 2, s 3 and s4
t2
0.02 0.01 z/c
t3
t6
0 -0.01 -0.02 -0.03
0
0.5
0
1
Hexagonal wing box model with optimized thickness and area ratios
Detailed FE Model which can be Analyzed by NASTRAN
t4
t5
-0.04
E q u i v a l e n t n o n -u n i f o r m beam model
Double Plate Model
0.2
0.4
0.6
0.8
1
•
Hexagonal Box Model Formulation
– Better simulate the airfoil geometry – Give bending and torsional stiffnesses
•
Lockheed Martin Hexagonal Box Model
– Assign only one independent variable – The other 13 variables are set as dependent variables – The aspect ratio between 14 variables are optimized offline
•
Present Equivalent Hexagonal Box Model
– Follows the Lockheed Martin model – The independent variable is obtained from the double plate model result – Three equivalent hexagonal box models are proposed
13
14
Wing Stiffness Distribution
Non-uniform Beam Finite Element Formulation
10
Hexagonal Model 1
EI(x) == EI
Karabalis and Beskos (1983)
1 EI (x 1e9 lb sq. ft)
Beam stiffness distribution function
Method
Hexagonal Model 2 0.1
EI double plate model
EI(x) == EI
EI(x) == EI
EI(x) ==
+
EI hexagonal model 3
Present
∑ i∑ ==1
pi
Mk 0.4 0.6 y/halfspan
0.8
EI
1
N/A
{{1 ++r (( x/L )) }}
0
m p mq
Nk
0.001 0.2
{{ 1 ++ r ((x/L )})} m
Np
p(x) == si x i 1 0 < x < Li==(full span)
∑ ∑
E Iy z = 0
Hexagonal Model 3
EI hexagonal model 2
0
0
E Iy z = 0
Baker (1996)
EI hexagonal model 1
0.01
N/A
E Iy z = 0
Banerjee and Williams (1986)
+
Load distribution
{{ 1 ++ r ((x/L )})} m
0
15
Present Formulation of the Stiffness Matrix
yz (x)
==
∑ ∑ i= =1
xi
mj
N
== EI c
∏ ∏ j ==1 ((
x −− c
j
))
M
i
q i x == EI
yz0
∏ ∏ j ==1 ((
x −− d
j
Np
p(x) == nj
s x
i ∑ i∑ ==1
i
0 < x1 < x < x 2 < L
))
16
Stiffness Matrix Formulation P
s1 [K s]
1 Da
s2
s1
s4
s3
s
s5
s1 Symm .
a o a 1 x a 2 x 2 ... a n x n Beam element with n o n-u n i f o r m E I d i s t r i b u t i o n
2
[ Ks ]
s2 s3
s1 s2 s1
Symm .
s4 s5 s4 s6
Construct the stiffness matrix u s i n g t h e f1 , f2 , a n d f3
L
s4 s6
E I (x )
s1
M
EI ( x ) S o l v e t h e E u l e r-B e r n o u l l i e q u a t i o n for the cantilever beam problem
s1 = f 3 s2 = f 2 s3 = f 1
; ; ;
s4 = s1 L − s2 s5 = s2 L − s3 s6 = s4 L − s5
f1 f2
r1 x
N mj 1
N
G x, P , M
17
v x,P,M v x, P, M
G x , P, M
d
jk
A j,k
jk
B j, k
j 1 k 1 N mj 1
a j B j , mj j 1
f1 = v (0 , 1, 0 ) f 2 = v ( 0, 0 , 1 ) f3 = v′ ( 0 ,1 , 0 )
... rn x n
a j A j,m j j 1
f2 f3
Construct the (2 x 2) flexibility matrix for the cantilever beam
r2 x 2
N
g x , P, M
D a = f 1 f 3 − f 22 [F ]
ro
g x , P, M G L , P, M
d j 1 k 1
g L , P, M x g x,P,M
L g L , P, M 18
Present Formulation of the Geometric Stiffness Matrix
Validation : Tapered Beam 4.0 Distributed
Nastran, tip deflection
moment
load
3.5
Taylor series expansion
Exact Solution
n2 n3
[ K g]
n1 n2 n1
Symm.
n1
t 2 L D1
n2
1 4 D1
t Cos t
3 Sin t
t Sin t 4 Cos t
2t t
2
4
;
n3
L D3 4 D2
;
n4
n2 L
n2 n4 n2 n3
Nastran, tip rotation
n3
(y / L )m }
EI ( y ) = EI 0 1 + r
Proposed FEM 2.5
P0 L2 EI
t =
{
3.0
Error (%)
n1
r = 8, m = 1
2.0 1.5
Method
δδ
Elements
1.0 Exact 0.5
26.19722
Present FEM
1
26.19722
1
26.15411
2
26.19047
4
26.19665
8
26.19720
0.0
•
Nastran
[ K g ] is as also a function of P 0 which can be updated to
-0.5
increase the accuracy •
0
For P 0 = 0 , [ K g ] is the same as that of Gallagher and Lee
2
4
6
8
Number of elements 19
20
Validation : Tapered Beam (2)
Validation : Frame Structure
Distributed load
Tip deflection of the beam
Present (for all n)
60 Error (%)
L
NASTRAN
80
n=1
EI ( y ) == EI 0
n=2
r = - 1/2, L = 8
z
• All prismatic beam elements, except the beam 26
100
{{ 1 ++ r (( y / L )) }}n
8
7 6
5 8
y
3
4
n=3
Tip deflection forn = 5
n=4
40
9
n=5
element #
0
Exact Present FEM
1
8192.00000 8192.00201
1
3052.33000
4
7567.30300
8
8023.43500
16
8148.61800
Nastran
-20 0
2
4
6 8 10 12 Number of Elements
14
16
NASTRAN
x
2 5
Deformations at Point 1
Method
20
1
Beam 26 1 2 4
Tx -0.3083 0.7394 0.9161
Ty -15.2706 -16.3194 -16.4962
Tz 5.8374 6.1429 6.1943
Rx 0.1786 0.1911 0.1932
Ry 7.0130 7.3576 7.4157
Rz 4.4850 4.8062 4.8605
8
0.9318
-16.5120
6.1989
0.1934
7.4208
4.8653
1
0.9326
-16.5127
6.1991
0.1934
7.4211
4.8656
Present FEM
21
22
Validation : Buckling Analysis
Validation : Frame Structure 1.2
1.01 1
1
0.99 0.8
P
T / Treference
Tx / Treference
0.98 0.6 0.4 0.2
Present FEM
0.97 0.96
Nastran, Ty
0.95
Nastran, Tz
0.94
Nastran, Rx
0.93
Nastran, Ry
5/16 L
0 Present FEM -0.2
0.92
Nastran, Tx
-0.4 2 4 6 Number of Elements of Beam 26
8
1/4 L
Nastran, Rz
0.91 0
7/16 L
Tapered beam EI = EIo (1+r x/L) r= 8
0
2 4 6 Number of Elements of Beam 26
8
23
NASTRAN Present FEM
P = 21.405 (16 CBEAM elements) P = 21.40492 ( 3 tapered beam elements) 24
PART II: Unsteady Aerodynamic Load
Integration Scheme Rodden’s New DLM scheme (1998) Control “r e c e i v i n g ”
Number of Integration Points
– Based on a boundary element method
5
5
5
x
5
5
– Two methods are developed to predict the unsteady aerodynamic load of multiple lifting surfaces:
5
5
5
5
5
5 5
5
• Doublet lattice method (DLM1) is similar to the works by Rodden et al. (1969, 1972, 1998), but with different doublet integration schemes to improve the accuracy and computational time
point
panel
5
5
5
Present Doublet-Lattice Method (DLM1) For one integration point, the present
• Doublet point method (DPM1) is similar to the work by Dowell et al. (1982) but is extended to non-planar surfaces and with sweep angle correction
1
3
1
1
7
3
1
1
Doublet lattice
1
1
DLM becomes DPM
3
3
1
1
5
3
“ sending” panel
25
26
Incomplete Cylindrical Functions: Analytical Solution
Kernel formulation
For X 0, Ueda (1982) derived the analytical solution as follows
• The velocity and pressure difference equation:
X
⌠ ⌠ e −−i k u v++1 / 2 == B v ,r ++ i B v,i 2 ⌡ ⌡ r ++ u 2
w = ∫ ∫ K ( r , R, X , X 1 , k, M ∞ ) ∆ Cp dx dy e −−i k x0 r2
K == −−
B v ,r =
(1 ++ u 2 )n++1/2
∫∫
du == r 2n
)
∑ ( −1 ) U n
v , 2n
B v ,r =
∞
2
∑ ( −1 ) U n=0
v
n +v ( − k / 2) n 1 1 k + + ∑ − γ − ln ∑ ( 2 v − 1)! ! m =1 m m = n+1 2 m 2
n
v , 2n +1
−
e ik u
−−∞ ∞
(r 2 ++ u 2 )n++1/2
U v ,m≠ 2 v =
du == r 2n B n 27
∫
∫ K dy = ∫ K X
K singular dx = 2 cos( γ r − γ s ) e
∑
n=0
dx +
∫K
n
)
Wn
•
The incomplete cylindrical function is approximated by Laschka or Desmarais series
•
No separation between regular and singular kernel functions
•
The numerator of the mixed singular/regular kernel function is approximated as a quadratic function (v. 1969) or a quartic function (v. 1998)
dx
regular
Regular Part: Gauss-Legendre Integration − i k x0
N
∫
K regular dx =
∑
n=1
Wn e
(K
1 ,reg
2n
−
2v −1
( kr ) U (m − 2 v ) m v ,m −2 28
Comparison of the methods
Singular Part : Analytical integration
singular
( kr / 2 )
n =0
2
k ( kX ) ( m − 2 v) m X 1
DLM of Rodden et al (1969, 1998)
(− ik tan Λ n!
∞
∑ n ! (n + v)!
Modification to the series above is used in the present method
Exact Separation of Singular and Regular Kernel Functions − i k x1 ∞
v
π ( − k / 2) 2 ( 2 v − 1)! !
m −1
X
e ik u
∫∫
−−X/r
(
2
∞
n=0
e −−i k x 0 M ∞∞r 4 i k X â 2 2R −− M ∞∞ X 2 e i k M ∞∞ ++ R ++ ++ 3I 2 T 2 X1 r2 R X 1
∞ ∞
I n ==
−−∞ ∞
M ∞∞r 2 i k X e ++ I 1 T 1 R X 1
T1 − K 2 ,reg T2
)
29
Present DLM1 • The incomplete cylindrical function is solved analytically using Ueda’s expansion series •
Exact separation between singular and regular kernel functions
•
The singular function is integrated analytically
•
The regular function is integrated using Gaussian quadrature 30
AGARD Wing
Nissim & Lottati Rodden Giesing Kalman
Panel # 4 x 4 (16) 6 x 6 (36) 8 x 8 (64) 6 x 3 (18) 8 x 4 (32) 10 x 5 (50) 11 x 12 (132) 14 x 12 (168) 18 x 12 (216) 22 x 12 (264)
C L H / ik 3.875 + 2.797 4.142 + 2.892 4.321 + 2.957 4.409 + 2.881 4.401 + 2.856 4.415 + 2.861 3.961 + 2.964 4.159 + 2.973 4.297 + 2.958 4.357 + 2.954
y 0.06
i i i i i i i i i i
Real
horizontal tail
Experiment Present DPM
0.04
Present DLM
x 0.03
S w e p t- w i n g w i t h p a r t i a l f l a p
AR = 2.94, TR=1, = 250 The flap is oscillated with
0.02
0.01
= 0.66 0 and k = 0.372 The wind tunnel data is due to Forsching et al. (AGARD CP8 0- 71, 1971)
0
C LH = C L due to heaving motions of the wing and tail k = b / U = reduced frequency = 1.5 M = 0.8
AGARD wing -
Imag .
0.05
∆ Cp ∆
Methods Present DLM
Wing with Flap
-0.01 0.4
31
0.5
0.6
0.7
0.8
0.9
1
32
Chordwise x at Section y = 0.87
YF-17 Wing Planform
Delta Wing, AR = 2 1.4
Methods
Paneling
CL
Present VLM (+ Gauss elimination)
7 x 3 (i/b) 4 x 15 (o/b)
2.944419
Present VLM (+ LU decomposition)
7 x 3 (i/b) 4 x 15 (o/b)
2.944418
Nissim-Lottati *
3 x (3 x 3)
2.953
Nissim-Lottati * (+ reverse flow theorem)
3 x (3 x 3)
2.941
V L M o f Giesing et al *
144
3.004
1.2
cl(y) c(y)/ CL cavg
(0.78 , 0.318)
Delta Wing
1 (1.1616 , 0.318)
0.8
(1.6346 , 1.26)
0.6
0.4
• Aspect ratio = 2
Present DPM
• Taper ratio = 0
• Angle of attack = 4.3
(2.007 , 0.0)
Present VLM
0.2
• Mach number = 0.13
(2.0066, 1.26)
VLM (NASTRAN)
0
0
*
0
0.2
0.4 0.6 spanwise y/b
0.8
Lottati, I. And Nissim, E., “3D Oscillatory Piecewise Continuou s
1
Kernel Function Method”, J. Aircraft, May 1981
33
34
UCAV Flexible Wing Twist Actuation*
Variable -Sweep Wing Planform
Lifting surface pressure distribution Wing tip twist angle ±10.5 ° NASTRAN = 1.57 in Present method =1.56 in
Paneling
CL
Present VLM (+ Gauss elimination)
5 x 4 (i/b) 2 x 16 (o/b)
3.019929
Present VLM (+ LU decomposition)
5 x 4 (i/b) 2 x 16 (o/b)
cp
c L = 0.48 0
- 0. 5
Z
Methods
c L = 0.28
-1 - 1. 5 0 2
3.019928
4 6
8
(1.53352, 0.49827)
c L = 0.0
0
1.38408)
3.02485
N/A
3.153
Lan (Quasi VLM)
4 x 28
3.007
- 0. 5 -1 - 1. 5 0
2 4
2
0
6
4 6 Y
8 8
10 12
10
- 0. 5
Z
(2.26086, 1.38408)
5 x 32
Truckenbrodt (Kernel Function)
Z
(2.04494,
(2., 0.)
Margason-Lamar (VLM)
X
14
12
16
14
-1
18 - 1. 5 0
*Taken from Kapania, R.K., Gern , F, …, “Current Research on Morphing Wing Technology,” AFDC Meeting, November 2000
35
2 4 6
8
36
Present SBW Wing 0.08
y=0.9875 y=0.9125 y=0.8125 y=0.6375
z/c
0.06 0.04
Comparison between Lamdes and Present VLM for the cambered wing, M = 0.85
M
0.8
0.02
g
0.7
0
0.6
0.2
0.4
0.6
0.8
0.16 0.14 0.12 0.1
Cl c(y) / cave
0
1
y=0.0125 y=0.1125 y=0.2625 y=0.3875
K, K g, M
K- P Kg, M
Present VLM (V6), v+w
0
0.6
0.8
Trim at 1.0 g cruise
Lamdes,v+w
0.1
0.2
0.4 0.6 y/half-span
0.8
1
f,
[ A(k) ]
Unsteady aerodynamic
0
0.4
FEM
[ A ]
Steady aerodynamic
0.3
0 0.2
EI, GJ, EA
0.4
Present VLM (V7), u+v+w
0
Hexagonal box model
0.5
0.2
0.08 0.06 0.04 0.02
t/c
double plate model
-0.02
z/c
PART III : Aeroelasticity With Compressive Force Effect u K P K q A u 0
mode shape
flutter V f,
1
37
38
x/c
SBW Flutter Speed with Compressive Force Effect
Validation of the Present Aeroelastic Code 0.5
Case 2
Case 3
Mach 0.678
Mach 0.90
Mach 0.95
Vf (fps) Wind Tunnel Test
f (Hz)
759.1
Vf (fps)
17.98
973.4
f (Hz) 16.09
V f(fps) 1008.4
f (Hz) 14.50
13.0 0.3 Frequency (Hz)
Case 1
Damping Coefficient g
21.966
Method
15.0
0.4
AGARD Standard 445.6 Wing*
0.2 0.1 0.0 -0.1 -0.2 -0.3
14.496
758.0
19.98
938.1
16.27
1158.2
15.79
ZONA6 (Linear)172
766
19.81
984
16.31
1192
16.18
ZTAIC (Nonlinear) 172
761
19.30
965.2
16.38
944.0
13.46
CAPTSD (Nonlinear) 172
768
*) 172
952
15.8
956
5.0
1.0 500
1000
1500
0
2000
200
12.8
400
600
800
1000
True Air Speed (fps)
True Air Speed (fps)
Present DLM 19.2
7.0
3.0
0
Method 30
9.0
-0.4 -0.5
Present DLM (Linear) *)
11.0
Structure Model
Flutter Speed (fps)
Flutter Frequency (Hz)
Beam
834.6
3.859
ZONA6
Beam
849.0
3.106
ZONA6
Detail*
873.2
3.109
NASTRAN
Detail*
902.0
3.93
LMAS
Detail*
1005
N/A
Yates, E.C., AGARD Report 765 )) Z ZA AE ER RO O T Th he eo o rr e e tt ii cc a a ll M Ma an nu ua a ll V V .. 5 5 .. 0 0 ,, Z ZO ON NA A T Te e cc h h .. ,, 2 20 00 00 0
39
40
*No compressive force effect
Variation of the Location of Strut Junctions
SBW Flutter Envelope 45
• Flutter speed of SBW
Wing without Strut Wing with Strut
40
– Structural damping is assumed to be zero – Calculated using the K method
35 1.6
1.8
LMAS Result
1.6
1.2
1.4
1 0.8 0.6
30
y / b = 0.69
Altitude (k ft)
1.4
V flutter/ Vref
Vflutter / Vref
y / b = 0.56
y / b = 1.0
1.2 1 0.8
Flutter Boundary
25 20
• No flutter occurs within the flight envelope
15
0.6
0.4 0.2
10
0.4
Wing-Strut with P [ Kg ] = 0
0.2
Wing-Strut with P [ Kg ] = 0
0
5
0 0
0.2
0.4 0.6 y / half span
0.8
1
0
1
2 3 4 Offset beam length h o f f s e t / h ref
5 41
0 0
0.25
0.5
0.75 1 1.25 Vflutter / Vref
1.5
1.75
2 42
Conclusion (1)
Conclusion (2)
• Aeroelastic analysis procedure for a strutbraced wing has been developed:
• Application of the procedure for aeroelastic analysis of the SBW with fuselage enginemounted configuration
– Accurate non -uniform beam finite element formulation – Improved doublet-lattice/doublet point unsteady aerodynamic methods – Compressive force effect on aeroelastic analysis procedure
– The compressive force tends to reduce the flutter speed – No flutter nor buckling occurs within the flight envelope of the present SBW – The effect of compressive force is significant near the wing tip 43
44
Conclusion (3)
Appendix
– The best location of the wing-strut junction is between 50% -80% of the wing span with the strut elastic axis attaches to the wing front spar – Flutter speed increases as the strut junction changes as shown:
45
SBW Lift Distribution for 2.5- g Maneuver
2
46
Optimum Configuration for SBW Fuselage-Mounted Engine
1.8 1.6 y
CL c / c ave
1.4 1.2 1
•
0.8 0.6 Rigid Wing
0.4
Flexible Wing with Strut
0.2
The wing with strut gives the lift distribution almost the same as the wing without strut, except near the wing root
0.2
0.4 0.6 0.8 Nondimensional half span
Strut chord (constant)
6.62 ft
Wing - strut junction point
74.52 ft
Wing root t/c
13.75%
Wing sweep (3/4 chord) Strut sweep (3/4 chord)
25.98 19.01
Wing tip t/c Strut t/c
Strut offset length
2.74 ft
Fuselage diameter
Wing root chord Wing tip chord
Flexible Wing without Strut
0 0
Wing half-span
1
47
108.44 ft
32.31 ft 6.77 ft
Wing flap area Wing reference area
6.44% 8.00% 20.33 ft 1411.02 ft2 4237.30 ft2 48
Computational Aeroelasticity
The First 6 Mode Shapes
• Flutter speed of the detailed model was calculated using NASTRAN • Assume 0 % structural damping
Z
Z
Y X
Y X
• The lowest flutter speed V
Mode 2 ( 2.21 Hz)
Mode 1 ( 1.72 Hz)
f
= 902 fps and frequency f = 3.93 Hz
0.6
14
Frequency (Hz)
Damping ratio g
0.2
Mode 4 ( 4.46 Hz)
Mode 3 ( 3.22 Hz)
Mode 2
12
Mode 3
Y X
Y X
Mode 1
Flutter Speed
0.4 Z Z
0 -0.2 -0.4
10
Mode 4 Mode 5
8
Mode 6 Mode 7
6
Mode 8 4
-0.6
Mode 9 Mode 10
Z Y X
Z
Mode 5 ( 6.93 Hz)
2
-0.8
Y X
-1
Mode 6 ( 7.83 Hz)
0 0
500
49
1000
1500
0
500
V (fps)
1000
1500
2000
50
V (fps)
Variations of the Spanwise and Chordwise Locations
Variation of Spanwise Locations of Strut Junctions
1.8 1.4
1.6 1.4 Vflutter / Vref
1.2
Reference Configuration 0.6
0.4
1.2 1 0.8 0.6
The strut support position
0.4
Front spar Elastic axis
0.2
Original Position
Rear spar
0
0.2
0
0.2
0.4
0
0.6
0.8
1
y / halfspan 0
0.2
0.4
0.6
0.8
1
51
52
y / halfspan
Variations of the Strut Root Location along Fuselage
Variations of the Offset Beam Length 1.4
1.6
1.2
1.5
1
Reference configuration
1.4
V /V flutter ref
Vflutter / V
0.8
Vflutter / Vref
ref
1
1.3 1.2 Original Position
1.1
0.8
0.6
0.4
1 0.2
Leading edge wing root
0.9
0
0.8
0
-60
-40 -20 0 20 40 60 80 Strut junction location at the fuselage (ft)
100
0.5
1
1.5
2 h o f f s e t/ h
53
2.5
3
3.5
4
offset-Ref.
54
Flutter Boundary 45
Wing+Strut, zero fuel Isolated Wing, full fuel Isolated Wing, zero fuel
35
3
ft)
30 25
Altitude (10
• Flutter speed of the SBW for
Wing+Strut, full fuel
40
– Failure case : strut is inactive
boundary 1.2 V
15
• During the positive-g maneuver the strut is stretched in tension.
– Nominal case
Flutter
20
SBW Buckling Analysis
d
• Calculated using the PK Method of NASTRAN
• Because the angle between strut and wing is shallow, the wing may subject to a high compressive load • For the present SBW configuration,
• No flutter occurs within the flight envelope
• Need to consider the change of the wing-strut stiffness due to the compressive force effect K mod = K - P B
– the strut force Fs
10 5 0 0
0.25
0.5
0.75
1
Vflutter
1.25
1.5
/ Vref
1.75
2
55
Dissertation : Motivation • Need consideration of the compressive force effect on the aeroelastic analysis – NASTRAN does not have flutter and buckling analyses in the same package – Indirect procedure using DMAP procedure in NASTRAN is complicated
800 kips during the + 2.5 g maneuver
56
Research Objective • To develop an aeroelastic analysis procedure with the compressive force effect of the strut-braced wing – To develop a FEM based code which c o m b i n e s b u c k l i n g a n d aeroelastic stability analyses
• The aeroelastic flutter analysis is performed off line
– T o e n h a n c e t h e f o r m u l a t i o n o f n o n -u n i f o r m beam finite elements to increase the accuracy and reduce computational time
– Challenge to combine NASTRAN and the MDO tool
– To improve vortex/doublet lattice methods
57
• Enhancement of the structural analysis
Application to SBW : Inner-Wing Buckling Analysis
Proposed Scheme M
u
C
u
K
P Kg
q A
u
10
F
Wing without strut Wing with Strut
9
M
P Kg u
u
K
Buckling
0
P Kg
8
u
Modal Analysis
0
K
P Kg
qc A
u
F
K
P Kg
qd
A
u
0
K
P Kg
Passive Load Alleviation Divergence
Buckling Load P / Prequired
K
y
7 6
u
C
u
q A
u
0
Wing buckling analysis for the 2.5 g maneuver case
•
The buckling load increases as the wing- strut junction moves inboard
Original position
4 3 2
Flutter
•
5
1
M
58
for more accurate prediction and fast
0 59
0
0.2
0.4 0.6 The junction location y / L
0.8
1
60
The Effect of Offset Length and Junction Position
The Effect of Offset Beam Position 8
Spanwiseposition of the junction 7
h
y = 0.65
6
y
y = 0.76
5
y = 0.80 4
y = 0.84
The stability margin increases with the decreased of the s p a n w i s e position of the strut junction y
2 1
y = 0.65
2
1
The junction moves inboard
Y
0 0
2
4
6
8
10
61
/ h a
y
The response surface of the buckling load as functions of the
0
lfs
pa
y = 0.84
3
3
y = 0.90
y = 0.97
strut offset length and position y = 0.97
y = 0.90 Original height
y = 0.71
4
y = 0.76
5
P / P required for 2.5 g maneuver
h
y = 0.71 6
y = 0.80
P / Preference for 2.5 g maneuver
7
se Off
t he
t h igh
62
n
Offset length h / hreference
Overview
Further Work • The present procedure can be used to analyze the aeroelastic stability of the SBW with underwing engine configuration. • For the SBW with wing -tip configuration, additional procedure should be added to take into account the effect of engine thrust force oscillation. • To accurately captured the transonic-dip effect in the flutter analysis, a higher order of unsteady aerodynamic load prediction should be used 63
Early Aeroelastic SBW : Keldysh SBW Problem
y
– Aeroelastic Sensitivity Analysis with respect to several design variables
• Conclusion
• Region ys= 0.47 L • Vdiv = V flutter
64
• Neglect buckling problem No geometric stiffness effect
Stable
• Region ys < 0.47 L • Vdiv > V flutter
region
ys
• V flutter is max at y s = 0.47 L
• Optimum location y s > 0.47 L. The strut increases Vcritical by 100%
• Neglect strut mass effect
ys = s p a n w i s e
• Present SBW Aeroelastic Analysis
Flutter Speed
• U n s w e p t rectangular high AR wing Used a simple strip theory approach • Neglect strut aerodynamic model No aerodynamic interference
ys
– Double Plate Model – Equivalent Hexagonal Box Model
• V div is constant
• Uniformly distributed wing stiffness Simplified the governing equations
ys x
– Keldysh SBW Problem
• Present SBW Finite Element Model
• Region ys > 0.47 L • Vdiv < V flutter
Imposed a simple support at wing-s t r u t j u n c t i o n
U
– MDO of Present Strut -Braced Wing (SBW) Aircraft
• Early Aeroelastic Analysis of SBW Wing Problem
Mailybaev-Seiranyan Solution
• I n f i n i t e l y -r i g i d s t r u t s t i f f n e s s z
• Introduction
ys ys = s p a n w i s e p o s i t i o n o f t h e s t r u t j u n c t i o n
position of the strut junction
65
66
DISKUSI
ADE JAMAL Penyederhanaan Model Elemen untuk Non uniform beam element method bertujuan untuk mempercepat waktu hitung. Apakah clustering computing bisa menjadi alternatif menggunakan Elemen Hingga Konvensional? ERWIN SULAEMAN Clustering computing tentu saja bisa dipakai untuk metoda elemen hingga konvensional. Pengaruh dari clustering computing bisa signifikan untuk pemecahan masalah elemen hingga yang ukuran degrees of freedomnya besar.
DODDY KASTANYA 1. 2.
Kalau kecepatan cruising dibatasi lebih rendah dari yang dipakai secara komersial saat ini, bagaimana menjual ide baru ini ke perusahaan pembuat pesawat? Berapa flutter speed dari pesawat komersial saat ini?
ERWIN SULAEMAN 1.
2.
Ide menggunakan strut untuk menopang sayap pesawat bisa saja digunakan untuk pesawat kecepatan rendah. Bahkan kemungkinan besar lebih baik karena pengaruh gaya drag dari strut tidak serumit dibandingkan dengan kecepatan tinggi. Kecepatan cruise biasanya diinginkan setinggi mungkin. Dengan penggunaan strut yang tepat posis inya, kecepatan crukenya bisa dinaikkan dibandingkan dengan sayap tanpa strut. Untuk tipe pesawat sejenis B777, biasanya sekitar Mcruise = 0,85. Kecepatan flutternya tergantung desain dari pesawat tersebut akan menurut peraturan harus lebih tinggi dari 1,15 MD .
BAKRI ARBIE Saya mengusulkan agar Batan, UI dan ITB mengadakan pertemuan reguler yang penting untuk masa depan bangsa. ERWIN SULAEMAN Tentang kerjasama, sebenarnya dalam PT DI lebih banyak dalam perancangan. Untuk bidang design sangat senang diundang untuk berdiskusi degan Batan, BPPT ataupun Universitas.
DAFTAR RIWAYAT HIDUP 1. Nama
: Erwin Sulaeman
2. Tempat/Tanggal Lahir
: Bandung, 31 Januari 1959
3. Instansi
: PT. Dirgantara Indonesia
4. Pekerjaan / Jabatan
: Head of Structure Analysis Group, NMX
5. Riwayat Pendidikan
:
• 1986, S1-ITB • 1993, S2-Univ. of Dayton • 2001, S3-Virginia Polytechnic Institute & State University 6. Pengalaman Kerja
:
• 1986-sekarang, PT Dirgantara Indonesia • 2001-2003, Zona Technology 7. Organisasi Profesional
: AIAA
8. Makalah yang Pernah disajikan : • Whirl Flutter • LCO (Limit Cycle Oscillation) • Strut Braced Wing Aeroelasticity • Boundary Element Method
Back
