524 IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 1 AEROELASTIC FLUTTER ANALYSIS OF SUPERSONIC WING WITH MULTIPLE EXTERNAL STORES NUR AZAM AND ERWIN SULAEMAN Department of Mechanical Engineering, Faculty of Engineering, International Islamic University Malaysia, P.O. Box 10, 50728, Kuala Lumpur, Malaysia. esulaeman@iium.edu.my ABSTRACT: Flutter may be considered to be one of the most dangerous aeroelastic failure phenomenon. The flutter characteristic differs for each aircraft type, and depends on the wing geometry as well as its operational region of subsonic, transonic or supersonic speeds. Prior to performing a flight flutter test, extensive numerical simulations and Ground Vibration Tests should be conducted where the structural finite element modes and the experimentation results should be matched, otherwise the numerical simulation model must be rejected. In this paper, the analysis of simulation of a supersonic wing equipped with external missiles loaded on the wing is presented. The structural mode shapes at each generated frequency are also presented. The analysis is carried out using MSC Nastran FEM software. The wing flutter with the external stores was simulated at different altitudes. The result shows that the flutter velocity is sensitive to the flight altitude. For this reason, the flutter analysis is conducted also for a negative altitude. The negative altitude is obtained by considering the constant equivalent speed- Mach number rule at the flutter speed boundary which is a requirement in standard transport aircraft regulations. ABSTRAK: Salah satu fenomena kegagalan aeroelastik yang paling membahayakan adalah kipasan (flutter). Ciri-ciri kegagalan kipasan (flutter) adalah berbeza untuk setiap jenis pesawat bergantung pada geometri sayap dan regim operasi sama ada subsonik, transonic atau supersonik. Sebelum melakukan ujian penerbangan kipasan ,simulasi berangka luas dan ujian getaran peringkat bawahan (darat) perlu dijalankan di mana struktur mod unsure terhingga dan keputusan eksperimen harus dipadankan, sebaliknya model simulasi berangka boleh ditolak. Dalam kertas kerja ini, simulasi sayap supersonic dilengkapi dengan beban luaran peluru berpandu di sayap telah dianalisis di daerah supersonic tinggi. Bentuk mod struktur pada setiap mod frekuensi yang dihasilkan juga dipersembahkan secara visual. Analisis ini dilakukan dengan menggunakan perisian “FEM” MSC Nastran. Kepakan sayap dengan kedai-kedai luar telah disimulasikan pada ketinggian yang berbeza. Hasil kajian menunjukkan bahawah alajukipasan(flutter) sensitive terhadap ketinggian penerbangan. Atas sebabini, analisis flutter dilakukan juga untuk ketinggian negatif. Ketinggian negative diperoleh dengan mempertimbangkan tetap bersamaan kelajuan Mach beberapa peraturan pada kelajuan sempadan debar sebagai keperluan dalam peraturan piawaiaan pesawat pengangkutan. KEYWORDS: aeroelasticity; flutter; finite element method; external store; mode shape; supersonic wing 1. INTRODUCTION Aeroelasticity can be defined as a branch of aeromechanics that deals with the interaction among inertial, aerodynamic and structural stiffness effects in air vehicle IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 2 design [1]. Flutter is a type of aeroelastic instability in which the structure extracts energy from the air stream and unstable self-excited oscillation results in catastrophic structural failure. There are several types of mathematical models used to represent the aircraft structure and it’s mass in aeroelasticity. These include the “Beam Like” model, “Box Like’ model, and“Box Like Condensed to a Beam Like” model [2]. In this paper, the model used for the wing structure is the “Box Like” representation where structural detail is explored and analyzed accordingly. The present flutter analysis is conducted based on certification standards for military fighter aircraft such as MIL-A-8861[3] which put the aeroelastic requirement under Subpart 3: Construction, Material and Design and MIL-A-8870 at Section 6.4 [4]. To determine the initial sizing of wing structure, a maximum aerodynamic load distribution during pitch up maneuver is assumed and a finite element approach has been used to set up the structural stiffness for the wing by recognizing its detailed ‘box like’ construction. The present wing structure data should be validated further by comparing to Ground Vibration Test (GVT) results where the main frequencies and mode shapes of the numerical model should be in agreement with the GVT results [5]. External stores can be defined as any object such as missiles, gun, and fuel tank which are attached to the wing outside the wing box model. Since the stores are attached to the wing, the missile launcher is designed based on the missile’s weight and there should be a clearance of at least 8 cm if the missiles are mounted near to each other [6]. Each external store center of gravity must be placed as far forward as possible with respect to the elastic axis to delay the occurrence of flutter, while satisfying other design constraints [7]. For this paper, the analysis is performed at varying altitudes where the flow is supersonic and the external stores attached to the wing are two types of missiles. The technical data for the missiles were taken from [8] and [9]. The aeroelastic simulation model used in the present work consists of two parts: structural model and aerodynamic model. The structural model of the wing box skin is modeled using quadrilateral shell elements, while the external store and launcher are modeled using bar elements. The analyses take into account the in-plane membrane characteristics as well as lateral bending and shear deformation of the shell element. The material used for the wing box is Aluminum which is used as a base line data for future work where the wing box skin is designed as a composite structure. The aerodynamic model of the wing for the supersonic region is based on the boundary element method where the unsteady aerodynamics is modeled using the ZONA51 method of MSC Nastran [10]. For the subsonic region, the doublet lattice method of MSC Nastran is used [12]. The wing lifting surface is divided into chordwise and spanwise directions using a number of trapezoidal panel elements. It is noted that even though the present wing is designed for the supersonic region, the flutter analysis is required to be conductedinthe subsonic region as well to ensure the instability does not occur in the complete flight design envelope. To ensure compatibility between the structural finite element and aerodynamic boundary element model, a surface spline is used such that the 6 degrees of freedom structural deformations at each point is related to each aerodynamic trapezoidal panel inclination angle. There are several methods to analyze the flutter speed based on frequency matching such as the K method and the PK method [10]. This paper presents the flutter speed results using the PK Method. Based on [2], the PKMethod is used to determine the aerodynamic stiffness and damping matrices as a function of reduced frequency. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 3 The SOL 103 module of MSC Nastran is used to simulate the normal modes of each frequency. The SOL 145 module is developed to obtain the damping and reduced frequency variation as a function of velocity. From this latter module, the flutter speed can be determined when the graph of velocity versus damping factor is plotted. This paper presents the first 10 fundamental normal modes of SOL 103 and determines the flutter speed of this supersonic wing with external stores by considering the variation of altitudes. 2. STRUCTURAL ANALYSIS 2.1 Supersonic Wing Characteristics The present work utilizes a wing platform with an aspect ratio of 5 and taper ratio of 0.5. The wing leading edge swept back angle is 30°. The airfoil of this supersonic wing is a double wedge shape as shown in Fig. 1. The wedge angle of the airfoil is 10°. Along the wing span, the airfoil is divided into three parts which are the main wing box and two control surfaces at the leading and trailing edge. The portions of the leading and trailing edge have been specified as 15% and 20% of the chord length, respectively. The performance of the selected airfoil uses the characteristics provided by [11] for higher supersonic region analysis. The present wing design is used as a baseline for further work where the wing geometry as well as wing composite structure is set as the sensitivity parameter to obtain an optimum supersonic wing design. Fig. 1: Double wedge airfoil. For the external store, Fig. 2 shows the configuration of the loaded missile on the wing. The external stores for each station of wing are specified in Table 1. There are two types of missiles used which are AMRAAM and Sidewinder. Fig. 2: External stores configuration of the wing. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 4 Table 1: External stores technical data. Station Missile Type Length [m] Diameter [m] Mass [kg] 1 AIM-9M 2.85 0.128 86.0 2 AIM-120 A 3.66 0.178 157.89 3 AIM-120 A 3.66 0.178 157.89 2.2 Wing Loading Based on [3], the load factor for the fighter aircraft is set atnz = 5.5. The load for this wing, as shown in Fig. 3, is assumed to be an elliptic load acting along the wing span wise (y-axis) direction and symmetric quadratic load along chord wise (x-axis) direction. With this load assumption, the sizing of the wing box can be conducted. Fig. 3: Wing loading estimation equation. The formula to calculate the load factor is given by Eq. (1) in which L is the lift and W is the weight of one side of the aircraft wing based on [11]; nz � LW (1) Here the lift can be calculated using Eq. (1) L � n�W (2) The span wise elliptic load can be formulated as: � y a� 2 �z bb � 2 � 1 (3) Where the value of parameter a is half the span length since it is the length of the major axis, and parameter b is the minor axis.The value of b can be calculated using (6). The chord wise quadratic load distribution is given by: � � �� ��� ���� (4) IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 5 The area of the quadratic load in chord wise direction can be calculated by integrating Eq. (4) acting along the x axis. Aquadratic � � zb �1‐ �xe�2# e ‐e dx (5) Then, the volume of the elliptic load can be found by integrating Eq. (5) along the y axis in Eq. (3) V� 43 � zb'y(e'y(dy a 0 (6) To find the minor axis of the elliptic equation, Eq. (3), equation (2) is divided by 2 since this is only applicable for the half wing, equal to the volume found in Eq. (3). This expression can be written as: V� L2 (7) The wing can be assumed as a beam along y axis to find the shear force Q and moment M of each section as denoted in Eq. (8) and Eq. (9), respectively. Q'y(� � dQ (8) M'y(� � dM (9) 2.3 Wing Sizing To calculate the skin thickness of the overall wing, the wing is divided into 4 regions along the span and the thickness in each region is calculated based on the maximum load in the respective region as shown in Fig. 4 and Fig. 5. Fig. 4: Top view of skin thickness division region. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 6 Fig. 5: Skin thickness segment in a region. The moment of inertia formula is given as: I� � h2 dA (10) The moment of inertia for the front and rear spar are calculated as a vertical segment in (11): Ixx� 112 th3 (11) The moment of inertia for the inclined segments which has an inclination angle of θ, can be derived using Eq. (10). This can be done by setting the limit for integration along z axis starting from 0 to the end of each inclination segment denoted as./. The final formula is given by Eq. (12): I � tcosθ � ' h2 xtanθ(2dz b' 0 (12) By assuming the thickness to be constant at every skin and spar, equation (12) reduces to form an equation to find the thickness in any region based on the associated moment Macting in that region as shown in Eq. (13): σy � M h middleI't( (13) where 7 is calculated based onEq. (9), 89:;;<= is the height of the inclination for the front spar of the wing only and >'?( is the summation moment of inertia of the wing box in terms of t as given in Eq. (11) and Eq. (12). 2.4 Safety Factor The safety factor FS for the structural strength analysis in the present work is set at 1.5. The shear stress at any location in any region in the y direction is denoted as @A and is calculated using Eq. (14) with Q as the shear force in that region calculated using Eq. (8). The FS can be computed using Eq. (15) or Eq. (16). τy� QSb I (14) IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 7 τallowτy E F. S (15) σallowσy E F. S (16) 3. UNSTEADY AERODYNAMICS For the prediction of flutter, unsteady aerodynamicloadacting on the wing surface, which is oscillating according to the structural dynamic mode shapes, is estimated using the boundary element method. For the subsonic region, the unsteady aerodynamic loads are calculated using the doublet lattice method (DLM). While for the supersonic region, a constant pressure method of ZONA51 is used. For both methods, the wing is modeled as a flat lifting surface and is discretized into a number of trapezoidal elements. Following [10], a set of aerodynamic influence coefficients in the form of a matrix equation is generated. The fundamental relationships between the lifting pressure and the dimensionless normal velocity induced by the inclination of the surface to the air stream can be formulated as in [10]: HwjJ�KAjjLMfj/qP (17) where w is the normal wash velocity, f is the aerodynamic pressure, q is the dynamic pressure and A is the aerodynamic influence coefficient matrix. The substantial differentiation matrix of the structural deflections to obtain the downwash is given by Eq.(18): HwjJ�KDjk1 i k Djk2 LMukP MwjgP (18) where k is the reduced frequency. The integration of the pressure to obtain forces and moments yields, MPkP�KSkjLMffP (19) The three matrices in Eqs. (17), (18) and (19) can be combined to give the aerodynamic influence coefficient matrix in Eq. (20): UQkkV�KSkjLKAjjL‐1UDjk1 ikDjk2 V (20) The ZONA51 and DLM theories compute the A matrix. Then, the matrix decomposition forward and backward substitutions are used in the computation of the Q matrix. Since this wing will be operating in the supersonic region, the ZONA51 of Nastran had been used. For this part, the Nastran coding development in view of aerodynamics will consider the outer part of the wing box structure including the control surface as shown in Fig. 1. 4. PK METHOD OF FLUTTER SOLUTION Following [10] the PK equation for modal flutter analysis can be formulated as in Eq. (21) IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 8 WMhh p2 YBhh‐ 14 ρc\VQhhIk ] p ^khh‐ 12 ρV2QhhR `a MuhP � 0 (21) where the circular frequency ω and the reduced frequency k are related to p as k� ω c / 2V (22) p � ω ' 2 g i( (23) The flutter solution is rewritten in the state space form as in Eq. (24) where A is complex . UA ‐ p IVMuhP� 0 (24) The eigen solution of Eq. (24) is in the form of a complex eigen value p for each mode, which in turn will give the structural damping g for the real part and frequency ω for the imaginary part. Note that the result is computed for each velocity which is embedded in the damping and stiffness matrix terms of Eqn. (21). 5. RESULTS 5.1 Normal Modes – SOL 103 of MSC Nastran For the FEM data, the boundary condition at the wing root is rigidly fixed, i.e. no deflection and no rotation in the x, y,and z directions at the front, middle and rear spars of the wing box. The first 8 normal modes of the wing structure are shown in Fig. 6 and Fig. 7. Note that the rigid body mode is not included in the list. The frequency and its associated shape are recorded in Table 2. Fig. 6: Normal mode 1 to 4. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 9 Fig. 7: Normal mode 5 to 8. Table 2: Flutter velocity and flutter Mach number at Mode 5. Altitude [ft] Flutter Speed [m/s] Flutter Mach number Flutter Frequency [Hz] -7943 518.5 1.48 9.289 0 581.5 1.71 9.282 10000 677.5 2.06 9.272 20000 799.5 2.53 9.261 30000 954.5 3.15 9.249 5.2 Flutter Solution – SOL 145 of MSC Nastran SOL 145 simulation is further carried out at different Mach numbers to find the match point velocity at 0 ft (sea level). The graph of structural damping versus velocity of every mode at 0 ft is plotted in Fig. 8. Besides that, a graph of velocity versus frequency at 0 ft in Fig. 9 is also plotted to show the rapid changes of frequency of the flutter mode. The flutter dominant mode most likely occurs at Mode 5 since the graph in Fig. 8 shows that the structural damping becomes zero when it approaches a velocity of 581.5 m/s. This is the match point for the graph since at this speed it gives the result of Mach number 1.71. The simulation for this solution is repeated for different altitudes. The flutter velocity and flutter Mach number for this variation is shown in Table 2. The result shows that the present wing flutter is sensitive to altitude. For this reason, the calculation is performed also for a negative altitude where hneg = - 7,943 ft. This negative altitude is derived analytically in [13] for transport aircraft and military UAV (Unmanned Air Vehicle) [14] as the result of all combinations of altitudes and speeds encompassed by the VDive or MDive versus altitude envelope enlarged at all points by an increase of 15 percent in equivalent airspeed at both constant Mach number and constant altitude. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 10 Fig. 8: V-gplot at 0 ft. Fig. 7: V-f plot at 0 ft. Fig .8: V-gplot at various altitudes. IIUM Engineering Journal, Vol. 15, No. 2, 2014 Azam and Sulaeman 11 6. CONCLUSION The present work establishes a procedure for the design of wings based on flutter analysis. The results show that the external stores have a significant effect on wing flutter. It is found that flutter speed of the present wing-external store configuration is sensitive to the flight altitude. The present result is used as a baseline for further work where aeroelastic tailoring is performed to find an optimum composite structure for the main wing box in order to find higher flutter speed and lighter wing structure. ACKNOWLEDGEMENT Special thanks to the Science and Technology Research Institute for Defense (STRIDE), Malaysia for allowing the authors to use MSC Nastran and MSC Patran in their facilities. The present research is funded by a research grant from the Ministry of Education Malaysia (Grant no. FRGS13-020-0261). REFERENCES [1] Erwin HJ.(1997)MSC Developments in Aeroelasticity.The MacNeal-Schwendler CorporationAerospace Users' Conference Proceedings (California, USA): 1-9. [2] Wright JR, Cooper JE. (2007)Aircraft Aeroelasticity and Loads. John Wiley & Sons Ltd. Manchester, United Kingdom. [3] MIL-A-M8861B(A) (1986) Military Specification :Airplane Strength and Rigidity Flight Loads. Lakehurst, New Jersey: Naval Air Engineering Center, Systems Engineering and Standardization. [4] MIL-A-8870C(AS) (1993)Military Specification: Airplane Strength and Rigidity Vibration, Flutter, Divergence. Lakehurst, New Jersey: Commanding Officer, Naval Air Warfare Center Aircraft Division. [5] Peeters B, Climent H, Diego RD, Alba JD, Ahlquist JR, Carreño JM, Wim H, Alain R, Gonzalo G, Jeroen D, Jan D.(2008)Modern solutions for Ground Vibration Testing of large aircraft. Proceeding of 26 th International Modal Analysis Conference (IMAC 26Orlando, USA) 1(4):110-123. [6] Raymer DP. (1999)Aircraft Design: A Conceptual Approach. Reston: American Institute of Aeronautics and Astronautics, Inc. [7] Janardhan S, Grandhi RV, Eastep F, Sanders B. (2005) Parametric Studies of Transonic Aeroelastic Effects of an Aircraft Wing/Tip Store.Journal of Aircraft42:253-263. [8] Approved Navy Trainning System Plan For The AIM-120 Advanced Medium Range Air to Air Missile. Nevada: Naval Strike and Air Warfare Center (NSAWC), 1988. [9] AIM -9 Sidewinder. Retrieved on May 2013 from Operation Guide: www.vojenskeletectvi.cz [10] Willian PR, Erwin HJ. (2004)MSC.Nastran version 68:Aeroelastic Analysis User’s Guide. MSC.Software Corporation, California USA. [11] Motte ED (2012) Gripen: When logic is part of the equation. Farnborough, UK: Gripen Export, Saab Aeronautics. [12] YuKH, Djojodihardjo H, Kadarman H. (2011) Acoustic effects on binary aeroelasticity model. IIUM Engineering Journal 12(2):123-126. [13] Sulaeman E. (2012) On the evaluation of negative altitude requirement for flutter speed boundary of transport aircraft and UAV. Applied Mechanics and Materials225:397-202 [14] UAV Systems Airworthiness Requirements (USAR) for North Atlantic Treaty Organization (NATO) Military UAV Systems – Standard Agreement STANAG 4671 (Edition 1), 2009.