AIAA 2008-6368

AIAA Modeling and Simulation Technologies Conference and Exhibit 18 - 21 August 2008, Honolulu, Hawaii

AIAA Modeling and Simulation Technical Conference 18 Aug 2008, Honolulu, Hawaii

AIAA-2008-6368

Development of a Pilot Training Platform for UAVs Using a 6DOF Nonlinear Model with Flight Test Validation Hou In, Leong * and Rylan Jager † The University of Kansas, Lawrence, KS, 66045 Shahriar Keshmiri ‡ The University of Kansas, Lawrence, KS, 66045 and Dr. Richard Colgren § The University of Kansas, Lawrence, KS, 66045

This paper describes the evaluation of three different modeling and simulation techniques to support unmanned aerial vehicle development for glacial ice research. A Cloud Cap Technologies Piccolo II autopilot is integrated within a 1/3 scale Yak-54. Two hardware-in-the-loop (HIL) simulation models of the Yak-54 are developed using tools provided with the Piccolo II. Conventional parametric modeling is utilized separate from the Piccolo system and two state space linear models are developed for both longitudinal and lateral dynamics. Open loop flight tests are conducted and the results are used to evaluate the accuracy of each modeling method. The best modeling data set is selected and used to develop a six-degree of freedom (6DOF) nonlinear model. A comparison of the nonlinear and linear model’s responses to the flight test data for each dynamic mode was performed. The study reveals that the 6DOF nonlinear model provides more accurate estimates of the dynamics of the vehicle than the linear model in every aspect. A pilot training platform is developed using the 6DOF nonlinear model in conjunction with a pilot joystick device and 3D visualization software for real-time simulation.

G

I. Introduction

lobal sea level rise has become one of the major environmental concerns of the world today. As the atmosphere warms one of the most significant contributors to the rise of the Earth’s oceans is the melting of the great ice sheets of Greenland and Antarctica. Current understanding of the behavior of the great ice sheets is limited as quality glacial data sets are not yet available. To meet this need to gather more data on the great ice sheets the National Science Foundation established the Center for Remote Sensing of Ice Sheets (CReSIS) (http://www.cresis.ku.edu) in 2005. CReSIS is a multi-university science and technology center whose main goal is to conduct and foster multi-disciplinary research that will result in technologies, new datasets, and models necessary to achieve a better understanding of the mass balance of the polar ice sheets (e.g., Greenland and Antarctica) and their contributions to sea level rise. In order to drive sensing technologies CReSIS has proposed the use of an Unmanned Aerial Vehicle (UAV) to carry in-house developed radar systems for the purpose of glacial ice depth measurement and analysis. This UAV, designed by the Aerospace Engineering Department at the University of Kansas and named the Meridian, is of medium size with a 26 ft. wingspan and weighing approximately 1,000 lbs. The Meridian will be equipped with an off-the-shelf autopilot as its primary guidance, navigation and control system. A through test and evaluation of the Piccolo II autopilot system, manufactured by Cloud Cap Technologies, *

Graduate Research Assistant, Department of Aerospace Engineering, [email protected], AIAA Student Member. Graduate Research Assistant, Department of Aerospace Engineering, [email protected], AIAA Student Member. ‡ Assistant Professor, Department of Aerospace Engineering, [email protected], AIAA Member. § Associate Professor, Department of Aerospace Engineering, [email protected], Fellow of the AIAA. 1 American Institute of Aeronautics and Astronautics †

020208 Copyright © 2008 by Copyright@2008_by_HIL. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

has been conducted. The Piccolo system is provided with a hardware-in-the-loop (HIL) simulation platform, called the Cloud Cap Simulator in this paper. Two parametric based modeling approaches, named Standard Cloud Cap Simulation (SCCS) model and Athena Vortex Lattice (AVL) model, are provided with the Cloud Cap Simulator to compute the aerodynamic derivatives. These modeling methods allow users to rapidly develop a vehicle dynamic model and simulate its response so that the performance of the system can be evaluated. In addition, the simulation is used to pre-tune the control system gains prior to first flight testing. Simulation also allows for a concrete gain tuning procedure to be developed by the users, as follow up gain tuning in flight will most likely be necessary. Because of the size and speed of the Meridian an accurate simulation that will minimize in-flight gain tuning is crucial. To evaluate the Piccolo flight control and simulation systems CReSIS use small scale R/C aircraft. Much research literature is available on the subject of modeling and simulation of UAVs but few address the details of modeling the aerodynamic derivatives for small UAVs. Wind tunnel testing1 can be utilized to determine the derivatives but it is a labor intensive task with high cost. Another approach is to use the geometric parameters of the vehicle to estimate the derivatives through computational methods. This approach has been applied in several UAV research programs.2,3 Ref. 4 conducted research similar to the work as presented in this paper, however, only the longitudinal dynamics were studied. This paper describes a complete six degree of freedom (6DOF) nonlinear simulation model developed using MATLAB/Simulink. The aerodynamic derivatives are obtained from the geometric data of the vehicle. The open loop dynamics are first calculated using methods laid out in Roskam5 and compared with the results of the Piccolo HIL simulator. Data reduction techniques6 are introduced and used to measure the damping ratio and frequency of the various dynamic modes from the time response data. Flight test data are presented and compared with each modeling method. From the comparison result, the best derivative data sets are selected and used to develop a 6DOF nonlinear model. Side-by-side comparison techniques are introduced and used to validate the accuracy of the 6DOF nonlinear and state space liner models using the flight test data.

II. Standard Cloud Cap Simulation (SCCS) Model The Piccolo II autopilot system is first integrated on 1/3 scale Yak54 RC (Figure 1) aircraft to validate and demonstrate the functionality and integrity of the flight control system. A pitot and static tube is also integrated with this aircraft to collect altitude and airspeed data. This aircraft incorporates several different safety features including dual battery power supplies with individual power switches for both the Piccolo and control surface actuation servos and an emergency engine kill switch controlled by a separate 72 MHz Futaba transmitter and console. A. Physical Geometry The modeling method used on the SCCS is mainly based on the Figure 1. 1/3 Scale Yak-54 RC Model geometric parameters of the aircraft which includes the wing, ailerons, horizontal tail, elevator, vertical tail, rudder and fuselage. The data were directly measured from the Yak-54. Some of the major geometric information are listed in Table 1.

Wing Horizontal Tail Vertical Tail

Wing Span (ft) 7.93 3.00 1.42

Table 1. Basic Geometry of Yak-54 Reference Aspect Taper Incidence Dihedral Area (ft2) Ratio Ratio Angle (deg) Angle (deg) 10.90 5.8 0.45 0.0 0.0 2.29 3.91 0.81 0.0 0.0 1.59 1.26 0.35 0.0 0.0

Quarter Cord Sweep Angle (deg) -2.0 12.60 12.68

B. Aerodynamic Data Aerodynamic tables of the lifting surfaces are required to complete the dynamic model of the aircraft. Each component is defined by three values as a function of angle of attack; life coefficient (CL), drag coefficient (CD), and pitching moment coefficient (CM). These coefficients can be obtained from equations (1) and (2).7 2 American Institute of Aeronautics and Astronautics 020208

CL = CL α × α CD = CD0 +

(1) 2

CL π × AR × e

(2)

Because the airfoils are symmetric, the value of CM for each lifting surface is approximated as zero for all angles of attack. The method to compute the values of CL α and CD0 is discussed in Section III. C. Moment of Inertia The moments of inertia for the Yak-54 are approximated using a component build-up method. The aircraft model is first disassembled in small components, left and right wings, left and right tails, wing and horizontal tail spars, propeller, spinner, engine, batteries, and the Piccolo control unit. Then each component is weighed individually. The position of each component is measured relative to the engine firewall datum. Those data are then used to calculate the moment of inertia8 of the Yak-54 about the X-axis, Y-axis and Z-axis of the Yak-54 in the body coordinate system. The Yak-54 is a symmetric aircraft and the weight distribution on the left and right wings is almost symmetric, thus the moment of inertia of XZ plane is zero. The results of the calculation are summarized in Table 2. D. Hardware-In-the-Loop Simulation Setup To implement the hardware-in-the-loop simulator9 an external CAN interface is utilized to connect the Piccolo avionics to the PC. The PC simulates the aircraft sensor data and sends it to the avionics through the CAN bus. The avionics will then react to the sensor data by generating the actuator command data. The whole system will close the loop by feeding the actuator commands back to the simulator PC while applying them to the aircraft model, computing the new sensor data and sending them back to the avionics through the CAN bus interface. To visualize the simulation response the simulator PC is networked to another PC, which runs the Flight Gear simulator software.10 Figure 2 shows the HIL simulation environment.

Table 2. Component Weights of Yak-54 Aircraft Components Weight (lbs) Fuselage Section 10.96 Right Wing 2.38 Left Wing 2.32 Right Horizontal Wing 0.62 Left Horizontal Wing 0.58 Wing Spar 0.43 Rear Spar for Horizontal Tail 0.08 Front Spar for Horizontal Tail 0.04 Propeller 0.42 Spinner 0.51 Engine 5.76 Servo Battery #1 0.44 Servo Battery #2 0.44 Piccolo Battery #1 0.32 Piccolo Battery #2 0.32 Ignition Battery 0.18 Piccolo Autopilot Unit 1.04 Total Weight 26.84 Moment of Inertia (slug-ft2) IXX (body axis) 1.0886 IYY (body axis) 2.1068 IZZ (body axis) 3.0382

Figure 2. Piccolo Hardware in the Loop Simulation Setup9 Once the required geometry and aerodynamic data are input by the user, the simulator software will compute the dynamics of the aircraft and simulate its response. In addition to visualizing the response on Flight Gear, the response is also plotted on a real time strip chart interface11 supported by the Piccolo system. One of the major disadvantages of the SCCS is the extreme black box nature of this simulator. It is very difficult to troubleshoot any mistakes in the modeling as no derivatives are output by the simulator. For this reason, mathematical analysis techniques cannot be applied for cross checking. The best way to evaluate the validity of the 3 American Institute of Aeronautics and Astronautics 020208

simulator therefore is to perform mode identification flight tests6. This is done by exciting the different modes of the aircraft using singlet or doublet control surfaces inputs. The response is then analyzed using data reduction methods6 to estimate the dynamic characteristics of each mode. E. Data Reduction Methods The data reduction methods described herein were developed by flight test engineers at the University of Kansas, based on standard industry practices. As simple “hand-calculation” type techniques to estimate frequency and damping. The details of these methods can be found in Ref. 6. Table 3. Applicability of Different Methods6 There are five different basic data reduction Range of Applicable Data Reduction Method techniques that are discussed in Ref. 6. Each of them is Damping Ratio designed for specific situations. All of them give an Transient Peak Ratio (TPR) -0.5 < ξ < 0.5 estimate of the damping ratio and natural frequency for a second order system, which is broken down from a Time Ratio (TR) 0.5 < ξ < 1.2 quadratic characteristic equation for the state space Maximum Slope (MS) 0.5 < ξ < 1.2 model. Only three of the methods are used in this article. The choice of method is mainly driven by the damping ratio of the response to be analyzed. Table 3. summarizes the applicability of these three methods for different ranges of damping ratio. F. SCCS HIL Simulation Results The simulation tests were performed for four different modes; Dutch Roll mode, roll mode, short period mode and Phugoid mode. The simulation responses are presented in Figure 3. The dynamics of each mode are measured from the time response traces using data reduction methods. For the Dutch Roll and Phugoid mode, the TPR method is used as they are lowly damped modes. For the short period mode the MS method is applied. Since the roll mode is a first order mode, finding the time constant (τroll) is the most important dynamic characteristic for this mode. According to the definition of time constant12 it is the time to reach 63.2% of the final constant value when the input is activated. Table 4 summarizes the dynamic analyses for each mode. Roll Mode Response - Simulation Test Date : Thu Dec 06 2007 Roll Rate P (deg/s) 1106

1106.5

1107

1107.5

1108

1108.5

1109

15

30 20 10 0 1295.5

1296

1296.5

1297

1297.5

1298

1298.5

1299

1297.5

1298

1298.5

1299

1297.5

1298

1298.5

1299

40 20 0 -20 -40 1295

10

Rudder (deg)

40

Time (sec)

Time (sec)

1295.5

1296

1296.5

1297

Time (sec) 2

LHS Aileron (deg)

5 0 -5 -10 -15

50

-10 1295

Bank Angle φ (deg)

Yaw Rate (deg/s)

Dutch Roll Mode - Date : Simulation Test Date: Dec 06 2007 100 80 60 40 20 0 -20 -40 -60 -80 -100 -120 -140

1106

1106.5

1107

1107.5

1108

1108.5

1.5 1 0.5 0 1295

1109

1295.5

1296

1296.5

1297

Time (sec)

Time (sec)

Phugoid Mode - Simulation Test Date : Thu Dec 06 2007

Short Period Response - Simulation Test Date : Thu Dec 06 2007 30

True Air Speed (Knots)

Pitch Rate (deg/s)

80 20 10 0

75

70

-10

65 -20 974

974.5

975

975.5

976

976.5

977

977.5

800

978

810

820

830

840

850

860

840

850

860

Time (sec)

Time (min)

-3.6

-2

-3.8

Elevator (deg)

Elevator δE (deg)

-4 -6 -8 -10 -12 974

-4 -4.2 -4.4 -4.6 -4.8

974.5

975

975.5

976

976.5

977

977.5

978

-5 800

810

Time (min)

820

830

Time (sec)

Figure 3. SCCS HIL Simulation Responses 4 American Institute of Aeronautics and Astronautics 020208

Table 4. Summary of Dynamics Analysis for SCCS HIL Simulation Responses Damping Natural Frequency Time Constant Period Mode T (sec) τroll (sec) Ratio ξ ωn (rad/sec) Roll ---------0.30 Dutch Roll 0.15 9.08 0.69 ---Short Period 0.67 12.67 0.50 ---Phugoid 0.16 0.36 17.45 ----

III. Advanced Aircraft Analysis (AAA) Model The Advanced & Aircraft Analysis (AAA) software is utilized to compute the stability and control derivatives, which are used to compose the state space models. The AAA based simulation model is constructed in the MATLAB environment. The methods described in Roskam5 are used to compute the aircraft dynamic model. Two aircraft state space models are constructed to predict the longitudinal and lateral dynamic response. In this technique, the user can gain more insight into the modeling details and have an improved modeling ability. A. Stability Derivatives and Control Derivatives Advance & Aircraft Analysis (AAA)13 was developed by DAR Corporation.14 It is primarily used for preliminary aircraft design purposes. AAA has a built in aerodynamic database for different size aircraft, mainly based on DATCOM.15 For any given aircraft geometry and flight condition, AAA interpolates based on the input geometry and the aerodynamic database to estimate the stability and control derivatives. Although this method cannot replace any wind tunnel test or flight test, it provides a rapid method to give an initial analysis of aircraft performance and stability. The same geometric data are used here as was used in the SCCS model, as described before in Table 1. The flight condition is set to the straight and level flight trim condition. Trim speed and altitude are captured from previous flight test data, which is 1200 ft ASL (above sea level) and 0.106 Mach. The resultant derivatives are summarized in Table 5. Table 5. AAA Model Derivatives Steady State Derivatives

CTx

CD1

0.1470

0.0422

CDu

CDα

CTx

CLu

CLα

CLα

CLq

Cmu

0.0011

0.0863

-0.1546

0.0017

4.5465

1.8918

5.5046

0.0002

Cmα

Cmα

Cmq

CmT

CmT

CDδ e

CLδ e

Cmδ e

-0.3937

-4.3787

-8.0532

0.00

0.0275

0.00

0.3792

-0.8778

1

Cm1

CmT

CL1

0.0515 0.0001 Longitudinal Derivatives (rad-1) u

u

1

0.0009

α

-1

Lateral Derivatives (rad )

C yβ

Cyp

C yr

Clβ

Cl p

Clr

Cnβ

CnT

-0.3602

0.0085

0.2507

-0.0266

-0.3819

0.0514

0.1022

-0.0045

Cn p

Cnr

C yδ a

C yδ r

Clδ a

Clδ r

Cnδ a

Cnδ r

-0.0173

-0.1270

0.00

0.1929

0.3490

0.0154

-0.0088

-0.0996

β

B. Longitudinal State Space Model The standard 6-DOF equations of motion for a conventional fixed wing aircraft can be found in many existing textbooks. Here, methods described in Roskam5 are applied to develop the dynamic model of Yak-54. The equations for the longitudinal dynamics are shown below. 5 American Institute of Aeronautics and Astronautics 020208

= − gθ cos Θ1 + X u u + X Tu u + X α α + X δ e δ e

u

U1α − U1θ = − gθ sin Θ1 + Z u u + Zα α + Zα α + Z qθ + Zδ e δ e

θ

(3)

= M u u + M Tu u + M α α + M Tα α + M α α + M qθ + M δ e δ e

For the straight and level flight condition, θ = q and θ = q . Rewriting Eq. (3) in state space format, yields:

0 ⎡1 ⎢0 (U − Z ) α 1 ⎢ ⎢0 − M α ⎢ 0 0 ⎣

0 0 ⎤ ⎡ u ⎤ ⎡ X u + X Tu ⎢ 0 0 ⎥⎥ ⎢⎢α ⎥⎥ ⎢ Zu = 1 0 ⎥ ⎢ q ⎥ ⎢ M u + M Tu ⎥⎢ ⎥ ⎢ 0 1 ⎦ ⎣θ ⎦ ⎣⎢ 0

Xα

0

Zα

U1 + Z q

M α + M Tα

Mq

0

1

− g cos Θ1 ⎤ ⎡ u ⎤ ⎡ X δ e ⎤ ⎢ ⎥ ⎥ − g sin Θ1 ⎥ ⎢⎢α ⎥⎥ ⎢ Zδe ⎥ δ + ⎥ ⎢ q ⎥ ⎢M ⎥ [ e ] 0 ⎥ ⎢ ⎥ ⎢ δe ⎥ 0 ⎦⎥ ⎣θ ⎦ ⎢⎣ 0 ⎥⎦

(4)

The definition of each variable used in Eq. (3) and Eq. (4) can be found in Ref. 5. The dimensional stability derivatives are calculated from the dimensionless derivatives, as listed in Table 5. The final state space model for the longitudinal dynamics used for simulation is presented in Eq. (5)

0.0000 −32.1554 ⎤ ⎡ u ⎤ ⎡ 0.0000 ⎤ ⎡ u ⎤ ⎡ −0.2374 12.4194 ⎢α ⎥ ⎢ −0.0042 −7.7904 0.9232 −0.0091 ⎥ ⎢α ⎥ ⎢ −0.6438 ⎥ ⎢ ⎥=⎢ ⎥⎢ ⎥+⎢ ⎥ [δ ] ⎢ q ⎥ ⎢ 0.0163 −19.3108 −9.0798 0.0299 ⎥ ⎢ q ⎥ ⎢ −105.5538⎥ e ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ 0.0000 1.0000 0.0000 ⎦ ⎣θ ⎦ ⎣ 0.0000 ⎦ ⎣θ ⎦ ⎣ 0.0000

(5)

By finding the eigenvalues of the dynamics matrix of the state space model, one can analyze the characteristics of the longitudinal dynamics. The result is shown in Table 6. Table 6. Longitudinal Dynamics Analysis from AAA Model Damping Frequency Period Eigenvalues Mode Ratio (rad/sec) (sec) -0.114 ± 0.248i 0.42 0.27 23.27 Phugoid -8.440 ± 4.118i 0.90 9.42 0.67 Short Period The analysis reveals that the Yak-54 model has two complex conjugated roots for its longitudinal dynamics, the Phugoid mode and short period mode. Both modes exhibit a highly damped response. This is especially true for the short period mode, which reacts as a first order response with no oscillation. This result will be validated by the flight test data and discussed in Section VI. C. Lateral State Space Model The perturbation equations of motion for the lateral dynamics are shown below:

U1β + U1ψ = gφ cos Θ1 + Yβ β + Ypφ + Yrψ + Yδ a δ a + Yδ r δ r

φ − A1ψ

= Lβ β + Lpφ + Lrψ + Lδ a δ a + Lδ r δ r

ψ − B1φ

= N β β + NTβ β + N pφ + N rψ + Nδ a δ a + Nδ r δ r

(6)

For the straight and level flight condition, p = φ ; p = φ and r = ψ ; r = ψ . Equation (7) shows the state space format of Eq. (6). The final numerical values of the lateral dynamics state space model are shown in Eq. (8). 6 American Institute of Aeronautics and Astronautics 020208

⎡ 1 ⎢ ⎢ 0 ⎢ 0 ⎢ ⎣ − B1

0 1

0 0 0 U1 0 0

− A1 ⎤ ⎡ p ⎤ ⎡ Lp ⎥⎢ ⎥ ⎢ 0 ⎥ ⎢ φ ⎥ ⎢ 1 = 0 ⎥ ⎢ β ⎥ ⎢ Yp ⎥⎢ ⎥ ⎢ 1 ⎦ ⎣ r ⎦ ⎢⎣ N p

⎡ p ⎤ ⎡ −16.6421 ⎢ φ ⎥ ⎢ 1.0000 ⎢ ⎥=⎢ ⎢ β ⎥ ⎢ 0.0005 ⎢ ⎥ ⎢ ⎣ r ⎦ ⎣ 0.0926

0 0 g cos Θ1 0

0.0000 −37.3608 0.0000 0.2722

0.0000 −0.6238

0.0000

46.4013

⎤ ⎡ p ⎤ ⎡ Lδ a ⎢ 0 0 ⎥⎥ ⎢ φ ⎥ ⎢ 0 ⎢ ⎥+ Yβ (Yr − U1 ) ⎥ ⎢ β ⎥ ⎢ Yδ a ⎢ ⎥ N r ⎥⎦ ⎢⎣ r ⎥⎦ ⎢ Nδ ( N β + NTβ ) ⎣ a Lβ

Lr

Lδ r ⎤ ⎥ 0 ⎥ ⎡δ a ⎤ Yδ r ⎥ ⎢⎣δ r ⎥⎦ ⎥ Nδ r ⎥⎦

2.3631 ⎤ ⎡ p ⎤ ⎡ 454.0359 22.8520 ⎤ 0.0000 ⎥⎥ ⎢⎢ φ ⎥⎥ ⎢⎢ 0.0000 0.0000 ⎥⎥ ⎡δ a ⎤ + −0.9854 ⎥ ⎢ β ⎥ ⎢ 0.0000 0.3341 ⎥ ⎢⎣δ r ⎥⎦ ⎥ ⎥⎢ ⎥ ⎢ −2.0395⎦ ⎣ r ⎦ ⎣ −14.0214 −46.9710 ⎦

(7)

(8)

Table 7 summarizes the results of the eigenvalue analysis for the lateral state space model. Table 7. Lateral Dynamics Analysis from AAA Model Damping Frequency Period Time Constant Eigenvalues Mode Ratio (rad/sec) (sec) (sec) 0.0115 ---------86.96 Spiral -16.7 ---------0.06 Roll -1.32 ± 6.75i 0.19 6.88 0.91 ---Dutch Roll The analysis shows that the Yak-54 has an unstable spiral mode with a long time constant, which means this mode will diverge very slowly. The roll mode has a very small time constant so the Yak-54 is very agile in roll. For the Dutch Roll mode, it shows a low damped oscillatory motion with a period close to one second. These results will be verified with the flight test data and presented in Section VI.

IV. Athena Vortex Lattice (AVL) Model Cloud Cap Technology offers a second modeling technique that addresses some of the issues associated with the SCCS model. This simulator makes use of the Athena Vortex Lattice (AVL)16 software, a vortex lattice code developed by Mark Drela and Harold Youngren at the Massachusetts Institute of Technology. The basic logic behind AVL is that the aircraft stability derivatives generated by its lifting surfaces can be estimated based on the steady vortex shedding of the surfaces at small angles of attack and sideslip. In AVL the fuselage effects and engine power are ignored. In the model each major airfoil section can be defined. Using an aircraft’s wing as an example, the center airfoil section, the airfoil section where the flap or aileron starts and ends, the airfoil section where the sweep angle changes, or the airfoil section where the dihedral changes, can all be defined so the geometry will exactly match the actual aircraft. This also allows for the creation of butterfly tail geometry models directly, which is not an option in the SCCS. A significant advantage of AVL is that the geometry can be graphically represented in three dimensions, allowing for increased troubleshooting ability. Figure 4 shows the geometry model for the Yak-54 used in the vortex analysis. The stability derivatives about the center of gravity are then calculated using the lifting surface geometry. The angle of attack of the aircraft is varied based on intervals set by the user and the stability derivatives are determined for each angular position. The stability derivatives are output to an .XML file for use by the Cloud Cap Simulator. The simulator uses the same propulsion and inertia models as the SCCS, substituting in the AVL aerodynamics for the Figure 4. Yak-54 AVL Model lifting surface tables. The .XML file allows greater insight into what goes into the Cloud Cap Simulator, allowing for better troubleshooting of the dynamics model. Ref. 17 contains a more detailed discussion of the AVL simulator. Table 8 shows the longitudinal and lateral derivatives generated by the vortex lattice simulation. The AVL model is trimmed at an angle of attack of α=1.0 degree. All the derivatives shown in Table 8 are given from this trim condition and are compared with the AAA model derivatives. 7 American Institute of Aeronautics and Astronautics 020208

Table 8. AVL Model Derivatives Steady State Derivatives AAA AVL (% Difference)

CL1

CD1

0.1470 0.1426 (-2.99)

0.0422 0.03245 (-23.10)

CTx

Cm1

CmT

0.0515 N/A

0.0001 0.00974 (9640)

0.0009 N/A

1

1

Longitudinal Derivatives (rad-1) AAA AVL (% Difference)

AAA AVL (% Difference)

CDu

CDα

CTx

CLu

CLα

CLα

CLq

Cmu

0.0011 N/A

0.0863 N/A

-0.1546 N/A

0.0017 N/A

4.5465 3.7007 (-18.60)

1.8918 N/A

5.5046 3.5927 (-34.73)

0.0002 N/A

Cmα

Cmα

Cmq

CmT

CmT

CDδ e

CLδ e

Cmδ e

-0.3937 -0.2850 (-27.61)

-4.3787 N/A

-8.0532 -4.3732 (-45.70)

0.00 N/A

0.0275 N/A

0.00 N/A

0.3792 0.3967 (4.61)

-0.8778 -0.7572 (-13.73)

u

α

u

Lateral Derivatives (rad-1) AAA AVL (% Difference)

AAA AVL (% Difference)

CnT

C yβ

Cyp

C yr

Clβ

Cl p

Clr

Cnβ

-0.3602 -0.2727 (-24.29)

0.0085 0.0194 (128.24)

0.2507 0.2531 (0.96)

-0.0266 -0.0314 (18.05)

-0.3819 -0.5858 (53.39)

0.0514 0.0743 (44.49)

0.1022 0.1052 (2.94)

-0.0045 N/A

Cn p

Cnr

C yδ a

C yδ r

Clδ a

Clδ r

Cnδ a

Cnδ r

-0.0173 -0.0387 (123.93)

-0.1270 -0.1156 (-8.98)

0.00 0.00 (0.0)

0.1929 0.2228 (15.50)

0.3490 0.3707 (6.22)

0.0154 0.0219 (42.47)

-0.0088 N/A

-0.0996 -0.1003 (0.70)

β

Due to the lack of a provision for an engine model within the AVL method, all thrust derivatives are not available. The AVL model also does not provide the derivatives related to drag and the rate of change of angle of attack. Because of this an accurate state space model cannot be constructed. However, the Phugoid and short period approximations (Ref. 5) were used to estimate the frequency and damping of each mode using the available stability derivatives from Table 8. For the lateral dynamics, the derivatives CnTβ and Cnδa are assumed to be zero. The lateral state space model is developed using the same techniques as described in Section III for eigenvalues mode analysis. As described before, the AVL model derivatives can be used on the Cloud Cap simulation platform to conduct hardware-in-the-loop simulation. Using the same simulation technique as presented in Section II, the dynamics of different modes are excited by a singlet or doublet input. The Dutch roll, short period and Phugoid modes responses conducted in the AVL HIL simulation are shown in Figure 5. Data reduction methods are then applied to extract the data of interest from the responses to identify the characteristics of each mode. 20

80

50

0 -20

75

-40 -60 468

468.5

469

469.5

470

470.5

0

-50

224

471

True Air Speed (Knots)

40

-80 467.5

100

q - Pitch Rate (deg/s)

Yaw Rate (deg/s)

60

224.5

225

225.5

226

226.5

227

65

60 370

Time (sec)

Time (sec)

70

380

390

400

410

420

430

440

420

430

440

Time (sec) 0.9

6

10

0.8

0

-5

-10 467.5

Elevator (deg)

Elevator (deg)

Rudder (deg)

4

5

2 0 -2

468

468.5

469

469.5

470

Time (sec)

Dutch Roll Mode

470.5

471

0.6 0.5 0.4 0.3

-4 -6 224

0.7

224.5

225

225.5

226

226.5

227

0.2 370

Time (sec)

Short Period Mode

380

390

400

410

Time (sec)

Phugoid Mode

Figure 5. AVL HIL Simulation Responses The results given from the approximation method, state space model eigenvalues analysis and the AVL HIL simulation responses are all summarized in Table 9. 8 American Institute of Aeronautics and Astronautics 020208

Table 9. AVL Model Mode Analysis and HIL Simulation Results Spiral

Modeling Method

Roll

τs

τr

(sec) 55.87 ----

AVL AVL-HIL

Dutch Roll

(sec) 0.04 ----

Short Period

ωn (rad/sec) 7.13 7.06

ζ 0.16 0.15

ζ 0.65 0.85

Phugoid

ωn (rad/sec) 7.50 12.89

ζ 0.15 0.17

ωn (rad/sec) 0.37 0.30

V. Flight Testing Results As the objective of this test is to identify the characteristics of the open-loop aircraft dynamics, the response of each dynamic mode needs to be excited. The test method used in flight testing is same as was used in the SCCS and AVL simulation test. As this is an actual flight test, disturbances exists due to atmospheric conditions. To ensure the quality of the flight test results, each test subject is conducted more than once to ensure sufficient data are available for analysis so that the uncertainties are minimized. The Dutch Roll mode test was conducted five times and three sets of data were chosen for analyses as shown below. The TPR method is used to estimate the characteristic of the mode. The final estimation is the average of the results from those three responses and is listed in Table 10. Dutch Roll Mode - Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

60

60 X: 1122.1 Y : 53.92

Yaw Rate (deg/s)

40

30

X: 1082.45 Y : 19.26

20 X: 1083.5 Y : 4.647

10 0 -10

X: 1082.9 Y : -9.586

-20

X: 1129.5 Y : 37.16

40

20

Yaw Rate (deg/s)

40 Yaw Rate (deg/s)

D utch R oll Mode - Flight Test D ate : Thu Dec 20 2007; Time : 13-29-35

Dutch R oll M ode - Flight Test D ate : Thu Dec 20 2007; Time : 13-29-35

X: 1081.5 Y : 41.79

50

X: 1123.05 Y : 1.736

0 X: 1123.75 Y : -5.649

-20 X: 1122.6 Y : -20.11

-40

X: 1128.6 Y : 47.21

20

X: 1130.8 Y : 4.515

0 -20

X: 1130.15 Y : -12.36

-40

-30 -60

-40 -50 1080.5

1081

1081.5

X: 1081.9 Y : -43.17

1082

1082.5

1083

1083.5

1084

1084.5

1085

1085.5

-80 1121

1086

-60

X: 1121.7 Y : -74.98

1121.5

1122

1122.5

Time (sec)

1123.5

1124

1124.5

-80 1128

1125

1128.5

1129

1128.5

1129

X: 1129.1 Y : -73.57

1129.5

T ime (sec)

8

8

7

6

6

4 3 2

2 0 -2 -4

1

1130.5

1131

1131.5

1132

1130

1130.5

1131

1131.5

1132

5 Rudder (deg)

5

1130 Time (sec)

10

4 Rudder (deg)

Rudder (deg)

1123

0

-5

-6

0 1080.5

1081

1081.5

1082

1082.5

1083

1083.5

1084

1084.5

1085

1085.5

-8 1121

1086

1121.5

1122

1122.5

Time (sec)

1123

1123.5

1124

1124.5

-10 1128

1125

1129.5

Time (sec)

Tim e (sec)

Figure 6. Dutch Roll Mode Flight Test Responses Table 10. Dutch Roll Mode Flight Test Analysis Results Test I Test II Test III Average 0.22 0.28 0.32 0.27 Damping Ratio ξ 6.04 5.95 5.10 5.70 Natural Frequency ωn (rad/sec)

Three sets of flight data were selected to analyze the dynamics of the roll mode. The responses are shown as follows. The time constant is measured from the response and is listed in Table 11. Roll Mode Response - Test Date : Thu Dec 20 2007; Time : 13-29-35

Roll Mode Response - Test Date : Thu Dec 20 2007; Time : 13-29-35

-20 X: 1300.25 Y : -46.08

-40 1299

1299.5

1300

1300.5 Time (sec)

1301

1301.5

X: 1358.1 Y : -42.06

1357.5

1358

1358.5

1359

1299

1299.5

1300

1300.5 Time (sec)

1301

1301.5

1302

-0.5 -1 -1.5 1300

1300.5 Time (sec)

Bank Angle φ (deg)

20 0

0.5

1301

1301.5

1302

1357.5

1358

1358.5

1359

1563

1563.5

1562

1562.5 Time (sec)

1563

1563.5

1562

1562.5

1563

1563.5

-20 -40 -60 -80 -100 1561.5

1359.5

X: 1357.45 Y: 0.2636

1

0 -0.5 -1 -1.5 -2 1357

X: 1562.2 Y : -81.23

1562.5 Time (sec)

Time (sec)

0

1299.5

1562

0

40

-20 1357

X: 1298.9 Y : 0.3209

1299

-50

-100 1561.5

1359.5

LHS Aileron (deg)

LHS Aileron (deg)

Bank Angle φ (deg)

-40

LHS Aileron (deg)

Bank Angle φ (deg)

0 -20

-2 1298.5

-40

60

20

0.5

-20

0

Time (sec)

40

-60 1298.5

0

-60 1357

1302

Roll Rate P (deg/s)

Roll Rate P (deg/s)

Roll Rate P (deg/s)

0

-60 1298.5

Roll Mode Response - Test Date : Thu Dec 20 2007; Time : 14-05-38 50

20

20

1357.5

1358

1358.5

1359

1359.5

X: 1561.6 Y : 0.3495

0 -1 -2 -3 -4 1561.5

Time (sec)

Figure 7. Roll Mode Flight Test Responses

9 American Institute of Aeronautics and Astronautics 020208

Table 11. Roll Mode Flight Test Analysis Results Test I Test II Test III -42.06 -81.23 Constant Roll Rate (deg/sec) -46.08 Time to Reach 63.2% 1299.57 1357.70 1561.99 of Constant Roll Rate (sec) Time for Initial Input (sec) 1298.85 1357.45 1561.60 0.72 0.25 0.39 Time Constant (sec) Three out of five sets of flight test data were chosen from the short period test. Note that no overshoot is observed from the flight test response. The TR method is used here to measure the damping ratio and natural frequency of this mode. The results are summarized in Table 12. Short Period Response - Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

Short Period Response - Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

60

Short Period Response - Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

50 40

20 0 -20

30 Pitch Rate (deg/s)

Pitch Rate (deg/s)

Pitch Rate (deg/s)

40

20 10 0 -10 -20 -30

-40

-40 -60 1210

1210.5

1211

X: 1211.2 Y : -52.62

1211.5

1212

1212.5

1213

1213.5

-50 1214

1214

1214.5

1215

X: 1220.55

X: 1215.2 Y : -46.11

1215.5

Time (sec)

1216

1216.5

1217

1217.5

1218

1219

2

-6 -8

0

Elevator (deg)

Elevator (deg)

Elevator (deg)

-4

-2 -4 -6

1212

1212.5

1213

1213.5

1221.5

1222

1222.5

1223

1221.5

1222

1222.5

1223

0 -2 -4 -6 -8

-8 1214

1214

1221

2

2

1211.5

1220.5Y : -51.64

4

0 -2

1211

1220

Time (sec)

4

1210.5

1219.5

Time (sec)

4

-10 1210

50 40 30 20 10 0 -10 -20 -30 -40 -50

1214.5

1215

1215.5

1216

1216.5

1217

1217.5

1218

Time (sec)

-10 1219

1219.5

1220

1220.5

1221 Time (sec)

Figure 8. Short Period Flight Test Responses Table 12. Short Period Flight Test Analysis Results Test I Test II Test III Average 0.97 1.00 1.03 1.00 Damping Ratio ξ 16.44 16.62 Natural Frequency ωn (rad/sec) 13.47 19.96 Two Phugoid mode tests were performed in the flight test. The mode was excited by an elevator input to gain altitude until the airspeed dropped about 10 knots while the throttle is maintained constant. As the plot reveals the test was not completed as the Phugoid mode has a very long period. Before the mode started to be excited, the pilot had to take over the control of the aircraft due to the limited visual range. The flight test results are shown in Figure 9. Phugoid M ode - Flight Test Date : Thu Dec 20 2007; Time : 14-05-38 80

70

70

True Air Speed (Knots)

True Air Speed (Knots)

Phugoid Mode - Flight Test Date : Thu Dec 20 2007; Time : 14-05-38 80

60 50 40 30 1449

1450

1451

1452

1453 1454 Time (sec)

1455

1456

1457

40 30 1496

500

1500

1502

400 300 1450

1451

1452

1453 1454 Time (sec)

1455

1456

1457

1508

1510

1504

1506

1508

1510

1504

1506

1508

1510

Time (sec)

200 1498

1500

1502

2 Elevator (deg)

4 2 0 -2 -4 -6 1449

1506

300

100 1496

1458

1504

400

6 Elevator (deg)

1498

500 Altitude AGL (ft)

Altitude AGL (ft)

50

1458

600

200 1449

60

1450

1451

1452

1453

1454

1455

1456

1457

1458

Time (sec) 0

-2

-4 1496

1498

1500

1502 Time (sec)

Time (sec)

Figure 9. Phugoid Mode Flight Test Responses The first and second tests lasted for about five and ten seconds respectively before the pilot reassumed control. If the estimation from the AAA result is correct, the period of the Phugoid mode response is approximately 26 10 American Institute of Aeronautics and Astronautics 020208

seconds. In order to observe the first peak of the oscillation response it takes half of the period, which is at least 13 seconds. Therefore, it is reasonable to anticipate that the actual period of the Phugoid mode is longer than 20 seconds and the flight test conducted did not last long enough to see the response due to the concerns of pilot controllable visual range.

VI. Comparison of Flight Test and Simulation Results Now the analysis of the Yak-54 open loop dynamics from the flight test data is completed. It would be useful to compare the actual flight dynamics results with the estimations predicted from different modeling and simulation methods. The summary of comparison is shown as below. Table 13. Comparison of Flight Test with Simulation Results Lateral Dynamics Longitudinal Dynamics Modeling Spiral Roll Dutch Roll Short Period Phugoid Method ωn ωn ωn τs τr ζ ζ ζ (rad/sec) (rad/sec) (rad/sec) (sec) (sec) 86.96 0.06 0.42 0.27 0.19 6.88 0.90 9.42 AAA ------0.15 9.08 0.67 12.67 0.16 0.36 SCCS-HIL 55.87 0.04 0.16 7.13 0.65 7.50 0.15 0.37 AVL ------0.15 7.06 0.85 12.89 0.17 0.30 AVL-HIL ------0.27 5.70 1.00 16.62 ------Flight Test

A. Comparison of the Dutch Roll Mode From the Dutch roll mode dynamics comparison, it reveals that all the simulation results under-estimate the damping ratio term slightly. The errors among the different modeling techniques are similar. For both the damping ratio and natural frequency, the AAA model provides the closest results to the flight test data among all the modeling methods. It is concluded that the AAA method gives the best estimation for the Dutch roll mode dynamics. B. Comparison of the Short Period Mode For the short period mode, Table 13 indicates that the AAA also provides the best estimation to the damping ratio data (≈10% difference), while the AVL-HIL gives the closest predication to the natural frequency term (≈30% difference). Keep in mind that according to the flight test results, the aircraft demonstrated very highly damped short period response with no oscillation. Therefore, the natural frequency data is not reliable and thus should not be used to evaluate the results. It is also seen that the AVL-HIL model shows a large difference in the short period mode as opposed to using the short period approximation in the AVL aerodynamics only model. This is most likely due to some power effects in the Cloud Cap simulator. According to the comparison of the damping ratio data, it is concluded that the AAA model also gives the best predication of the short period mode dynamics. C. Comparison of the Phugoid Mode There were, however, large differences between the models in Phugoid mode damping. The frequency predicted by the AAA model and the AVL-HIL were very close but the damping ratio of the AAA model Phugoid mode is a full 0.25 higher. Because of the incomplete flight data, direct comparison with the Phugoid mode could not be made. However, plots of the various modal predictions were a useful tool in evaluating the Phugoid dynamics. For this reason, MATLAB/Simulink was used to simulate the AAA state space model presented in Section III. Using similar elevator inputs to those used in flight test as shown in Figure 9, a 10 second time period plot was made and is compared with the AVL-HIL Phugoid mode simulation result as shown in Figure 10. The damped airspeed response from the AAA model was slightly faster than the results from the flight test, with a half period peak occurring about 6 seconds into the response. It should also be noted that with these large airspeed changes the validity of the model tends to degrade. The AVL-HIL simulation shows a faster initial response with less damping. In any event, the Phugoid response in flight test appeared to be slower than either simulator. This indicates that the actual dynamics fit closer to the higher damped Phugoid predicted by the AAA model. 11 American Institute of Aeronautics and Astronautics 020208

True Air Speed (Knots)

True Airspeed (Knots)

70

70 60 50 40

0

1

2

3

4

5

6

7

8

9

60

50

40

30

10

401

402

403

404

405

Time (sec)

406

407

408

409

410

407

408

409

410

Time (sec) 2

Elevator (deg)

0.5

1.5

Elevator (deg)

0 -0.5 -1 -1.5 -2

1 0.5 0 -0.5 -1

0

1

2

3

4

5

6

7

8

9

-1.5

10

401

402

403

404

405

406

Time (sec)

Time (sec)

AVL-HIL Simulation Responses

AAA MATLAB/Simulink Simulation Responses

Figure 10. Comparison of Phugoid Mode Dynamics Simulation Responses D. Comparison of the Roll Mode The time constant for the roll mode should not be compared directly. The time constant calculated from eigenvalue analysis refers to a singlet aileron step input, which is an instant input with no time delay. However, the actual pilot inputs shown from the flight test takes about 0.5 seconds to establish a constant aileron input, which is 8 to 12 times higher than the calculated time constant given from AAA and AVL. The roll time constant for the Yak-54 is very small, making it highly affected by the elapsed time of the aileron input. The only possible way to truly evaluate the results is to make a side-by-side comparison using identical aileron input. This can only be done on the MATLAB/Simulink platform, and only the AAA and AVL state space models can be implemented within MATLAB/Simulink for simulation. The pilot commands used in the roll mode test are first extracted from the flight test data. These commands are then placed into both the AAA and AVL state space models in MATLAB/Simulink to simulate the responses. The simulation results are then compared with the flight test responses so a side-by-side comparison can be achieved based on identical inputs. Note that the trim values for the control surfaces are essential as they will affect the dynamics responses of the aircraft from the trim condition. The pilot commands shown in the flight test data reflect the actual trim values set by the pilot. In the simulation, the aircraft state space model was trimmed at predefined condition. Therefore any control commands placed on the simulator are considered as perturbation signals. To compensate, the trim value used for simulation input needs to be reset. The final comparison results for roll mode test I, II and III are shown in Figure 11. Roll Mode Simulation Response - AAA Model vs AVL Model Compare with Flight Test Response. Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

Roll Mode Simulation Response - AAA Model vs AVL Model Compare with Flight Test Response. Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

40

-40 -60 -80 0.5

1

1.5

2

2.5

3

3.5

4

30 20 10 0 -10 -20

4.5

0.5

1

Flight Test Data AAA State Space Model AVL State Space Model

Roll Rate - p (deg/sec)

-40

1.5

2

-80 -100

3

0

0.5

1

1.5

2.5

3

3.5

4

0

Flight Test Data AAA State Space Model AVL State Space Model

-20 -30 -40 -50

4.5

Aileron (deg)

Flight Test Data Simulation Input

3.5

-80 -100 -120

0

0.5

1

1.5

2

2.5

Flight Test Data AAA State Space Model AVL State Space Model

-60

3

0

0.5

1

1.5

2

2.5

3

3.5

Time (sec) 1

0

0

-1 -1.5

3

0

-40

0.5

-0.5

2.5

-20

Time (sec)

0.5

2

Time (sec) 20

Time (sec)

Aileron (deg)

2.5

-10

-60

1

2

0

-0.5

Aileron (deg)

Roll Rate - p (deg/sec)

0

0.5

-60

10

-20

0

-40

Time (sec)

Time (sec) 20

-60

1.5

Flight Test Data AAA State Space Model AVL State Space Model

-20

-120

0

Roll Rate - p (deg/sec)

0

0

Flight Test Data AAA State Space Model AVL State Space Model

40

Bank A ngle - φ (deg)

Bank Angle - φ (deg)

Bank Angle - φ (deg)

0 -20

-100

Roll Mode Simulation Response - AAA Model vs AVL Model Compare with Flight Test Response. Flight Test Date : Thu Dec 20 2007; Time : 13-29-35

50

Flight Test Data AAA State Space Model AVL State Space Model

20

Flight Test Data Simulation Input

-1 -1.5

-1

Flight Test Data Simulation Input

-2 -3

-2

-2 -2.5

-2.5

0

0.5

1

1.5

2

2.5

3

3.5

Time (sec)

Figure 11.

4

4.5

-4

0

0.5

1

1.5

2

2.5

3

0

0.5

1

1.5

Time (sec)

2

2.5

3

3.5

Time (sec)

AAA vs. AVL vs. Flight Test Roll Dynamics Comparison

The comparison clearly shows that the roll dynamics predicted by the AVL model is more accurate than the AAA model. In fact, if only the time constant characteristics are considered, both models demonstrate similar performance with very fast responses that match closely with the flight test response. This is consistent with the time constant values calculated from the previous state space model analysis. The amplitude of the roll rate response, however, is overestimated by the AAA model and thus results in larger bank angle values than the flight test results. 12 American Institute of Aeronautics and Astronautics 020208

E. Summary of Comparison Results Through the comparison study on different modes analysis, the following summary can be made: ƒ ƒ ƒ ƒ ƒ

Among the three different modeling methods, the AAA model provides the closest predictions to the flight test results of the Dutch roll, short period and Phugoid mode dynamics. The lack of an engine model on the AVL method makes its approximation results less accurate for the short period and Phugoid mode estimates. The missing longitudinal derivatives in the AVL model makes it less useful for further longitudinal simulation activities. Side-by-side comparison techniques are introduced to truly compare the roll dynamics. These techniques provide great benefits in the validation process but these can only be done on the MATLAB/Simulink platform. According to the side-by-side comparison, the roll dynamic given by the AVL model is more precise than the AAA model.

From this study, it is concluded that overall the AAA model does the best job of modeling the aircraft mode dynamics. However, there are still drawbacks to this technique. First, this state space model is a simplified linear model. Second, it assumes these is no coupling between longitudinal and lateral dynamics. Third, it uses an earth fixed coordinate system, which does not consider the rotation of the earth. Finally, the state space model does not provide full state system outputs, which is required by the Flight Gear software if a visualization of simulation response is demanded. Overall, this state space simulation model does not provide real-time simulation capabilities with an external pilot joystick device, which is a capability provided by the Piccolo simulation platform. Due to the limitations of this linear model, it has limited potential for full state simulation activities and the HIL simulation. To develop a good simulation platform, it is essential to provide full access and flexibility for users to conduct troubleshooting on modeling inputs. In addition, a high fidelity simulation model that allows HIL simulation in realtime is also important. Currently, none of the modeling methods described can satisfy all the requirements. A solution for this is to use a 6DOF nonlinear simulation model, which would potentially allow real-time simulation in conjunction with external joystick commands and 3-D visualization on aircraft responses. Based on the comparison results, the team decided to choose the longitudinal derivatives given from the AAA model and the lateral derivatives provided by the AVL model for the 6DOF nonlinear model development. Table 14 shows the list of derivatives used to develop the 6DOF nonlinear model. Note that all the power relevant derivatives, i.e. CTxu, are not shown here as they are not required for the 6DOF nonlinear model. Table 14. Stability Derivatives used for the 6DOF Nonlinear Model Coefficients for Zero Angle of Attack (rad-1)

CL0

CD0

Cm0

0.0000 0.0513 0.0020 Longitudinal Derivatives from AAA Model(rad-1)

CDu

CDα

CLu

CLα

CLα

CLq

Cmu

Cmq

0.0011

0.0863

0.0017

4.5465

1.8918

5.5046

0.0002

-8.0532

Cmα

Cmα

Cmq

CDδ e

CLδ e

Cmδ e

-0.3937

-4.3787

C yβ

Cyp

C yr

Clβ

Cl p

Clr

Cnβ

Cn p

-0.2707

0.0194

0.2531

-0.0314

-0.5858

0.0743

0.1052

-0.0387

Cnr

C yδ a

C yδ r

Clδ a

Clδ r

Cnδ a

Cnδ r

-0.1156

0.00

0.2228

0.3707

0.0219

-0.0088

-0.1003

0.00 0.3792 -8.0532 Lateral Derivatives from AVL Model (rad-1)

-0.8778

13 American Institute of Aeronautics and Astronautics 020208

VII. Development of the 6DOF Nonlinear Model using MATLAB/Simulink In this section, the development of a 6DOF nonlinear model for a 1/3 scale Yak-54 aircraft is presented. This 6DOF nonlinear model is constructed using the derivatives shown in Table 14. This 6DOF model is developed based on the MATLAB/Simulink platform, as it supports many handy tools that allows for the rapid development of this comprehensive model. A. Coordinate System The “6DOF ECEF (Quaternion)” blockset provided from the MATLAB/Simulink Aerospace Blockset18 is chosen to implement the equations of motion for the 6DOF nonlinear aircraft dynamics model. This model utilizes an Earth-Centered Earth-Fixed (ECEF) coordinate system19 which can be illustrated in Figure 12. The origin of the ECEF coordinate is considered at the center of the earth. This ECEF frame is rotating about the Earth-Centered Inertial (ECI) reference frame. The representation of the ECEF frame from the ECI frame is simplified by only considering a constant rotation about the Z axis. The body of interest is assumed to be rigid and is represented in the ECEF frame.

Figure 12. ECEF Coordinate System18

B. Servo Dynamics Module The servo dynamics module simulates the dynamics characteristics of the actual servo movements from a given command. This refers to the response delay time between the command signals and the actual response time of the servos. To quantify the time delay of the servo, a servo delay tester was designed for this purpose. The servo delay tester is a device to measure the amount of time from when the command is sent through the joystick to the time when the servo actually reaches to the commanded position. The tester consists of a 555 analog circuit timer that is triggered by the joystick command. The tested servo is mounted on the device and is attached with a small paddle that moves in the same direction the servo is to travel. This paddle makes contact with two micro switches which can be adjusted to set at the desired travel angles required for the test. When a joystick command is given, it triggers the timer circuit and starts the timing cycle until the micro switch on the end side is actuated. The timer signal is displayed on an oscilloscope so the command triggering time and the servo traveling time can be compared and quantified. On the Yak-54 platform, Hitec HS-5985MG servos were used to actuate all control surfaces. They are coreless digital servos driven by metal gears for high torque delivery. Because this type of servo is designed and built for fast response, it is reasonable to assume that the dynamics of the system behaves as a first order system. For a first order system, the time constant12 is used to describe the system response characteristics. The micro switches were set 10 degrees apart so that the total elapsed time for the 10 degree travel required for the Hitec servo was measured. According to the test results, the average traveling time for a 10 degrees movement is 200 ms. Note that the result measured in this test is the total elapsed time to reach the final state, which is different than the definition of the time constant. According to Ref. 12, for a first order system, the response time to reach the final value is about four times the time constant. Thus the time constant for this servo can be obtained as: 200/4 = 5 ms. Using this result, a first order transfer function, as shown in Eq. (9), can be developed to simulate the dynamics of the servo response.

T (s) =

1 1 20 = = TS + 1 0.05S + 1 S + 20

(9)

C. The 6DOF Equations of Motion The core of the 6DOF nonlinear model are the equations of motion (EOM) used in the simulation model. The EOM system is implemented using the standard six degree of freedom nonlinear differential equations for a conventional aircraft. Details discussion of the 6DOF equations can be found in Ref. 5, 20, 21 and 22. The 6DOF nonlinear differential equations are:

14 American Institute of Aeronautics and Astronautics 020208

m ( u − vr + wq ) = mg x + FAx + FTx Force Equations:

m ( v + ur − wp ) = mg y + FAy + FTy

(10)

m ( w − uq + vp ) = mg z + FAz + FTz I xx p − I xz r − I xz pq + ( I zz − I yy ) rq = LA + LT Moment Equations:

I yy q + ( I xx − I zz ) pr + I xz ( p 2 − r 2 ) = M A + M T

(11)

I zz r − I xz p + I xz qr + ( I yy − I xx ) pq = N A + NT Quaternion methods19, 20, 23, 24 are used on this 6DOF model so that the singularity problem at a pitch angle of 90 degree can be avoided. In fact, it also provides the best performance for simulations in terms of computational efficiency,23, 25, 26 which is a great advantage when real-time simulation is required. The transformation matrix and the kinematics equations using quaternions are expressed in equations (12) and (13). Transformation Matrix using quaternions:

⎡ 2 2 2 2 ⎡ xb ⎤ ⎢e0 + ex − ey − ez ⎢y ⎥ = ⎢ 2 e e −e e ⎢ b⎥ ⎢ ( x y z 0) ⎢⎣ zb ⎥⎦ ⎢ 2 e e − e e ⎣ ( x z y 0)

2 ( ex ey + ez e0 ) e02 − ex2 + ey2 − ez2 2 ( ey ez − ex e0 )

2 ( ex ez + ey e0 ) ⎤ ⎡ x ⎤ ⎥ f ⎢ ⎥ 2 ( ey ez + ex e0 ) ⎥ ⎢ y f ⎥ ⎥ ⎢ ⎥ e02 − ex2 − e y2 + ez2 ⎦⎥ ⎣ z f ⎦

(12)

Note the subscript “b” refers to body-fixed frames, and the subscript “f” refers to the Earth-fixed frame. Kinematics equations using quaternions: ⎡ e0 ⎤ ⎡ 0 − p − q −r ⎤ ⎡ e0 ⎤ ⎢ e ⎥ ⎢ r −q ⎥⎥ ⎢⎢ ex ⎥⎥ ⎢ x⎥ = 1 ⎢p 0 ⎢e y ⎥ 2 ⎢ q −r 0 p ⎥ ⎢ey ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥  e − 0 r q p ⎣ ⎦ ⎣ ez ⎦ ⎣ z⎦

(13)

The Euler angle representation can be calculated from the quaternions using the following expression:

⎡atan2 ⎡ 2 ( e0 ex + ey ez ) , ( e02 + ez2 − ex2 − e2y ) ⎤ ⎤ ⎣ ⎦⎥ ⎡φ ⎤ ⎢ ⎢θ ⎥ = ⎢ ⎥ asin ⎡⎣ 2 ( e0 ey − ex ez ) ⎤⎦ ⎢ ⎥ ⎢ ⎥ ⎣⎢ψ ⎦⎥ ⎢atan2 ⎡ 2 ( e e + e e ) , ( e 2 + e 2 − e 2 − e2 ) ⎤ ⎥ x y z ⎦⎥ 0 ⎢⎣ ⎣ 0 z x y ⎦

(14)

The function atan2 in Eq. (14) is a two argument arctangent function that determines the result in the proper quadrant by examining the signs of both arguments. From a mathematical viewpoint, Eq. (10) and (11) form a set of six differential equations with six unknown variables (u, v, w, p, q, r). Each variable is presented in different equations and interacts with the others. Therefore, the equations cannot be solved individually. The total solution of the system can be obtained only by applying the numerical integration to each equation for each given time step. For the forces equations, the inputs for the system are the three major forces applied on the aircraft. Those are the aerodynamic forces, thrust forces and gravity forces. These forces are nonlinear and time variant. The components of each force will be broken down in detail and presented shortly. The inputs for the moment equation are the three moment terms that apply to the three aircraft body axes. These moments are contributed by the forces, as described above, with respect to the center of gravity (C.G.) of the 15 American Institute of Aeronautics and Astronautics 020208

aircraft. As the gravity forces are defined at the C.G. point, no moment is generated by the gravity forces. In addition, when modeling the engine thrust forces for the Yak-54, assumption is made that the thrust line passes through the C.G. point. Therefore, the moment terms that contribute to the moment equations are only associated with the aerodynamic forces. These aerodynamic moment terms will be broken down one-by-one and presented later in this section. When the variables p, q, r are available from the 6DOF equations, they will then be applied to the quaternion kinematic equations (13) to update the quaternion terms. The updated quaternion terms are then used for vector coordinate transformations through Eq. (12). The latest orientations of the aircraft are then provided from Eq. (14). D. Aerodynamics Module As mentioned, there are three major forces that apply to the aircraft system. Aerodynamic forces are one of the major forces, and these forces create aerodynamic moments that contribute to the moment equations. The main purpose of this aerodynamics module is to provide the values of the aerodynamic forces and moments at the X, Y, Z axis in the aircraft body frame. The aerodynamics forces compose of three forces which are the lift, drag and sideforces. These forces generate moments with respect to the center of gravity along the X, Y, Z axis and are described as the rolling moment, pitching moment and yawing moment. The component build-up method5, 20 is used here to generate the forces and moments . The total forces and moments which acting on the aircraft are simply the summation of the forces and moments contributed by each component. This method has been widely used and found to be acceptable. The aerodynamics forces and moments are first described as dimensionless coefficients, which are associated with the stability and control derivatives given from the AAA and AVL results. The equations used to implement the aerodynamics forces and moments coefficients are presented one-by-one as follows:

(

)

u c + CLδ e δ e + CLq q + CLα α × V∞ 2V∞

Lift Coefficient:

CL = CL0 + CLα α + CLu ×

Drag Coefficient:

CD = CD0 +

Sideforce Coefficient:

CY = CYβ β + CYδ a δ a + CYδ r δ r + CYp p + CYr r ×

Pitching Moment Coefficient:

Cm = Cm0 + Cmα α + Cmu ×

Rolling Moment Coefficient:

b 2V∞ b Cn = Cnβ β + Cnδ a δ a + Cnδ r δ r + Cn p p + Cnr r × 2V∞

Yawing Moment Coefficient:

CL2 u + CDα α + CDu × + CDδ e δ e π eAR V∞

(

(15) (16)

)

b 2V∞

(

(17)

)

u c + Cmδ e δ e + Cmq q + Cmα α × V∞ 2V∞

(18)

(

)

(19)

(

)

(20)

Cl = Clβ β + Clδ a δ a + Clδ r δ r + Clq p + Clr r ×

where e is the Oswald coefficient and AR is the main wing aspect ratio The final values of the forces and moments are then computed from the dimensionless coefficients using the expressions below:

L = CL qS D = CD qS FAy = CY qS

LA = Cl qSb M A = Cm qSc N A = Cn qSb

(21) (22) (23)

(24) (25) (26)

where q = 1 ρ∞V∞2 is the dynamic pressure

2

The forces components shown in Eq. (21), (22) and (23) are defined in the stability axis. To make it useable for Eq. (10) which is defined at body axis, a vector transformation is needed to transfer the stability axis forces to body axis forces using the equations (27).

16 American Institute of Aeronautics and Astronautics 020208

⎡ FA ⎤ ⎡C os (α ) 0 − Sin (α ) ⎤ ⎡ − D ⎤ ⎢ x⎥ ⎢ ⎥ ⎢ ⎥ 1 0 ⎢ FAy ⎥ = ⎢ 0 ⎥ × ⎢ FAy ⎥ ⎢ ⎥ ⎢ Sin α ( ) 0 C os (α ) ⎦⎥ ⎢⎣ − L ⎥⎦ ⎢⎣ FAz ⎥⎦ ⎣

(27)

E. Engine Dynamics and Thrust Force Module Another major force that applies to the aircraft system is the thrust force. In this section, the modeling techniques of the engine dynamics and thrust forces are discussed. The objective of this module is to calculate the thrust forces for any given throttle position input. The Yak-54 is driven by a 26”x10” propeller with a gas powered, single 3W 80cc XI CS engine. Since no power measurement device is available to measure the engine power, a simplified method is used that assumes the engine RPM is only a function of the throttle position in any conditions. An experiment was conducted on ground to measure the thrust and RPM values with respect to different throttle position; the results are shown in Figure 13. Keep in mind that the thrust value shown here is the static thrust. In order to simulate the thrust forces in flight conditions, it is necessary to use dynamic thrust rather than static thrust. The variation of dynamic thrust depends not only on the rotational speed of propeller (RPM) but also on the forward flight speed (relative airspeed). For this reason, the advance ratio6, 21 term J is introduced. The advance ratio is a dimensionless term defined as a function of the forward airspeed (Vx) and the propeller rotational speed per second (n) as expressed below. The notation d refers to the diameter of the propeller.

J=

V n×d

(28)

Using the blade element theory6, 21, the thrust coefficient (CT) and power coefficient (CP) are defined and used to calculate the thrust forces and absorbed power of the propeller as a function of advance ratio. The thrust and power coefficients curves mainly depends on the four design parameters of the propeller, which are its size, the number of blades, the airfoil shape and the pitch angle values. For industry standard products, these coefficient data are available from the propeller manufacture. Wind tunnel testing is used to measure and calculate the CT and CP coefficients experimentally. CFD methods could also be used to simulate and estimate those coefficients for the propeller. However, either the wind tunnel or the CFD methods are both costly and time intensive. As the propeller used on the Yak-54 is a hobby grade product, no engineering data is available. Using a wind tunnel or a CFD method will not meet the required schedule and the available budget and these methods are thus beyond the scope of this project. For these reasons, a computer program, JavaProp27, is used to estimate these coefficients data based on the configuration of the four propeller design parameters. The results are illustrated in Figure 13. Propeller Coefficients and Efficiency

RPM Thrust (lbs)

7000

45.00

90.00

0.0700

80.00

25.00

3000

20.00

Thrust & Pow er Coefficient

30.00

4000

70.00

0.0600

35.00 5000

Thrust (lbs)

Engine RPM

0.0800

40.00

6000

60.00

0.0500

50.00 0.0400 40.00 0.0300

30.00

15.00

2000

0.0200

20.00

10.00 1000 0 0

20

40 60 Throttle Position (%)

80

5.00

0.0100

0.00

0.0000 0.0000

100

Cp [-] Ct [-] Effec [%]

Propeller Effeciency (%)

Engine RPM and Thrust Curve

10.00

0.2000

0.4000

0.6000

0.8000

1.0000

1.2000

0.00 1.4000

Ad cance Ratio (J)

Figure 13. Engine Performance Data The propeller coefficient curve shows a typical performance chart for a fixed pitch propeller. As the advance ration increases, the thrust coefficient decreases and thus less thrust force is generated. The generated thrust will become zero when it reaches to a combination of airspeed and RPM value. Beyond this point, the sign of the thrust 17 American Institute of Aeronautics and Astronautics 020208

changes to the opposite direction and creates reverse thrust. This is known as the windmill effect. The performance chart also reveals that this propeller provides the best efficiency when operated at the 0.85 advance ratio value. The logic of the engine module is described as follows. For a given throttle input, the corresponding RPM output are determined according to the RPM curve. This RPM value is then used to calculate the advance ratio with the current airspeed feedback from the simulation model. The thrust coefficient can then be obtained from the thrust coefficient curve according to the calculated advance ratio value. Finally, the thrust force is computed using Eq. (29). The overall modeling method of the engine module is illustrated in Figure 14. The notation ρ∞ used below refers to the air density in current altitude and is provided from the atmosphere module, which will be discussed shortly.

T = ρ ∞ n 2 d 4CT

(29)

V∞ δT

LUT

δT vs RPM

J=V/n*d

ρ∞

LUT

J vs CT T = ρ∞n2d4CT

Thrust (lbs)

Figure 14. Block Diagram of the Engine and the Thrust Model F. Gravity Module Gravity is the last major force that is applied to the aircraft system. In most simple simulation models, gravity is always assumed to be constant regardless of the aircraft position and altitude with respect to the Earth. In this 6DOF nonlinear model, the World Geodetic System28 (WGS84) is used to calculate the Earth’s gravity at a specific location given from the inputs. The concept of the World Geodetic System was first started during the late 1950s through the effort of United States Department of Defense, aided by scientists from other institutions and countries. It was motivated by the need for a unified world geodetic system that could be used in geodesy and navigation worldwide. It was essential for global maps navigation, aviation and geography. The WGS model was first developed in 1960 and was called WGS60. Since then, major updates have been made to improve the fidelity and the accuracy of the model using the latest technology and mathematic model. The WGS72 was updated in 1972. The latest model is the WGS84,28 which was published in the early 1980s. The WGS84 model is currently being used as the reference system for the Global Positioning System (GPS). The “WGS84” model18 is available from the Aerospace Blockset in MATLAB/Simulink. The model takes the input of the position in geodetic format (latitude, longitude and altitude) and provides a three vectors output for the gravity in the North-East-Down (NED) coordinate system.19 Transformation is required from the NED frame to the body frame before applying the gravity values to Eq. (10) for the forces equations. This transformation is done through Eq. (12) using the quaternion method. The construction of the WGS84 in MATLAB/Simulink is shown in Figure 15. G. Atmosphere Module During the simulation process, updated atmospheric data are frequently required to give feedback to some of the modules to provide the necessary information for computation. For this reason, an atmosphere module is required to provide the latest atmospheric data for the current altitude. The 1976 COESA29 (Committee on Extension to the Standard Atmosphere) atmospheric model is used here to fulfill this requirement. The work of the U.S. COESA was established in 1953, and major revisions were made in 1958, 1962, 1966 and 1976 accordingly. Many U.S. government organizations contributed to this work, including NOAA (National Oceanic and Atmospheric Administration), NASA (National Aeronautics and Space Administration), U.S. Air Force, industry, industries, research institutions, and universities. The COESA atmosphere model is identical to the Standard Atmosphere of the International Civil Aviation Organization (ICAO) up to 32 Km, and the ISO (International Organization for Standardization) standard up to 50 Km. The COESA 1976 is an idealized, steady-state representation of the earth’s atmosphere from the surface to 1000 km for moderate solar activity. The “COESA Atmosphere” model18 is available in MATLAB/Simulink Aerospace 18 American Institute of Aeronautics and Astronautics 020208

Blockset. This MATLAB/Simulink COESA model implements the mathematical representations of the 1976 COESA values and provides absolute temperature, pressure, density and speed of sound for a given geodetic altitude. The use of this COESA atmosphere model in MATLAB/Simulink is shown below.

Figure 15. COESA Atmosphere and WGS84 Gravity Model in MATLAB/Simulink18

VIII. Validation of 6DOF Nonlinear Simulation Model In this section, the flight test data are used to validate the 6DOF nonlinear model simulation responses. This study also includes the simulation results from the state space linear model for the purpose of comparing the accuracy between the nonlinear model and the linear model. The side-by-side comparison technique, as introduced in Section VI, is used here throughout the validation process. A. Validation of the Roll Mode Response Body Rate Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

Roll Mode Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

10

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

40

Bank Angle - φ (deg)

Roll Rate - p (deg/sec)

50

30 20 10 0 -10

0 -10

-30 -40 -50

0

0.5

1

1.5

2

2.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-20

3

0

0.5

1

Pitch Rate - q (deg/sec)

Roll Rate - p (deg/sec)

0 -10 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-20 -30 -40

0.5

1

1.5

2

5

2.5

0

3

0

0.5

1

Yaw Rate - r (deg/sec)

Aileron (deg)

0 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-1 -1.5 -2

2

2.5

3

5

0

-5 0

0.5

1

1.5

2

2.5

Flight Test Data P6DOF Nonlinear Model State Space Linear Model

10

3

0

0.5

1

Euler Angles Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

3

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0.6

X-Acceleration - G

30 20 10 0

0.5 0.4 0.3 0.2 0.1

0

0.5

1

1.5

2

2.5

0

3

0

0.5

1

Time (sec)

1.5

2

2.5

3

Time (sec)

12

0.2

10

Y-Acceleration - G

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

11

9 8 7

0

0.5

1

1.5

2

2.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0.1

0

-0.1

-0.2

3

0

0.5

1

Time (sec)

1.5

2

2.5

3

Time (sec)

360

-0.4 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

270

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-0.6

Z-Accelerati on - G

315

225 180 135 90

-0.8 -1 -1.2 -1.4 -1.6

45 0

2.5

0.7 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

40

6

2

Acceleration Simulation Responses - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

50

-10

1.5

Time (sec)

Time (sec)

Bank Angle - φ (deg)

1.5

Time (sec) 15

-0.5

Pitch Angle - θ (deg)

3

10

-5 0

0.5

Heading Angle -ψ (deg)

2.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

15

Time (sec)

-2.5

2

20

10

-50

1.5

Time (sec)

Time (sec)

0

0.5

1

1.5

2

2.5

-1.8

3

0

0.5

1

Time (sec)

1.5

Time (sec)

Figure 16. Validation of the Roll Mode Responses

19 American Institute of Aeronautics and Astronautics 020208

2

2.5

3

Figure 16 shows a set of comparisons from the roll dynamics test. Both the 6DOF nonlinear and the linear model demonstrate almost identical performance in the roll dynamics response. Examination of the pitch rate response, however, indicates that the 6DOF model matches closer to the flight test data then does the response of the linear model. In fact, the pitch rate response from the linear model has an opposite sign. As a result, the pitch angle drifted in the opposite direction with a large error, while the pitch angle response from the 6DOF model provides trend similar to the flight test data. This demonstrates that the 6DOF model does a better job of simulating the coupling effect between the lateral and longitudinal dynamics. According to the Newton’s second law, a comparison of the acceleration data provides a great insight into the process for validating the forces components estimated by the simulation model. Recall that the force components are the major inputs into the 6DOF EOM system. This comparison provides important information which is very helpful when improving modeling accuracy. For the Y-axis acceleration data, the 6DOF model matches closely with the flight test data, while large amplitude oscillations were seen in the linear model’s response. This indicates that the sideforce estimated by the nonlinear model is fairly accurate. The Z-axis acceleration data given from the 6DOF model also tracks the flight test data very well. Recall that the Z-axis forces in the body frame were computed from Eq. (27) which is a function of the angle of attack (AOA). This AOA value is fed back from the simulation outputs throughout the simulation process. However, an initial condition is required to start the computational process. Since the AOA was not available from the flight test data, an assumption was made that the initial AOA for simulation is the same as the trimmed AOA given from the AAA model. Keep in mind that this trim value is only valid when the aircraft is in a perfect steady level flight condition at 70 knots. Once the simulation process starts, this AOA is being updated and corrected. Therefore the accuracy of the Z-axis acceleration data is being improved through the simulation process, as shown on the acceleration plot in Figure 16. The thrust force has a major contribution to the X-axis force. From the X-axis acceleration data, it is revealed that the X-axis force is overestimated in the 6DOF model. This is mostly due to the simplified engine model that produces extra thrust than normal. B. Validation of the Dutch Roll Mode Response Body Rate Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

Dutch Roll Mode Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response

20

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

350

Roll Rate - p (deg/sec)

Heading Angle -ψ (deg)

360

340

330

320

0 -10 -20 -30

0

0.5

1

1.5

2

2.5

3

3.5

4

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

10

4.5

0

0.5

1

1.5

2

Pitch Rate - q (deg/sec)

Yaw Rate - r (deg/sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

50

0

-50

4

4.5

0.5

1

1.5

2

2.5

3

3.5

4

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0 -5 -10 -15

0

4.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec) 100

10 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

5

Yaw Rate - r (deg/sec)

Rudder (deg)

3.5

5

Time (sec)

0

-5

-10

3

10

100

-100

2.5

Time (sec)

Time (sec)

0

0.5

1

1.5

2

2.5

3

3.5

4

0

-50

-100

4.5

Flight Test Data P6DOF Nonlinear Model State Space Linear Model

50

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Time (sec)

Figure 17. Validation Results of Dutch Roll Mode Test Responses Figure 17 shows the dynamics responses from a Dutch roll mode test. The result indicates that the damping ratio is underestimated by both the nonlinear and the linear model, which is consistent with the conclusion from the eigenvalues analysis discussed in Table 13. From the roll rate response, the coupling between yaw and roll can be examined clearly. It is also found that both simulation models overestimate the coupling effect on the roll dynamics due to yaw. A derivative tuning technique is now introduced. As it is known that some of the derivatives have a significant affect on a specific dynamics response, it is possible to tune the values of those derivatives to adjust the dynamic response till performance goals are satisfied. The damping ratio of the Dutch roll mode is most affected by the yawdamping derivative Cnr. The current used value of the Cnr derivatives is -0.1156. To increase the damping, a bigger magnitude for Cnr is desired so that a higher yawing moment is generated in the opposite direction to quickly correct the oscillation in the yaw rate. The power control derivative Cnδr also has some effect on the yaw rate oscillations when the rudder is activated. In fact, different aspects of the yaw rate response is affected by the Cnr and Cnδr 20 American Institute of Aeronautics and Astronautics 020208

derivatives, and thus an iterative process is required to tune each of the derivatives until satisfactory performance is achieved. During a Dutch roll mode test, a sideslip angle is introduced along the oscillation in the yaw direction. The roll response is generated due to this sideslip angle effect. Therefore, the roll and yaw coupling dynamics is mainly determined by the dihedral derivative, Clβ. After several iterative tuning cycles, a satisfactory result was achieved from the simulation results that gives the best match with the flight test data using values of Cnr increased by 150%, Cnδr increased by 40% and Clβ decreased by 30%. The change is reflected in both the nonlinear and linear models and the new simulation responses are shown in Figure 18. From the new simulation response, it clearly shows that the damping and the coupling responses have been improved significantly. In the pitch rate data, two abrupt responses were seen in the flight test. It is confirmed that no elevator input was commanded during this test so it is most likely caused by the wind. Ideally, the pitch rate response in a Dutch roll mode test should be remained closely to the trim condition. The 6DOF model had successfully demonstrated this performance with a better pitch angle estimation. The 6DOF model also provides better estimations on acceleration data than the liner model. The overestimated thrust force issue is observed again in the X-axis acceleration data. The error of the Y-axis acceleration data given from the linear model is almost 400%. Dutch Roll Mode Response - 6DOF Nonlinear Model vs State Space Linear Model Cnr increase 150% ; Cnδr increase 40% ; Clβ decrease 30%

Body Rate Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Cnr increase 150% ; Cnδr increase 40% ; Clβ decrease 30% 20

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

350

Roll Rate - p (deg/sec)

Heading Angle -ψ (deg)

360

340

330

320

0

0.5

1

1.5

2

2.5

3

3.5

4

10 5 0 -5 -10 -15

4.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

15

0

0.5

1

1.5

2

Time (sec)

Pitch Rate - q (deg/sec)

Yaw Rate - r (deg/sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

20 0 -20 -40 -60 0

0.5

1

1.5

2

2.5

3

3.5

3

3.5

4.5

4

0

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-10 -15

4.5

0

0.5

1

1.5

2

4

4.5

Time (sec) 60

Yaw Rate - r (deg/sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

5

Rudder (deg)

2.5

4

-5

Time (sec)

0

-5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Flight Test Data P6DOF Nonlinear Model State Space Linear Model

40 20 0 -20 -40 -60 -80

0

0.5

1

1.5

2

Time (sec)

Euler Angles Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Cnr increase 150% ; Cnδr increase 40% ; Clβ decrease 30%

X-Acceleration - G

-4 -6

0.5

1

1.5

2

3.5

4

4.5

2.5

3

3.5

4

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0.4

-2

0

3

Acceleration Simulation Responses - 6DOF Nonlinear Model vs State Space Linear Model Cnr increase 150% ; Cnδr increase 40% ; Clβ decrease 30%

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0

-8

2.5

Time (sec)

2

Bank Angle -φ (deg)

3.5

5

10

-10

3

10

40

-80

2.5

Time (sec)

60

4.5

0.3 0.2 0.1 0 -0.1

Time (sec)

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec) 20 2

Y-Acceleration - G

Pitch Angle - θ (deg)

3

15 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

10

5

1 0 -1

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-2 -3 -4

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

-5

Time (sec)

Z-Acceleration - G

Heading Angle -ψ (deg)

340

330

0

0.5

1

1.5

2

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

350

320

0

Time (sec)

360

2.5

3

3.5

4

4.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-0.5

-1

-1.5

-2

0

0.5

1

Time (sec)

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Figure 18. Validation Results from the Dutch Roll Mode Test Incorporating new Derivatives Values C. Validation of Short Period Mode Response Figure 19 shows the comparison for a short period mode test. From the pitch rate data, all three sets of data demonstrated very high damping. The damping from the 6DOF nonlinear model, however, was not as high as seen in the other responses. In addition, the peak of the first oscillation in the response to the elevator input was slightly higher than the others. As a result, an overshoot in pitch angle was introduced. After the elevator input was released, the pitch angle response slowly drifted away from the flight test data due to its slowly damped response. 21 American Institute of Aeronautics and Astronautics 020208

The derivative tuning technique is applied here to improve the short period mode dynamics. In this case, the pitch damping derivatives Cmq and the control power derivatives Cmδe are adjusted, as the pitch damping response is most affected by Cmq. The derivative Cmδe is also required for tuning because the pitch rate response is affected by both the Cmq and Cmδe derivatives, which is similar to the case of Cnr and Cnδr. The new simulation responses are shown in Figure 20 with Cnr increased 100% and Cnδr increased 40%. Using the new derivatives, the pitch rate response given from the 6DOF model is now more closely matched to the flight test data. The damping action now responds as fast as the real dynamics. As the accuracy of the pitch rate response improves, so does the pitch angle responses. The error in the final pitch angle values are reduced by utilizing the new derivatives in both simulation models.

10 5 0 -5 0

0.5

1

1.5

2

2.5

3

3.5

4

Pitch Angle - θ (deg)

10 5 0 -5 0

0.5

1

1.5

2

Pitch Rate - q (deg/sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

-20 -40

0

0.5

1

1.5

2

3.5

4

4.5

Elevator (deg)

10 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

5

0

-5

-10

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Figure 19. Validation Responses for the Short Period Mode Body Rate Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Cmq increase 100% ; Cmδe increase 40% Flight Test Data 6DOF Nonlinear Model State Space Linear Model

5 0 -5 -10 -15 -20

4.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

Pitch Rate - q (deg/sec)

Pitch Rate - q (deg/sec)

3

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

20 0 -20 -40

0

0.5

1

1.5

2

2.5

3

3.5

4

20 0 -20 -40 -60

4.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

40

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec) 4

Yaw Rate - r (deg/sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

5

Elevator (deg)

2.5

60

40

0

-5

0

0.5

1

1.5

2

2.5

3

3.5

4

Flight Test Data P6DOF Nonlinear Model State Space Linear Model

2 0 -2 -4 -6

4.5

0

0.5

1

1.5

2

2.5

Time (sec)

3

3.5

4

4.5

Time (sec)

Acceleration Simulation Responses - 6DOF Nonlinear Model vs State Space Linear Model Cmq increase 100% ; Cmδe increase 40%

Euler Angles Simulation Response - 6DOF Nonlinear Model vs State Space Linear Model Cmq increase 100% ; Cmδe increase 40%

0.6

2 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0 -2 -4

X-Acceleration - G

4

Bank Angle - φ (deg)

4.5

0

Time (sec)

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0.4

0.2

0

-6 -8

-0.2

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Time (sec) 0.1

25 Flight Test Data 6DOF Nonlinear Model State Space Linear Model

20 15

Y-Acceleration - G

Pitch Angle - θ (deg)

4

Time (sec)

10

10 5 0

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

0.05 0 -0.05 -0.1

-5 -10

-0.15

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

0.5

1

1.5

2

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

15

10

5

0.5

1

1.5

2

3.5

4

4.5

2.5

3

3.5

4

4.5

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

1 0 -1 -2 -3

0

3

2

Z-Acceleration - G

20

0

2.5

Time (sec)

Time (sec)

Heading Angle -ψ (deg)

3.5

20

-60

60

-10

3

40

Time (sec)

-60

2.5

Time (sec)

Roll Rate - p (deg/sec)

Pitch Angle - θ (deg)

15

15

60

10

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

20

Flight Test Data 6DOF Nonlinear Model State Space Linear Model

20

-10

Short Period Mode Response - 6DOF Nonlinear Model vs State Space Linear Model Cmq increase 100% ; Cmδe increase 40%

25

-10

Short Period Mode Response - 6DOF Nonlinear Model vs State Space Linear Model Compare with Flight Test Response. 25

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Time (sec)

Time (sec)

Figure 20. Validation of the new Derivatives Values using Short Period Test Responses A yaw rate oscillation in the flight test data was seen that interacted with the roll dynamics. It was confirmed that no aileron and rudder inputs were commanded during this test. Therefore, it is reasonable to infer that the lateral responses were most likely caused by the wind, which won’t be shown in the simulation responses as the current simulation model does not simulate any disturbances. Basically, if one assumes a no wind condition, the 22 American Institute of Aeronautics and Astronautics 020208

aircraft should remain in a level and straight flight condition, which is consistent with the current responses given by both simulation models. In the comparison of acceleration data, the Z-axis response demonstrates that the 6DOF simulation response is almost identical with the flight test results. This points out that the estimation of lift forces is very accurate. The tendency of the X-axis response is similar to that seen in flight test, with some offsets seen at the beginning. As discussed before, the error could be introduced by the thrust model and the uncertain values of AOA. The bottom line clearly shows that all the acceleration data calculated from the state space liner model are not validated. D. Validation of the Phugoid Mode Response An incomplete Phugoid mode flight test was conducted. Although a complete cycle of the Phugoid mode was not available, the flight test data is used in an attempt to validate the airspeed and altitude simulation responses over a long simulation time. Figure 21 shows the Phugoid mode comparison results and the corresponding acceleration responses. Accelerations Simulation Response - 6DOF Nonlinear Model Simulation Response Compare with Flight Test Response

Phugoid Mode Simulation Response - 6DOF Nonlinear Model Simulation Response Compare with Flight Test Response

0.6

Flight Test Data 6DOF Nonlinear Model

80 70 60 50

0.4 0.3 0.2 0.1

40 30

Flight Test Data 6DOF Model Response

0.5

X-Acceleration - G

True Air Speed - (knots)

90

0 -0.1 0

2

4

6

8

10

12

14

0

2

4

6

Y-Acceleration - G

ASL Altitude (ft)

Flight Test Data 6DOF Model Response

1700 1600 1500 1400

12

14

Flight Test Data 6DOF Model Response

0.1 0.05 0 -0.05 -0.1 -0.15 -0.2

0

2

4

6

8

10

12

14

0

2

4

6

10

12

14

1

8

Z-Acceleration - G

Flight Test Data 6DOF Model Response

6 4 2 0 -2 -4 -6 0

8

Time (sec)

Time (sec)

Elevator (deg)

10

0.15

1800

1300

8

Time (sec)

Time (sec)

2

4

6

8

10

12

-2 -3 -4

14

Flight Test Data 6DOF Model Response

0 -1

0

2

4

6

8

10

12

14

Time (sec)

Time (sec)

Figure 21. Validation Results from the Phugoid Mode Tests Overall, the simulation responses tend to follow the trend of the flight test results with larger errors when compared to the lateral-directional response. The main reasons for this is that the current 6DOF model does not provide sufficient accuracy in predicting the airspeed and altitude change as explained as follows: ƒ

ƒ

The variation in airspeed during the Phugoid mode test about the trim speed (70 knots) is too high. Recall that the 6DOF nonlinear model was built using the derivatives calculated at the 70 knots trim speed condition. The simulation model is no longer an accurate way to predict the aircraft dynamics when the flight condition is beyond the trim condition. In this Phugoid mode test, the airspeed varied from 70 knots to 40 knots. The variation in airspeed and altitude has a coupling effect, and these dynamics are mainly controlled by the elevator and throttle inputs. Currently, a simplified modeling technique is applied to the engine model in response to the throttle input. Because the engine dynamics have highly nonlinear characteristics, it is not trivial to develop an engine model with high fidelity. As the inaccurate predication of airspeed will couple with the altitude response and vice versa, this error will build with time, making the simulation response less accurate.

E. Summary of Validation Results Examining the validation results for each mode’s dynamics, the following summary can be made: ƒ ƒ ƒ

The 6DOF nonlinear model shows very promising performance in duplicating the roll mode, Dutch roll mode and short period mode dynamics. The 6DOF model provides more accurate simulation responses than the linear mode in every aspect. The 6DOF model also gives a very accurate predication of the coupling effect between the longitudinal and lateral dynamics, which the state space linear model does not consider. 23 American Institute of Aeronautics and Astronautics 020208

ƒ

ƒ ƒ ƒ

The derivative tuning technique is introduced and its capability to improve the accuracy of the simulation model successfully demonstrated by adjusting the values of the relevant derivatives. This points out one of the great advantages of using this 6DOF nonlinear model, which the SCCS cannot achieve. The lift force and sideforce component estimates are shown to be very good in every dynamic mode comparison. It is observed that the X-axis force is overestimated due to the overestimated thrust forces. This is mostly caused by the simplified engine model. Due to the lack of available stability derivatives at different trim speeds, and the low fidelity of the engine model, the current 6DOF nonlinear model does not provide sufficient accuracy to predict the variation in the airspeed and altitude responses.

Five derivatives are adjusted during the Dutch roll mode and short period mode validation process. The changes to derivatives are summarized in Table 15. These new derivatives are applied to update the state space linear model to recalculate the mode dynamics results using the eigenvalues analysis. The new mode analysis result is then compared with the original results as is shown in Table 15 and 16. Table 15. Summary of Derivatives Change Derivatives WAS NOW % Increment Cnr -0.1156 -0.2890 150 -0.1003 -0.1404 40 Cnδr -0.0314 -0.0220 -30 Clβ Cmq -8.0532 -16.1064 100 -0.8778 -1.2289 40 Cmδe Table 16. Comparison of Mode Analysis with New Derivatives Model Lateral Dynamics Longitudinal Dynamics Spiral Roll Dutch Roll Short Period Phugoid Model ω ω ωn τs τr n n ζ ζ ζ (rad/sec) (rad/sec) (rad/sec) (sec) (sec) Original Model 0.90 9.42 0.42 0.27 55.87 0.04 0.16 7.13 Results New Model 0.36 7.04 0.98 11.6 0.53 0.22 202.83 0.04 Results ------0.27 5.70 1.00 16.62 ------Flight Test The new damping ratio for the Dutch roll mode is now identified as 0.36, which is 125% higher than before, and is 33% higher than determined from the flight test data. From Figure 18, it is already proven that the new simulation response is closer to the true dynamics. This points out that the flight test data reduction method involves some level of error due to the use of graphical techniques. In fact, research using system identification techniques30 is conducted using the same set of flight test data as present in this paper. According to the system identification analysis results31, the damping ratio for the Dutch roll mode is determined as 0.30 with 77% accuracy, which concurs with the final result present here. The newly determined damping ratio for the short period mode is 0.98, which is 7.2% higher than before. This is very reasonable as the actual dynamics from flight test demonstrated a very highly damped response, which could be considered as a critically damped response. For a critically damped system, its damping ratio is equal to 1.0.

IX. Pilot Training Platform It is shown in the previous study that the 6DOF nonlinear model provides good performance. To extend its advantages and applications, a pilot training platform is developed using this MATLAB/Simulink based simulation model. The pilot training platform is constructed by connecting the Futaba pilot console with the MATLAB/Simulink model so that direct simulation control inputs can be given by pilots through a pilot console. The outputs of the simulation are then visualized using the Flight Gear software10. Figure 22 illustrates the pilot training platform. 24 American Institute of Aeronautics and Astronautics 020208

Futaba Pilot Console

MATLAB/Simulink

Flight Gear

6DOF Model Figure 22. Pilot Training Platform

A. Interface Setup The Futaba T9CAP console serves as the input unit, which is identical to the control device used by the pilot in flight testing. A specific USB interface cable for the Futaba pilot console from MileHighWings32 is used. It provides a physical connection between the console and the desktop computer where the 6DOF model is executed. The Flight Gear platform, which has been widely used in various simulation environments2, 33, 34, is chosen here for visualization purposes. This was mainly due to the available tools supported by MATLAB/Simulink that allows data from MATLAB/Simulink to be transmitted to Flight Gear. Options are provided that the data transmission can either be accomplished on the same computer machine or between two different computer platforms through a network connection. After appropriately assigning the required transmitted data, Flight Gear receives the simulation outputs from the 6DOF model at a fixed update rate and constantly visualizes the response within a graphical based aircraft model in real-time. Figure 23 illustrates a screen shot taken while the pilot training platform was running on MATLAB/Simulink. Using this simulation platform, flight training lessons can be given to the pilot for a specific aircraft platform. Other aircraft dynamics models can be implemented if desired by easily replacing those aerodynamics derivatives and relevant geometric data. The Meridian UAV dynamic model will be implemented in this platform for pilot training purpose in the near future.

Figure 23. Screen Shot of the MATLAB/Simulink Pilot Training Platform

X. Summary and Conclusions Three different aircraft dynamic modeling techniques were developed for a 1/3 scale Yak-54. Open loop dynamics flight tests were conducted and results were used to evaluate the accuracy of each modeling method. The best modeling data set was then selected and used to develop the 6DOF nonlinear model using MATLAB/Simulink. A through validation process was conducted to compare the simulation performance between the nonlinear and linear models with the flight test data. The study reveals that the 6DOF model provides better accuracy than the linear model. A pilot training platform has been set up using this MATLAB/Simulink based 6DOF model. This platform allows for pilot training using a definable aircraft dynamic model.

Acknowledgments This research work is supported by the Center of Remote Sensing of Ice Sheets (CReSIS) funded by NSF Grant number AST-042-4589. 25 American Institute of Aeronautics and Astronautics 020208

References 1

Ly, L. and Higashino, S., “Development of a UAV Flight-Test Vehicle at the University of Washington,” 2nd AIAA “Unmanned Unlimited” Systems, Technologies, and Operations – Aerospace Conference, San Diego, California, 15-18 September 2003. 2 Jung, D. and Tsiotras, P., “Modeling and Hardware-in-the-Loop Simulation for a Small Unmanned Aerial Vehicle,” AIAA Infotech Aerospace Conference, Rohnert Park, California, 7-10 May, 2007. 3 Sadraey, M. and Colgren, R., “UAV Flight Simulation: Credibility of Linear Coupled vs. Nonlinear Coupled Equations of Motion,” AIAA modeling and Simulation Technologies Conference, San Francisco, California, 15-18 August 2005. 4 Jodeh, N. M., Blue, Paul A. and Waldron, A. A., “Development of a Small Unmanned Aerial Vehicle Research Platform: Modeling and Simulation with Flight Test Validation,” AIAA Modeling and Simulation Technologies Conference, Keystone, Colorado, 21-24 August 2006. 5 Roskam, Jan, “Airplane Flight Dynamics And Automatic Flight Controls Part I,” DAR Corporation, Lawrence, Kansas, 2003 6 Ward, Donald T. and Strganac, Thomas W., “Introduction to Flight Test Engineering,” 2nd ed, Kendall/Hunt, Dubuque, Iowa, 2001. 7 Anderson, John D., “Fundamentals of Aerodynamics,” 3rd ed., McGraw Hill, 2001 8 Norton , Robert L., “Machine Design: An Integrated Approach,” Prentice Hall, 3rd ed., 2005. 9 Vaglienti, B., Niculescu, M. and Becker, J., “HIL/SIL Simulator for the Piccolo Avionics,” Cloud Cap Technology, Hood River, OR, 5 April 2007. 10 “Flight Gear Flight Simulator,” http://www.flightgear.org/. 11 Becker, Jon, ”Piccolo Dev Interface”, Cloud Cap Technology, Hood River, OR, 2 March 2007 12 Katsuhiko, Ogata, “Model Control Engineering,” 4th ed., Prentice Hall, New Jersey, 2002. 13 “Advanced and Aircraft Analysis”, Version 3.1, DAR Corporation, Lawrence, Kansas, 14 “DAR Corporation”, http://www.darcorp.com/ 15 Willams, J. E. and Vukelich, S. R., “The USAF Stability and Control DATCOM,” Tech. Rep., McDonnell Douglas Astronautics Company, St Louis, CO, October 1979, AFFDL-TR-79-3032. 16 “Athena Vortex Lattice,” http://web.mit.edu/drela/Public/web/avl/. 17 Becker, J., “Creating Vortex Lattice Models for the Piccolo Simulator with AVL,” Cloud Cap Technology, Hood river, OR, 18 April 2007. 18 “Aerospace Blockset 3 User’s Guide,” MATLAB Simulink, The MathWorks, Natick, Maryland, 2007. 19 Rogers, Robert M., “Applied Mathematics In Integrated Navigation Systems,” 3rd edition, AIAA Education Series, Reston, Virginia, 2007 20 Stevens, Brian L. and Lewis, Frank L., “Aircraft Control and Simulation,” 2nd edition, John Wiley & Sons, Hoboken, NJ, 2003. 21 Phillips, W. E., “Mechanics of Flight” Wiley, Hoboken, New Jersey, 2004. 22 McFarland, Richard E., “A Standard Kinematic Model for Flight Simulation at NASA-Ames,” NASA CR-2497, January, 1975. 23 Phillips, W. E., Hailey, C. E., “Review of Attitude Representations Used for Aircraft Kinematics,” AIAA Aerospace Sciences Meeting and Exhibit, Journal of Aircraft, Reno, Nevada, 2000. 24 Pamadi, Bandu N., “Performance, Stability, Dynamics, And Control of Airplanes,” 2nd edition, AIAA Education Series, Reston, Virginia, 2004 25 Cooke, J. M., Zyda, M. J., Pratt, D. R., and McGhee, R. B., “NPSNET: Flight Simulation Dynamic Modeling Using Quaternions,” Presence, Vol. 1, No. 4, pp. 404-420, January 25, 1994. 26 Icker, B. P., “A New Method for Performing Digital Control System Attitude Computations Using Quaternions,” AIAA Guidance, Control, and flight Dynamics Conference, Pasadena, California, August 12-14, 1968. 27 Hepperle, Martin, “JavaProp – Design and Analysis of Propellers,” http://www.mh-aerotools.de/airfoils/javaprop.htm 28 NIMA TR8350.2, 3rd edition, “Department of Defense World Geodetic System 1984, Its Definition and Relationship with Local Geodetic Systems,” 3 January 2000. 29 NASA-TM-X74335, “U.S. Standard Atmosphere, 1976,” National Oceanic and Atmospheric Administration, National Aeronautics and Space Administration, United States Air Force, Washington, D.C., October 1976. 30 Jategaonkar, Ravindra, “Flight Vehicle System Identification: A Time Domain Methodology,” 1st edition, AIAA Education Series, Reston, Virginia, 2006. 31 Keshmiri, Shahriar, Leong, Edmond, Jager, Rylan, and Hale, Richard, “Modeling and Simulation of the Yak-54 Scaled Unmanned Aerial Vehicle Using Parameter and System Identification,” AIAA Atmospheric Flight Mechanics Conference and Exhibit, Honolulu, Hawaii, 18-21, Aug, 2008. 32 “Mile High Wings,” http://milehighwings.com/. 33 Perry, A. R., “The Flight Gear Flight Simulator,” USENIX Annual Technical Conference, UseLinux SIG sessions, Boston, MA, June-July 2004. 34 Dimock, G. A., Deters, R. W., and Selig, M. S., “Icing Scenarios with the Icing Encounter Flight Simulator,” the AIAA 41st Aerospace Sciences Meeting and Exhibit, AIAA-2003-23, Reno, NV, Jan. 2003.

26 American Institute of Aeronautics and Astronautics 020208