 # Appendix B. Nonlinear Six-Degree-of-Freedom Aircraft Simulation

## Supplement for Flight Dynamics, a book by Robert F. Stengel FLIGHT.m provides a six-degree-of-freedom simulation of an aircraft, as well as trimming calculations and the generation of a linearized model at any flight condition chosen by the user. It is a tutorial program, heavily commented to aid interpretation. Changes to aircraft control histories, initial conditions, flag settings, and other program control actions are made by changing the numbers contained in the code; there is no separate user interface. The code has been designed for simplicity and clarity not for speed of execution, leaving a challenge to the interested reader to find ways to make the program run faster. Numerous additions could be made to the code, including implementation of feedback control logic, simulation of random turbulence or microburst wind shear, and interfaces for real-time execution. No explicit or implicit warranties are made regarding the accuracy or correctness of the computer code.

The mathematical models given below are motivated by full-scale aircraft; however, they should not be construed as accurate representations of any specific aircraft, as considerable liberties have been taken in the derivation of the models. Specifically, Mach effects are estimated using the methods of Section 2.4, rotary derivatives are calculated with the equations of Section 3.4, power effects are limited to a simple model of thrust (with no modeling of associated aerodynamic effects), and moments and products of inertia are estimated using simplified mass distributions.

The aircraft dynamic models are presented as MATLAB functions that are called by the program FLIGHT.m. FLIGHT.m illustrates the essential features of flight analysis and simulation using a modular framework. The program provides the option of computing rotation matrices using either Euler angles or quaternions. Aircraft characteristics are described in the AeroModel.m function. Any aircraft other than the examples given here can be simulated by a user-defined AeroModel.m function.

### B.1 Main Program for Analysis and Simulation (FLIGHT.m)

FLIGHT.m is the executive file (or script) that calls program functions. Initial conditions are defined here, the three primary features (trim, linearization, and simulation) are enabled, and output is generated. Initial perturbations to trim state and control allow transient effects to be simulated. As shown, trimming for steady, level flight is accomplished by first defining a cost function, J, that contains elements of the state rate, then minimizing the cost using the Downhill Simplex (Nelder-Mead) algorithm contained in fminsearch. The longitudinal trimming parameters are stabilator angle, throttle setting, and pitch angle. The linear model is generated by numjac, a numerical evaluation of the Jacobian matrices associated with the equations of motion. The linear model is saved to disk files in the variables Fmodel.mat and Gmodel.mat. MATLAB's ode23, ode45, or ode15s integrate the equations of motion to produce the state history. The state history is displayed in time plots, with angles converted from the radians used in calculation to degrees.

The user can readily change the units of plotted quantities or add additional plots through minor modifications to the code. Any result (e.g., numerical values of the state history) can be displayed in the MATLAB Command Window simply by removing the semi-colon at the end of the line of code. The MODEL flag selects either the low-angle-of-attack, Mach-dependent model for BizJet A, the high-angle-of-attack, low-subsonic model for Bizjet B, or a user-defined specification of inertial, geometric, aerodynamic, and thrust properties. The QUAT flag selects either Euler-angle or quaternion calculation of the rotation (direction cosine) matrix.

### B.2 Low-Angle-of-Attack, Mach-Dependent Model for BizJet A (AeroModelMach.m)

This version of AeroModel.m uses aerodynamic and dimensional data contained in Ref. B-1 with estimates of inertial properties of the generic business jet. All of the constants defined in this function are specific to BizJet A, and they would be different for other airplane types. Details of the configuration, such as sweep and aspect ratio of the wing and tail, are used in the estimates of Mach effects. Wind tunnel results were reviewed for all test conditions and control settings, then reduced to simple linear and quadratic coefficients that describe significant aerodynamics for low angles of attack and sideslip (within about +/-8 degrees of zero). Increments for landing gear, symmetric spoilers, and flaps are added, as appropriate. Estimates of Mach effects are based on the Prandtl factor (Section 2.4) or the modified Helmbold equation (eq. 2.4-9). Thrust available is assumed to be proportional to the air density ratio raised to a power that fits both static thrust at sea level and the thrust required for flight at maximum Mach number (0.81) at an altitude of 13,720 m (45,000 ft).

### B.3 High-Angle-of-Attack, Low-Subsonic Model for Bizjet B (AeroModelAlpha.m)

The model for BizJet B is derived using handbook methods for estimating geometric, inertial, and aerodynamic characteristics (Sec. 2.4, 2.5, 3.4, and 3.5). The model is first built using GeoMassAero.m, which saves three files describing the airplane: InerGeo.mat, DataTable.mat, and RotCont.mat. GeoMassAero.m uses a "lumped-parameter" approach to define the model; the dimensions, locations, and shapes of components (e.g., wing, horizontal tail, vertical tail, ventral fins, fuselage, landing gear, nacelles, engines, fuel, and various systems) are used to estimate mass, moments and products of inertia, center of mass, and non-dimensional aerodynamic coefficients.

AeroModelAlpha.m loads the .mat files for use in FLIGHT.m:

The angle-of-attack range extends from -10 to 90 deg based on conventional low-alpha and Newtonian high-alpha estimates, with a quadratic interpolation between the two regimes. Linear interpolation is used between tabulated data points. The controls are elevator, ailerons, rudder, and thrust. No Mach, landing gear, spoiler, or flap effects are considered.

### B.4 User-Defined Model for Advanced Jet Trainer (AeroModelUser.m)

This version of AeroModel.m is a much simplified Mach-dependent function that is based on a component buildup, as described in the previous section. Incompressible-flow aerodynamic coefficients are tabulated, with Mach-dependency specified by the subsonic Prandtl factor. The model could easily be extended to supersonic flight by an extension of the Mach multiplier. This version serves as a template for modeling other aircraft types.

### B.5 Supporting Functions

#### B.5.1 Equations of Motion

The equations of motion for functions EoM.m and EoMQ.m are written using the flat-earth assumption, as described in Section 3.2 of Flight Dynamics. Linear and translational rates are expressed in body axes, linear position is expressed in earth-relative axes, and Euler angles (EoM.m) or quaternions (EoMQ.m) describe the orientation of the body frame with respect to the earth frame. Input and output orientations are expressed as Euler angles. If quaternion calculation of the rotation matrix is selected, transformations to and from quaternions are internal to the program. The lift and drag coefficients produced by AeroModel.m are transformed to body-axis coefficients; the remaining coefficients and thrust produced by the function are expressed in body axes. Effects of the wind are calculated by WindField.m prior to calling AeroModel.m; hence, they are accounted for in the aerodynamic calculations. With the Euler-angle option, an ad hoc limit on the cosine of the pitch angle is imposed to prevent singular calculations in near-vertical flight. This expediency introduces a small error when the pitch angle is near +/-90 deg. The quaternion option is non-singular throughout the flight envelope and does not require the approximation.

The function event.m specifies a stopping condition that terminates the simulation before the final time if the altitude goes below zero.

#### B.5.2 Cost Function for Aerodynamic Trim (TrimCost.m)

Trim control settings are calculated by minimizing the quadratic function of longitudinal accelerations (i.e., rates of change of axial velocity, normal velocity, and pitch rate) contained in TrimCost.m. As shown, equal weight is given to each component in the cost function, J, because the defining matrix, R, is an identity matrix. The parameters that are adjusted to minimize J are the stabilator angle, throttle setting, and pitch angle (which equals the angle of attack in steady, level flight). TrimCost.m calls EoM.m to generate the needed accelerations. x(1) and x(3) are varied to maintain constant airspeed with varying angle of attack.

#### B.5.3 Rotation Matrices (DCM.m and RMQ.m)

DCM.m computes the direction-cosine (or rotation) matrix as a function of Euler angles. RMQ.m computes the rotation matrix with quaternions. The rotation matrices transform vectors from the earth-relative frame of reference to the body-axis frame.

#### B.5.4 Linear System Matrices (LinModel.m)

The function LinModel.m calls function EoM.m to generate f(x,u) and appends seven dummy elements to account for the controls. Thus, xdotj is a function of [xT uT], and numjac calculates the Jacobian matrix evaluated at the nominal values of state and control. The stability matrix, F, is contained in the upper-left (12 x 12) partition, and the control-effect matrix, G, is contained in the upper-right (12 x 7) partition of the Jacobian.

#### B.5.5 Wind Field (WindField.m)

WindField.m produces a three-component wind vector as a function of altitude, with linear interpolation between tabulated points. The first point is tabulated at negative altitude to assure that no computational problems occur at zero altitude. The wind vector is rotated to body axes for application in FLIGHT.m and EoM.m.

#### B.6.6 Atmospheric State (Atmos.m)

1976 U.S. Standard Atmosphere air density, air pressure, air temperature, and sound speed are generated as functions of altitude by Atmos.m. The first point is tabulated at negative altitude (-1000 m) to assure that no computational problems occur at zero altitude.

#### Reference

B-1 Soderman, P. T., and Aiken, T. N., "Full-Scale Wind-Tunnel Tests of a Small Unpowered Jet Aircraft with a T-Tail," NASA TN D-6573, Washington, DC, Nov. 1971.  http://www.stengel.mycpanel.princeton.edu/FDcodeB.html
Last updated on November 21, 2018.