The code has been designed for simplicity and clarity, leaving a challenge to find ways to make the program run faster or more efficiently. Possible additions include simulation of "frozen" random turbulence or microburst wind shear, flight in aggressive maneuvering, and reproduction of accident scenarios. The early phase of rocket launch trajectories or the late phase of atmospheric entry can be simulated by declaring the proper initial conditions and vehicle model.
V/STOL aircraft flight simulation can be accommodated within the program framework. Aircraft mass is assumed to be constant during the simulation. Modeling mass variation with thrust, atmospheric state, and control setting requires the addition of at least one element to the state variable. With additional programming, the suite of control effectors can be expanded, and closed-loop control can be incorporated. The flat-Earth dynamic equations can be replaced by round, rotating earth equations (Section 3.3 of Flight Dynamics) for 6-DOF launch, orbit, and reentry simulation.
The MATLAB code is written in chapter format as follows:
TABLE OF CONTENTS
The program provides the option of computing rotation (direction cosine) matrices using either Euler angles or quaternions. FLIGHTv2.m and its associated functions could be translated to similar open-source languages, such as Python, Scilab, or GNU Octave. No explicit or implicit warranties are made regarding the accuracy or correctness of the computer code. The code is subject to continuing change, and user suggestions for improving the code are welcome.
The twelfth-order mathematical models are motivated by actual aircraft; however, they should not be construed as accurate representations of any particular aircraft. Aerodynamic coefficients are either angle-of-attack-dependent or Mach-number-dependent functions. Inertial properties (total mass and moments of inertia) are constant during a simulation run. The aircraft dynamic model is specified by ACtype (aircraft type, see FLIGHTv2.m) and MODEL (angle-of-attack-dependent or Mach-dependent math model). The user determines which aircraft to simulate by setting these flags in the program. Aircraft characteristics are described in a function such as BusinJetA.m or BusinJetM.m. ModelBuild.m is an adjunct program for creating mathematical models using handbook methods, as described in Section B.4. Other aircraft types can be simulated in user-defined m-files.
FLIGHTv2.m is the executive file (or script) that calls program functions. Initial conditions are defined, the three primary features (trim, linearization, and simulation) are enabled, and computed output is generated. The nonlinear equations of motion take the vector form, dx/dt = f[(x(t), u(t),w(t)], x(0) given; x, u, and w are the state, control, and disturbance vectors. Primary features are described as follows:
The user can change the units of plotted quantities or add additional plots. 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 ACtype flag selects either a low-angle-of-attack, Mach-dependent model or a high-angle-of-attack, low-subsonic model. The RM flag selects either Euler-angle or quaternion calculation of the rotation matrix.
The aircraft model must be determined by running a script such as BusinessJetA.m, Business JetM.m, or ModelBuildJT.m prior to running FLIGHTv2.m (Sections B.2.1, B.2.2, B.4.2).
The inertial, geometric, and aerodynamic characteristics of the aircraft type are specified by the parameter ACtype, which selects a function such as BusinJetA.m or BusinJetM.m to download the tabulated parameters. The specified parameters are saved in a mat-file that is loaded by the run-time function AlphaModel.m or MachModel.m, as selected by MODEL. The run-time function is called by the equations-of-motion function EoMv2.m or EomQv2.m for use in trimming and simulating aircraft motions and for calculating a linearized model at the trim condition.
Trim conditions for the angle-of-attack-dependent and Mach-dependent models of an aircraft type are approximately the same at one flight condition. For the business jet example, this condition is an altitude of 3,050 m and Mach number of 0.3. For all other flight conditions, trim values, LTI models, and simulated transient responses are different. Large-amplitude responses for Mach-independent aerodynamic coefficients are computed using AlphaModel.m. Small-amplitude responses for aerodynamic coefficients that vary with Mach number are computed using MachModel.m. If aerodynamic data for concurrent angle of attack and Mach number variations are available, an alternative function that follows the AlphaModel.m and MachModel.m structures could be scripted. (The scalar arguments for table interpolation in each function would be replaced by vector arguments of angle of attack and Mach number.)
The BusinJetA.m script expresses the inertial, geometric, aerodynamic, and thrust properties of the aircraft in incompressible flow (approximating M < 0.5). The angle-of-attack ranges from -10 to +90 deg, and the sideslip-angle range is about +/-8 degrees. Low-angle-of-attack static coefficients up to stall (angle of attack = 14 deg) are based on the wind tunnel data given in [B-1]. High-angle-of-attack coefficients (i.e., beyond 35 deg) are Newtonian estimates based on the aircraft planform. The low- and high-angle data are connected by cubic curve fits to the values and slopes of the low- and high-alpha coefficients at the matching conditions. Control surface effects are linear in their respective deflections and are tapered to zero at 90-deg angle of attack. Static lateral-directional effects are linear in sideslip angle. Dynamic coefficients are assumed to be functions of the static coefficients, as presented in Section 3.4 of the book. Thrust available is 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) and altitude of 13,720 m (45,000 ft). It is independent of angle of attack. BusinJetA.m is called prior to FLIGHT.m initialization. Tabulated parameters are saved in AlphaModelData.mat, which is loaded by AlphaModel.m.
The BusinJetM.m script uses aerodynamic and dimensional data contained in [B-1] with estimates of inertial properties of a generic business jet. Mass, inertial, geometric, aerodynamic, and thrust properties are tabulated in BusinJetM.m. Details of the configuration are extrapolations of the subsonic wind-tunnel data based on sweep and aspect ratio of the wing and tail. Estimates of Mach effects are based on the Prandtl factor or the modified Helmbold equation (Section 2.4 of Flight Dynamics). Thrust is modeled as in the alpha-dependent case. BusinJetM.m is called during FLIGHT.m initialization. Tabulated parameters are saved in the mat-file MachModelData.mat, which is loaded by MachModel.m.
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 AlphaModel.m or MachModel.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 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.
Trim control settings are calculated by minimizing a quadratic function of longitudinal accelerations (i.e., rates of change of axial velocity, normal velocity, and pitch rate) contained in TrimCost.m. Equal weight is given to each component in the cost function, J; the defining matrix, R, is an identity matrix. The stabilator angle, throttle setting, and pitch angle (which equals the angle of attack in steady, level flight) are adjusted to minimize J.
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.
LinModelv2.m calls function EoMv2.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.
NatFreq presents complex-valued eigenvalues, natural frequencies, damping ratios, and time constants that identify oscillations, convergence, and divergence in the system's modes of motion.
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.
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.
The ModelBuild.m script uses handbook methods to estimate inertial and aerodynamic properties of a subject aircraft from available aircraft specifications. The specific constants and equations contained in the script depend on the configuration being modeled; for example, tailed and tailless aircraft have different control surfaces, thrust models are different for reciprocating and jet engines, and thrust offset effects are different for single- and multi-engine models. The script first computes AlphaModelData.mat and then MachModelData.mat, expanding on and combining the calculations of Sections B.2.1 and B.2.2. ModelBuild.m is run prior to executing FLIGHTv2.m. A prototype is presented in ModelBuildJT.m.
TABLE OF CONTENTS
The prototype ModelBuild.m script, presents a model for a generic jet trainer with supersonic capability. The script ModelBuildJT.m generates both alpha-dependent and Mach-dependent models, with output stored as AlphaModelData.mat and MachModelData.mat. The incompressible-flow model, identified as ACType = JetTrainA, is the first to be computed, and its output provides input for AlphaModel.m; hence, MODEL = Alph. The compressible flow model is specified by JetTrainM and Mach, with Mach-dependent coefficients extrapolated from M = 0 via Prandtl-Glauert-Goethert-Ackeret and Helmbold functions contained in MachSimFcnJT.m.
The generic jet trainer is similar to the USN Super Hornet, although it is smaller and has a single engine. Aerodynamic and inertial properties of the trainer are estimated from three-view drawings, which are scaled according to the wingspan ratio of the two aircraft. Scaling allows areas, chord and fuselage lengths, taper ratios, and sweep angles to be estimated. Masses and moments of inertia of components are estimated using conventional structural analysis of simple geometric shapes and the parallel-axis theorem, as in [R-1]. For this conceptual design, static and dynamic aerodynamic coefficients are estimated with handbook methods, following Sections 2.4 and 3.4 of Flight Dynamics, Second Edition and [F-1]. Thrust for the trainer model is based on estimated characteristics of a current powerplant.
The aerodynamic coefficients are tabulated separately as functions of angle of attack (-10 deg to 90 deg) and Mach number (0 to 1.5). The alpha-dependent data are stored in AlphaModelData.mat, which serves as input to the simulation function AlphaModel.m. The Mach-dependent data are stored in MachModelData.mat, which serves as input to the simulation function MachModel.m.
This function provides a smooth approximation in a region between scalar aerodynamic coefficients that are tabulated for low and high values of the independent variable. Angle-of-attack-dependent coefficients may be tabulated for partially attached or fully separated upper-surface flows at low and high angles. For the alpha-dependent model, the cubic fit provides estimates in the intermediate range from stall to post-stall that match the values and slopes of the coefficients at the interior points. For the Mach-dependent case, the cubic fit provides estimates in the transonic regime, matching the values and slopes of the coefficients at the interior points.
Flight Dynamics, Second Edition.
http://www.stengel.mycpanel.princeton.edu/FDcodeB.html